Reproducing 'pspectrum' using 'periodogram'

Hello,

I am trying to understand how pspectrum works. After reading the documentation, I thought I understood how it uses a kaiser window and developed the below code to try to reproduce results. But unfortunately I still seem to be missing something.

Also, curious why the magnitudes of power are also different? By like a factor of ~500. Any insight into this would also be helpful.

Would greatly appreciate to see an exact code of how to reproduce and any other pre-processing that is needed to match the two methods. Thank you in advance!

Cheers,

Edwin

% Method 1, Power Spectrum using 'pspectrum'.
[p,f] = pspectrum(x,fs,'Leakage',0.5);

% Method 2, Power Spectrum using 'periodogram'
lk = 0.5;
M = length(x);
g = kaiser(M,40*(1-lk));
[p,f] = periodogram(x,g,length(x),fs);

[Edwin]

After looking into the code of 'pspectrum' I believe I see now how it's calculating power via Matlab signal analyzer.

Documentation for signal analyzer says:

1. Divide the signal into equal-length segments. The segments must be short enough that the frequency content of the signal does not change appreciably within a segment. The segments may or may not overlap.
2. Window each segment and compute its spectrum to get the short-time Fourier transform.
3. Display segment-by-segment the power of each spectrum in decibels. Depict the magnitudes side-by-side as an image with magnitude-dependent colormap.

Which now this makes sense seeing the other options for 'pspectrum' such as 'NumPowerBins' and 'OverlapPercent'. So, thus, it seems periodogram cannot be used to replicate pspectrum but need to use a STFT method. Other useful info can be found here.

Cheers

Answer 1

[Nayan]

As I understand it you would like to understand the working of "pspectrum(x)". I would suggest you to note down the following points for a better understanding.

To compute signal spectra, pspectrum finds a compromise between the spectral resolution achievable with the entire length of the signal and the performance limitations that result from computing large FFTs:

• If possible, the function computes a single modified periodogram of the whole signal using a Kaiser window.
• If it is not possible to compute a single modified periodogram in a reasonable amount of time, the function computes a Welch periodogram: It divides the signal into overlapping segments, windows each segment using a Kaiser window, and averages the periodograms of the segments.

Hence, there are two different techniques in which pspectrum is evaluated depending on the situation at hand. I would suggest you go through the following MATLAB documentation link for a detailed understanding https://www.mathworks.com/help/releases/R2023a/signal/ref/pspectrum.html?s_tid=doc_ta#d124e149267

I would suggest you to share the exact data over which you are evaluating the power-spectrum using "pspectrum(x)" and "periodogram(x,window,f,fs)" for comparision(since the method used for evaluating the power spectrom very much depends on the data)

You can go through the following link to know more about "periodogram(x,window,f,fs)" and its working.

https://www.mathworks.com/help/releases/R2023a/signal/ref/periodogram.html?s_tid=doc_ta#btuc2ow-1

Also find the code snippet that produces the same power spectrum with both methods.

% Generate a test signal with sampling frequency 100 Hz and length 1000 

Fs = 100; % Sampling frequency (Hz) 

t = (0:999)/Fs; % Time vector (sec) 

x = cos(2*pi*10*t) + 0.5*sin(2*pi*20*t) + randn(1,1000); % Test signal 

% Compute the PSD using the periodogram method 

[pxx,f] = periodogram(x,[],[],Fs);  

% Visualize the PSD 

plot(f,10*log10(pxx)) 

xlabel('Frequency (Hz)') 

ylabel('Power/frequency (dB/Hz)') 

title('Periodogram PSD estimate') 

% Zoom in to show only frequencies up to 50 Hz 

xlim([0 50]) 

 


% Generate a test signal with sampling frequency 100 Hz and length 1000 

Fs = 100; % Sampling frequency (Hz) 

t = (0:999)/Fs; % Time vector (sec) 

x = cos(2*pi*10*t) + 0.5*sin(2*pi*20*t) + randn(1,1000); % Test signal 

 
 

% Compute the PSD using the "pspectrum" function with the default settings 

[pxx,f] = pspectrum(x,Fs); 

 
 

% Visualize the PSD 

plot(f,10*log10(pxx)) 

xlabel('Frequency (Hz)') 

ylabel('Power/frequency (dB/Hz)') 

title('PSD estimate using pspectrum') 

 
 

% Zoom in to show only frequencies up to 50 Hz 

xlim([0 50]) 

Hope this helps!

[Karl Warnick]

I think this answer is a bit misleading, as pspectrum and periodogram are scaled differently. pspectrum gives the total power in each frequency bin, whereas periodogram defaults to spectrumtype = 'psd' and the output is scaled to power per Hz, rather than per bin.