Converting spectrum.pwelch//psd to equivalent pwelch call

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);

[Mathieu NOE]

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

Answer 1

[William Rose]

@Ron Fowler,

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.

[William Rose]

@Ron Fowler,

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.