Time-Series Analysis - MATLAB Recipes for Earth Sciences - page 161

where the function nextpow2 computes the next power of two closest to
the length of the time series x(t) . h is code allows the sampling frequency
to be modii ed and the dif erences in the results to be explored. We
can now compare the results with those obtained using the function
periodogram(x,window,nfft,fs) .
[Pxx,f] = periodogram(x,[],1024,1);
h is function allows the windowing of the signals with various window
shapes to overcome spectral leakage. However, we use the default rectangular
window by choosing an empty vector [] for window to compare the results
with the above experiment. h e power spectrum Pxx is computed using an
FFT of length nfft=1024 , which is the next power of two closest to the length
of the series x(t) and which is padded with zeros to make up the number
of data points to the value of nfft . A sampling frequency fs of one is used
within the function in order to obtain the correct frequency scaling for the
f -axis. We display the results by typing
plot(f,Pxx), grid
xlabel('Frequency')
ylabel('Power')
title('Auto-Spectrum')
or alternatively
plot(1./f,Pxx), axis([0 100 0 1000]), grid
xlabel('Period')
ylabel('Power')
title('Auto-Spectrum')
h e graphical output is almost identical to our Blackman-Tukey plot and
again shows that there are three signii cant peaks at the positions of the
original frequencies (or periods) of the three sine waves. h e same procedure
can also be applied to the noisy data:
[Pxx,f] = periodogram(xn,[],1024,1);
plot(f,Pxx), grid
xlabel('Frequency')
ylabel('Power')
title('Auto-Spectrum')
Let us now increase the noise level by introducing Gaussian noise with a
mean of zero and a standard deviation of i ve.
rng(0)
n = 5 * randn(size(x));
xn = x + n;