
Power Spectrum On 1D Time Series Data
I'd like to generate a power spectrum plot for a time series data set using the Foundation version of PVWave on a Unix platform.
So if I have a data set x, which is a 1D, doubleprecision array of 100000 elements, I'm attempting to create the power spectrum with:
n = n_elements(x)
y = fft(x, 1)
t=1800 ;the time interval between samples, in seconds
fs = 1/t ;sample frequency
y = fft(x, 1)
power = y * conj(y) / n
frequency = double(n)
for i=0, n1 do frequency(i)=i
frequency = frequency*(fs/n)
plot, frequency, power
Am I close, or is there a better way to do this? Also, how can I scale the yaxis as log(10)?
Thanks!
Last edited by solargus; 07152014 at 12:26 PM.
Reason: Oops. left out a step.

Hi Solargus;
You were very close.
Below is a working example for you.
Please note the following changes from your original code:
(a) t is a double so integer division is not used to get fs;
(b) the fft is shifted so the negative frequencies show up on the plot;
(c) the frequencies are computed to match (b);
(d) log scaling is used on the yaxis;
Thanks;
allan
n = 100000
t = 1800d
fs = 1 / t
u = t * dindgen(n)
x = cos(2*!pi*fs/4*u) + sin(2*!pi*fs/8*u) + 1d6*randomn(s,n)
y = fft( x, 1 )
y = shift( y, n/2 )
power = y * conj(y) / n
frequency = dindgen(n)  (n/2)
frequency = (fs/n) * frequency
plot, frequency, power, ytype=1

Thank you very much, Allan!

There was one more thing I was hoping to accomplish and that is to have Wave print out the dominant frequencies. For that, I thought of using the EXTREMA function, but as it turns out, you cannot exercise that function on the complex data type. Any suggestions?
Thanks again!

Hi Solargus;
EXTREMA would not work anyway since it gives local extrema, not a sorted list.
Continuing with the previous example, you can instead do
pm, frequency((sort(double(power)))(0:3))
which prints the four dominant frequencies.
allan

Hi Allan,
With most of the power being at the Nyquist in my particular signal, what I am getting from the SORT command is a list of frequencies at and near that frequency. This is the case even when you expand the range from 0:3 to hundreds. Perhaps the better way to do this is to calculate the mean then perform an NSigma filter.
Solargus

sounds good to me;
thanks for the information!
allan
Tags for this Thread
Posting Permissions
 You may not post new threads
 You may not post replies
 You may not post attachments
 You may not edit your posts

Forum Rules