# Thread: Power Spectrum On 1-D Time Series Data

1. ## Power Spectrum On 1-D Time Series Data

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!  Reply With Quote

2. 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  Reply With Quote

3. Thank you very much, Allan!  Reply With Quote

4. 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!  Reply With Quote

5. 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  Reply With Quote

6. 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  Reply With Quote

7. sounds good to me;
thanks for the information!
allan  Reply With Quote

#### Tags for this Thread

fft, power spectrum #### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•