pspectrum power spectrum or power density scaling

Hi everybody !

I have a matter with the pspectrum power scaling. In help documentation, it is noticed that "pspectrum scales the spectrum so that, if the frequency content of a signal falls exactly within a bin, its amplitude in that bin is the true average power of the signal".

Ref.: https://fr.mathworks.com/help/signal/ref/pspectrum.html#d123e133339

I may have a trouble for understanding the operation underlying this explanation. I've coded periodograms computations (via |DFT|², periodogram and pwelch implementations, not indicated below) and all my results are matching with periodogram and pwelch functions results. But not with pspectrum outputs. Please find below an adapted code for testing my observation :

Ts = 0.001; % [s]
Fs = 1/Ts; % [Hz]
Duration = 5; % [s]
Ns = floor(Duration/Ts); % [-]
NoisePower = 0.5^2; % [unit²]
v_time = [0:(Ns-1)]'*Ts;
v_signal = sqrt(NoisePower)*randn(Ns,1);
energy_viatime = Ts * sum(abs(v_signal).^2);
power_viatime = energy_viatime/((Ns-1)*Ts);

% Windowing specifications
% The window used by pspectrum is a Kaiser window with BETA parameter = 40*(1-Leakage)
% where Leakage value is [0;1] (1 = rectwin)
beta = 0; leak = 1-beta/40; % here, i chose a rectwin, i.e. leak = 1

% pspectrum call
[v_PSP_2S, v_fPSP_2S] = pspectrum(v_signal,v_time, ...
    'Leakage',leak,'TwoSided',true,'Reassign',false,'power');
[v_PSP_1S, v_fPSP_1S] = pspectrum(v_signal,v_time, ...
    'Leakage',leak,'TwoSided',false,'Reassign',false,'power');

