Interpreting frequency using pwelch function

30 次查看（过去 30 天）
Chris
Chris 2012-5-10
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.

请先登录，再回答此问题。

回答（4 个）

Wayne King
Wayne King 2012-5-10
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)');
  3 个评论
显示 1更早的评论
Dr. Seis
Dr. Seis 2012-5-10
Wouldn't it just be 26.0893(samples/year)*150(years) = 3913.4 samples?
Wayne King
Wayne King 2012-5-10
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

请先登录，再进行评论。

Dr. Seis
Dr. Seis 2012-5-10
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".
Honglei Chen
Honglei Chen 2012-5-10
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')
Chris
Chris 2012-5-10
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)
  1 个评论
Dr. Seis
Dr. Seis 2012-5-10
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).

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by