Getting wrong periodogram

Question

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

Answer 1

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.

Answer 2

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)