fprintf('Relative error of power estimation [PSP-2S]: %2.1f%%\n',abs(sum(v_PSP_2S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PSP-1S]: %2.1f%%\n',abs(sum(v_PSP_1S) - power_viatime) / power_viatime * 100);
fprintf('\n');

% Report #1: sum(v_PSP_2S) > power_viatime && sum(v_PSP_1S) > power_viatime
% -> Power creation is impossible ! v_PSP_*S is not the expected power spectrum.

% Report #2: sum(v_PSP_1S) / sum(v_PSP_2S) ~ 2
% -> Seems to be related to the frequency interval of respective output

% Report #3: mean(diff(v_fPSP_1S)) / mean(diff(v_fPSP_2S)) ~ 0.5
% -> As a matter of fact, frequency discretization is not the same whereas sizes are equal. OK.
% -> Still a power spectrum or rather a power density ???

Can you confirm my understanding, please ?

'Input arguments' section about type as 'power' mentions the power spectrum (which is different from the PSD).

No reassignment scaling is explicitly ordered.

'More about' section about Spectrum computation does not answer about that aspect.

But here, it seems that I have a beginning of an answer (just like a clue ...)

% ... trying to implement a bandwidth scale factor ... 
dF_window_enbw_2S = enbw(kaiser(size(v_PSP_2S,1),beta),Fs);
sum(v_PSP_2S)*dF_window_enbw_2S % -> may be the right one if *2. Yes, but enbw is already 2-sided.
dF_window_enbw_1S = enbw(kaiser(size(v_PSP_1S,1),beta),Fs);
sum(v_PSP_1S)*dF_window_enbw % -> may be the right one ! Yes, but enbw is 2-sided. To be single-sided ?
% => pspectrum outputs a density, not a power spectrum (in the case of white noise spectrum).

enbw documentation helps to be aware about the energy carried by each frequency bin.

There's something that I don't catch ! Several days already spent on that matter ...

Thank you for your explanation.

Answer 1

[Star Strider]

Look carefully at the units, specifically for periodogram pxx, and for pwelch pxx. Compare these with pspectrum p.

I prefer pspectrum exactly for this reason.

.

[Christophe]

Hello Star Strider,

Thank you for answering my post.

Your remark is right about the units : PSD are expressed in u²/Hz, where as Power spectrum are in u² (u for general unit). This is the reason why I have also used these functions periodogram and pwelch in two manner (with the spectrumtype input - here for periodogram, here for pwelch - set as either 'psd' or 'power').

Here some computation code :

% Windowing specifications
% The window used by pspectrum is a Kaiser window with BETA parameter = 40*(1-Leakage)
% where Leakage value is [0;1] (1 = rectwin)
beta = 0; leak = 1-beta/40;
% leak = 1; beta = 40*(1-leak);

% Schuster periodogram
% -------------------------------------------------------------------------
Nf_PER = Ns; % Ns or 2^nextpow2(Ns)
dF_PER = Fs/Nf_PER;
% Windowing initialization
% v_window_PER = rectwin(Ns);
v_window_PER = kaiser(Ns,beta);
% wvtool(v_window_PER)
energy_window_PER = Ts*sum(v_window_PER.^2);
power_window_PER = energy_window_PER / (Ns*Ts);

% Computation of PSD by periodogram function
[v_PERPSD_2S,v_fPERPSD_2S] = periodogram(v_signal,v_window_PER,Nf_PER,Fs,'centered');
[v_PERPSD_1S,v_fPERPSD_1S] = periodogram(v_signal,v_window_PER,Nf_PER,Fs,'onesided');
% Computation of Power Spectrum by periodogram function
[v_PERPSp_2S,v_fPERPSp_2S] = periodogram(v_signal,v_window_PER,Nf_PER,Fs,'centered','power');
[v_PERPSp_1S,v_fPERPSp_1S] = periodogram(v_signal,v_window_PER,Nf_PER,Fs,'onesided','power');

% Welch periodogram
% -------------------------------------------------------------------------
Ms = 4096; % Segment length for pwelch estimate
Rs = floor(Ms*0.50); % Overlap length
number_of_segments = floor((Ns-Rs)/(Ms-Rs));
Nf_PWE = Ms; % Ms or 2^nextpow2(Ms)
dF_PWE = Fs/Nf_PWE;
% Windowing initialization
% v_window_PWE = rectwin(Ms);
v_window_PWE = kaiser(Ms,beta);
% wvtool(v_window_PWE)
energy_window_PWE = Ts*sum(v_window_PWE.^2);
power_window_PWE = energy_window_PWE / (Ms*Ts);

% Computation of PSD by pwelch function
[v_PWE_2S,v_fPWE_2S] = pwelch(v_signal,v_window_PWE,Rs,Nf_PWE,Fs,'centered');
[v_PWE_1S,v_fPWE_1S] = pwelch(v_signal,v_window_PWE,Rs,Nf_PWE,Fs,'onesided');
% Computation of Power Spectrum by pwelch function
[v_PWEPSp_2S,v_fPWEPSp_2S] = pwelch(v_signal,v_window_PWE,Rs,Nf_PWE,Fs,'centered','power');
[v_PWEPSp_1S,v_fPWEPSp_1S] = pwelch(v_signal,v_window_PWE,Rs,Nf_PWE,Fs,'onesided','power');

And here some comparison code :

(for checking coherence with the power amount identified in the time signal) :

fprintf('Relative error of power estimation [PER-2S]: %2.1f%%\n',abs(sum(v_PERPSp_2S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PER-1S]: %2.1f%%\n',abs(sum(v_PERPSp_1S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PER-2S]: %2.1f%%\n',abs(sum(v_PERPSp_2S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PER-1S]: %2.1f%%\n',abs(sum(v_PERPSp_1S) - power_viatime) / power_viatime * 100);
fprintf('\n');

fprintf('Relative error of power estimation [PWE-2S]: %2.1f%%\n',abs(sum(v_PWEPSp_2S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PWE-1S]: %2.1f%%\n',abs(sum(v_PWEPSp_1S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PWE-2S]: %2.1f%%\n',abs(sum(v_PWEPSp_2S) - power_viatime) / power_viatime * 100);
fprintf('Relative error of power estimation [PWE-1S]: %2.1f%%\n',abs(sum(v_PWEPSp_1S) - power_viatime) / power_viatime * 100);
fprintf('\n');

When you sum all the elements of the output vector of periodogram or pwelch, you can check that this sum equals the power observed during the time evolution (according to the Parseval's theorem).

pspectrum seems to not comply with this simple coherence check, where periodogram and pwelch does.

I performed some modifications since yesterday and my comments in my initial post would have to be modified, especially about the effect of the window on the ratio between sum(v_PSP) and power_viatime. But the observation is the same.

pspectrum is creating power (!!!) if I strictly consider the definition given by the help for p output for 'power' input argument : "p contains the power spectrum estimate of each channel of x.".

Thank you again for your care !

[Star Strider]

I cannot run that, so I cannot comment on it.

Power is generally considered to be amplitude² and I believe that is what psepectrum displays. I doubt that it is creating anything.

[Christophe]

I fully agree with you. Lavoisier said it : "Nothing is lost, nothing is created, everything is transformed".

So I am rather suspecting a PSD [u²/Hz] or Power spectrum [u²] scaling issue, that may be ambiguous in the documentation and which is not so obvious when inspecting the code of pspectrum.m or ZoomSpectrum or ZoomEstimator files of the Digital Signal Processing Toolbox.

Thanks a lot nonetheless for your conversation.

Maybe any other aware contributors ?

[Star Strider]

My pleasure.

[Paul]

@Christophe

Can you show or post a link to the definition of power spectrum?

Answer 2

[Christophe]

Hello Paul,

Yes, of course : https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density

I may understand that it's a little bit confusing, one is a density (expressed as unit²/Hz), the other is a power (unit²).

Power spectral density (PSD) :

==> It is the quantity returned by periodogram or pwelch when input argument spectrumtype is 'psd' (default).

What I call power spectrum :

==> It is the quantity returned by periodogram or pwelch when input argument spectrumtype is 'power'.

The difference, or more exactly the ratio, between these two objects is a scale factor linked to the "weigth" of the elementary bandwidth, that is to say each bin. I may have found an explanation on this website : here (DSP Stack Exchange), and especially this note as a link to a paper from academics in Switzerland.

pspectrum may be a too advanced, too integrated or too automated function, with convenient pre-settings, so powerful but less-tunable. pspectrum does not allow to specify if 'psd' or 'power' scaling is desired. It has to be understood as 'power' (in u²), and not as 'psd' (in u²/Hz), which is not so practical when you analyze wide-band spectrum (such as noises). For peaks, it is easy-to-use and ready-to-read.

Your opinion is welcome.

[Paul]

Hi Christophe,

I was happy to have found this thread because, coincidentally, I've been reading about this topic on and off for the last few weeks or so (and have had that wikipedia page as an open tab for quite a while). So I don't have an opinion, just questions at that this point.

The wikipedia page seems to use the terms power spectrum and power spectral density interchangeably. Quoting from the top of the page: "More commonly used is the power spectral density (or simply power spectrum), which applies ..." So I'm still not sure if the PSD and PS are really two different things. Maybe some technical communities use the terms to mean different things and others do not?

I don't understand how Pbandlmited can be called a power spectrum. As defined, it's an integral and therefore evaluates to a number (assuming the integral exists). I don't see how it's a function of frequency as I would expect for something called a spectrum.

[Jason]

Check this video out: https://www.mathworks.com/videos/understanding-power-spectral-density-and-the-power-spectrum-1707740400025.html

I think it explains some concepts very well!