Difference between theoretical power spectral density and pwelch/periodogram ARMA models
8 次查看(过去 30 天)
显示 更早的评论
Hello,
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1587711/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1587716/image.png)
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);
end
Bn = repmat([1 cell2mat(EstMdl.MA)]', [1 length(theta)]);
Bz= sum(Bn.*epj_b_theta);
else
Bz =1;
end
%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);
end
An = repmat([1 -cell2mat(EstMdl.AR)]', [1 length(theta)]);
Az = sum(An.*epj_a_theta);
else
Az =1;
end
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')
0 个评论
采纳的回答
William Rose
2024-1-10
When estimating the transfer function with pwelch(), remember to compute the denominator.
Estimate the PSD of the input (innovation) signal, e(t).
[x,e,~]=simulate(mdl,1000);
[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=numPSD./denPSD;
H2est should be a good match to the theoretical squared magnitude of the transfer function.
3 个评论
William Rose
2024-1-11
编辑:William Rose
2024-1-11
[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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588691/image.jpeg)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588696/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588701/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588706/image.jpeg)
When the sampling rate is not specified, Matlab's periodogram() and pwelch() assume
radians per unit time. Then
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588711/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588716/image.jpeg)
This is why you observed a mean PSD of approximately 1/3 for an input signal with unit variance.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!