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

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:

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

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

Answer 1

[William Rose]

@Simona Turco,

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.

[William Rose]

@Simona Turco,

Here is a simple script that estimates the frequency response of a simple AR model. The model is

whose z-transform is

therefore

therefore

The script below simulates the model, with a random Gaussian input signal. Then it uses the simulated input and output to estimate , and it computes the theoretical squared transfer function magnitude.

N=1000;                     % signal duration (points)
ar=0.5;                     % autoregressive coefficient
x=randn(1,N);               % Gaussian input signal
y=zeros(1,N);               % allocate array for output
for i=2:N
    y(i)=x(i)+ar*y(i-1);    % compute output
end
% Estimate transfer function
[denw,w]=pwelch(x);                 % denominator PSD (radians)
[numw,w]=pwelch(y);                 % numerator PSD (radians)
psd=numw./denw;                     % PSD
% Compute theoretical transfer function
Hw2=(abs(1./(1-ar*exp(-1i*w)))).^2; % theo |H(w)|^2
% plot results
figure; plot(w,psd,'r.',w,Hw2,'-b');
grid on; xlabel('\omega (radians/sample)'); ylabel('|H(\omega)|^2')
legend('PSD','H^2_t_h_e_o')

The estimate matches the theoretical value.

[Simona Turco]

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]

@Simona Turco,

[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.