Converting spectrum.pwelch//psd to equivalent pwelch call

5 次查看(过去 30 天)
I have successfully converted a number of my code blocks from the older method to pwelch. Each of those explicitely called out a window type and segment length.
I can not convert the simplest case because I do not understand what the defaults for spectrum.pwelch & psd are. I do not think I am specifying the correct segment length or window to match the defaults in the original code. I have tried a number of window and segment lengths to try and get close to original psdEst spectrum ... but not even close thus far.
Original Code
h=spectrum.welch;
psdEst = psd(h, timeSeries, 'FS',44100); timeSeries is a 20000 sample sinewave type signal
Replacement Code
psdEst = pwelch(timeSeries,[],[],[],44100);
  1 个评论
Mathieu NOE
Mathieu NOE 2024-9-5
Each of those explicitely called out a window type and segment length.
I don't see that neither in the original code or replacement code that you posted . The old and new pwelch functions may have different default settings so we have to use more parameters to make sure you are still having matching results

请先登录,再进行评论。

回答(1 个)

William Rose
William Rose 2024-9-5
The defaults for spectrum.welch are listed here.
The defaults for pwelch are listed toward the bottom of this page.
As you can see, the use of spectrum.pwelch is not recommended, and the use of pwelch is recommended.
  1 个评论
William Rose
William Rose 2024-9-5
I like your choice of psdEst for the output, because it is a good reminder that when we compute a spectrum, we are computing an estimate of the true but unknown spectrum.
You can observe default values of spectrum.pwelch by omitting the semicolon when defining h:
h=spectrum.welch
h =
EstimationMethod: 'Welch' SegmentLength: 64 OverlapPercent: 50 WindowName: 'Hamming' SamplingFlag: 'symmetric'
Create a signal: two sinusoids plus noise.
fs=44100; % sampling rate (Hz)
t=0:1/fs:0.1; % time vector (0.1 sec)
f1=7789; f2=10103; % frequencies (Hz)
x=cos(2*pi*f1*t)+sin(2*pi*f2*t)+randn(size(t));
Find spectrum by two methods:
psdEst1=psd(h,x,'FS',fs);
[psdEst2,f2]=pwelch(x,64,32,[],fs);
Plot the two spectra:
figure;
subplot(211), plot(psdEst1.Frequencies,psdEst1.Data,'-b.')
grid on; ylabel('PSD'); title('spectrum.welch')
subplot(212), plot(f2,psdEst2,'-r.')
grid on; xlabel('Frequency (Hz)'); ylabel('PSD'); title('pwelch')
The spectra look the same. I chose the values for pwelch() arguments window and overlap, based on the values displayed when I defined h.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Spectral Estimation 的更多信息

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by