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

solargus

07-15-2014, 11:25 AM

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, 03: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, 05:41 AM

Thank you very much, Allan!

solargus

07-18-2014, 08: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, 05: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, 05: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, 08:23 PM

sounds good to me;

thanks for the information!

allan

Powered by vBulletin® Version 4.2.3 Copyright © 2019 vBulletin Solutions, Inc. All rights reserved.