Interpreting frequency using pwelch function

Hello,

I'm just looking for some advice in interpreting the frequency axis from a Power Spectral Density estimate by Welch's method. I'm very much new to this function. I have a time series with data roughly 150 years in length, and reported in timesteps of two week intervals. Working from the example given in the help function, I used

pwelch(y,[],[],[],1/26.0893,'onesided')

trying to account for leap years. My goal is to interpret the resulting power in units of cycles/year. The frequency axis from the result is by default giving me information in mHz, and I'm struggling to relate that to anything physical in the dataset.

If I'm doing something wrong, I'd appreciate a pointer in the right direction. Thanks.

Answer 1

[Wayne King]

If you sample every two weeks, your frequencies range from 0 cycles/(2 weeks) up to 1 cycle/(4 weeks) your Nyquist frequency.

You can simply scale your frequency axis to display in cycles/year

   x = randn(15600,1);
   % Fs - sampling frequency is 1 sample/2 weeks (1/2)
   [Pxx,Fxx] = pwelch(x,[],[],[],1/2);
   Fxx = Fxx.*52;
   plot(Fxx,10*log10(Pxx)); xlabel('Cycles/Year');
   grid on; ylabel('dB (Power)');

[Wayne King]

yes, I realized that after that is 15,600 samples

[Dr. Seis]

Wouldn't it just be 26.0893(samples/year)*150(years) = 3913.4 samples?

[Wayne King]

Oh, yea, sorry for some reason I gave him two samples/week instead of 1 sample every two weeks :) Should have been 150*52/2

Answer 2

[Dr. Seis]

1/26.0893 years per sample equals roughly 1208771.41 seconds per sample, which is roughly equivalent to an Fs of 8.27286276e-07 samples per second. Matlab assumes that the units of Fs are samples per second... so if you were to give an Fs with units of samples per year (26.0893 in your case, not 1/26.0893), then labeling the axis in terms of Hz (which it would do automatically) would be nonsensical.

Giving "pwelch" that Fs=8.27e-07 would at least allow the frequency axis to be defined properly (in terms of Hz) as it is automatically displayed, but you would have to do some mental math to convert cycles/sec to cycles/year. Alternatively, you can set Fs equal to 26.0893 and then just ignore the "Hz" label or manually change the label to "cycles per year".

Answer 3

[Honglei Chen]

Hi Chris,

pwelch assumes that the sampling frequency is specified in terms of seconds. When you specify 1/26.0893, you are saying that the sampling inteval is 1/26.0893 years. I would suggest you to change it to 26.0893 which means that you sample about 26 times a year. If you use that number, in the display, you can assume that 1 Hz = 1 cycle/year. Of course, you can also plot it yourself, such as:

[Px,f] = pwelch(rand(1000,1),[],[],[],26.0893,'onesided')
plot(f,abs(Px))
xlabel('Cycles/year')

Answer 4

[Chris]

Thanks everyone,

So if I understand correctly, there seems to be agreement that for a dataset with time of length t and corresponding data y,

pwelch(y,[],[],[],26.0893,'onesided')

would tell me what I want? (even if the axis is not in Hz, this seems fine)

[Dr. Seis]

Yes, it should... but you should also consider manually changing the x-axis label as I suggested and Honglei demonstrated in his example (just to keep units straight).