Clear Filters
Clear Filters

Difference between theoretical power spectral density and pwelch/periodogram ARMA models

14 views (last 30 days)
I am calculating the power spectral density (PSD) of ARMA models in two different ways:
1) By simulating the ARMA model and using Welch's periodogram
2) By using the theoretical PSD form of an ARMA model:
Strangely, enough the two functions seems to be off always by a factor 3 in amplitude, no matter the length of the simulated signal input to the periodogram, nor the order of the ARMA model.
Where is this factor 3 coming from?
Here is a plot of the resulting PSD
Below the code I am using for testing:
%simulate ARMA signal
N = 1000;
mdl = arima('Constant',0,'Variance',1, 'MA',{0.8}, 'AR',{-0.75, - 0.25});
x = simulate(mdl,1000);
EstMdl = mdl;
%calculate periodogram
[psd0, w] = pwelch(x);
figure, plot(w,psd0)
theta = linspace(0,pi,1000);
% Calculate PSD by ARMA theoretical PSD
%form numerator
if ~isempty(EstMdl.MA)
bn = 1:size(EstMdl.MA,2);
epj_b_theta = ones(size(EstMdl.MA,2)+1,length(theta));
for n=1:size(EstMdl.MA,2)
epj_b_theta(n+1,:) = exp(-1j*n*theta);
Bn = repmat([1 cell2mat(EstMdl.MA)]', [1 length(theta)]);
Bz= sum(Bn.*epj_b_theta);
Bz =1;
%form denominator
if ~isempty(EstMdl.AR)
an = 1:size(EstMdl.AR,2);
epj_a_theta = ones(size(EstMdl.AR,2)+1,length(theta));
for n=1:size(EstMdl.AR,2)
epj_a_theta(n+1,:) = exp(-1j*n*theta);
An = repmat([1 -cell2mat(EstMdl.AR)]', [1 length(theta)]);
Az = sum(An.*epj_a_theta);
Az =1;
Hz = Bz./Az;
psd2 = EstMdl.Variance.*abs(Hz).^2;
% Plot
hold on, plot(theta, psd2)
hold on, plot(theta, psd2/3)
xlabel('\theta [rad]')
legend('pwelch', 'PSD_A_R_M_A', 'PSD_A_R_M_A/3')

Accepted Answer

William Rose
William Rose on 10 Jan 2024
When estimating the transfer function with pwelch(), remember to compute the denominator.
Estimate the PSD of the input (innovation) signal, e(t).
[numPSD, w] = pwelch(x);
[denPSD, w] = pwelch(e);
Compute the squared transfer funciton estimate as the ratio of the output to the input PSD estimates:
H2est should be a good match to the theoretical squared magnitude of the transfer function.
Simona Turco
Simona Turco on 11 Jan 2024
Edited: Simona Turco on 11 Jan 2024
Thank you for your answer. As you explained, dividing for the spectrum of the input signal does the trick, but I am trying to wrap my head around why that is the case. Theoretically, the input is white noise, so its spectrum is "flat" at the value of the variance. So this factor should be captured by the variance multiplying the squared magnitude of the transfer function in the theoretical formula.
When I plot only the denominator PSD (the PSD of the input), I would expect it to be fluctuating at around 1, as the input variance in the ARMA model is set to 1, but instead is fluctuating at around 1/3. Why is that?
William Rose
William Rose on 11 Jan 2024
Edited: William Rose on 11 Jan 2024
[edit: Fix typo: ponts -> points]
When the input sequence is real, Matlab's periodogram() and pwelch() return a one-sided power spectral density (PSD). If a rectangular window is used, the 1-sided PSD is
where fs =sampling rate, N=sequence length in points, and X=fft(x). For a white signal, the expected value of X(k) is a constant: . From Parseval's relation, the expected value is . Therefore the expected value of the one-sided PSD of a random sequence is
When the sampling rate is not specified, Matlab's periodogram() and pwelch() assume radians per unit time. Then
This is why you observed a mean PSD of approximately 1/3 for an input signal with unit variance.

Sign in to comment.

More Answers (0)




Community Treasure Hunt

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

Start Hunting!