0
$\begingroup$

I'm driving a white noise signal $e(n)$ through an AR random process:

$$ x(n)=e(n)+\sum\limits_{k=1}^N \alpha_k x(n-k) $$

I got five values of $\alpha_k$ and plotted its power spectrum in Matlab for $\omega$ from $0$ to $\pi$.

Then I generated the realizations of $x(n)$ with 10000 samples

Now I'm trying to get the estimate of the spectrum from the expectation data using the periodogram formula

$$ S_{per}(\omega)= \frac{1}{K} \left|\sum\limits_{k=0}^{K-1}x[k]e^{-j\omega k} \right|^2 $$

with $K=256$ the number of first values of $x(n)$. When I run my code what I get is a vector with the same value in every position, so I'm plotting a straight line. I think that just applying the formula I should get the estimate, so what am I doing wrong? :-S

This is the section of the code I'm using to find the estimates

var_e=0.01;                                  % noise variance
b_1=1.4261;                                  % coefficients
b_2=-.7634;
b_3=-.9002;
b_4=1.2548;
b_5=-.5707;    
K=256;
wn=rand(1,10000)-0.5;                        % white noise
samp=1:1:1e4;
T=10000;
x=zeros(1,length(samp));
% here I generate the expectations
for n=6:T+5
    x(n)=b_1*x(n-1)+b_2*x(n-2)+b_3*x(n-3)+b_4*x(n-4)+b_5*x(n-5)+wn(n-5);
end
x=x(6:end);
figure
plot(samp,x)
    for i=1:length(w)                         % w is the frequency, has 1000 
        for k=1:K
        per2(i)=sum(x(k)*exp(-j*k*i));        % sum of exponentials
        pe2(i)=per2(i)*conj(per2(i));         % squared
        end
    end

Thanks in advance for your help

$\endgroup$
0
$\begingroup$

I think that is because your exponentials are wrong.
There is no variable named w in your code, although you refer to it, so I believe it should be something like:

per2(i)=sum(x(k)*exp(-2j*pi*k*w(i))));

But consider using fft instead:

pe2 = abs(fft(x)).^2;

or

pe2 = abs(fft(x(1:K))).^2;

if you insist on using only $K$ samples for some reason (why generate 10000 sample then??).

Perhaps you should consider the Welch method, with length 256

pe2 = pwelch(x,256,128,256);

BTW - you can generate AR process using filter:

x = filter(1,a,wn);

where a is the vector of AR coefficients.

$\endgroup$
0
$\begingroup$

I think the mistake is in the for loop
instead put it as pe2=zeros(K,1); for i=1:K for l=1:K per2(k)=(x(k)*exp(-1i*double(k*l*2)*pi/K)); end pe2(i)=(abs(sum(per2))^2; end plot(pe2); this gives the PSD plot for x(1:256)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.