Geoscience Reference
In-Depth Information
We now drop the redundant part in the power spectrum and use only the
i rst
(1+nfft/2)
points. We also multiply the power spectral density by two
to keep the same energy as in the symmetric spectrum, except for the i rst
data point.
Pxx = [Pxx2(1); 2*Pxx2(2:512)];
h e corresponding frequency axis runs from
0
to
Fs/2
in
Fs/(nfft-1)
steps,
where
Fs/2
is the Nyquist frequency. Since
Fs=1
in our example, the frequency
axis is
f = 0 : 1/(1024-1) : 1/2;
We then plot the power spectral density
Pxx
in the Nyquist frequency range
from
0
to
Fs/2
, which in our example is from
0
to
1/2
. h e Nyquist frequency
range corresponds to the i rst 512 or
nfft/2
data points. We can plot the
power spectral density over the frequency by typing
plot(f,Pxx), grid
h e graphical output shows that there are three signii cant peaks at the
positions of the original frequencies of the three sine waves (1/50, 1/15,
and 1/5). Alternatively, we can also plot the power spectral density over the
period by typing
plot(1./f,Pxx), axis([0 100 0 1000]), grid
where we observe the three periods 50, 15, and 5, as expected. Since the
values on the
x
-axis of this plot are not evenly spaced (in constrast to those
on the frequency axis), we i nd the long periods poorly resolved and a broad
peak at a period of 50 in this graphics. h e code for the power spectral
density can be rewritten to make it independent of the sampling frequency,
Fs = 1;
t = 1/Fs :1/Fs : 1000/Fs; t = t';
x = 2*sin(2*pi*t/50) + sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
nfft = 2^nextpow2(length(t));
Xxx = fft(x,nfft);
Pxx2 = abs(Xxx).^2 /Fs /length(x);
Pxx = [Pxx2(1); 2*Pxx2(2:512)];
f = 0 : Fs/(nfft-1) : Fs/2;
plot(f,Pxx), grid
axis([0 0.5 0 max(Pxx)])
