PDA

View Full Version : Power Spectrum On 1-D Time Series Data



solargus
07-15-2014, 12:25 PM
I'd like to generate a power spectrum plot for a time series data set using the Foundation version of PV-Wave on a Unix platform.

So if I have a data set x, which is a 1D, double-precision 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, n-1 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 y-axis as log(10)?

Thanks!

allan
07-16-2014, 04:36 PM
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 y-axis;
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) + 1d-6*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

solargus
07-18-2014, 06:41 AM
Thank you very much, Allan!

solargus
07-18-2014, 09:01 AM
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!

allan
07-19-2014, 06:05 PM
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

solargus
07-21-2014, 06:40 AM
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 N-Sigma filter.

Solargus

allan
07-21-2014, 09:23 PM
sounds good to me;
thanks for the information!
allan