Results 1 to 7 of 7

Thread: Power Spectrum On 1-D Time Series Data

  1. #1
    Junior Member
    Join Date
    Jul 2010
    Posts
    14

    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!
    Last edited by solargus; 07-15-2014 at 11:26 AM. Reason: Oops. left out a step.

  2. #2
    Member
    Join Date
    Sep 2005
    Posts
    41
    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

  3. #3
    Junior Member
    Join Date
    Jul 2010
    Posts
    14
    Thank you very much, Allan!

  4. #4
    Junior Member
    Join Date
    Jul 2010
    Posts
    14
    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!

  5. #5
    Member
    Join Date
    Sep 2005
    Posts
    41
    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

  6. #6
    Junior Member
    Join Date
    Jul 2010
    Posts
    14
    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

  7. #7
    Member
    Join Date
    Sep 2005
    Posts
    41
    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
  •