dsp2.m

和先前传的程序是配套的

%  姓名:齐家红 学号：200520203020
clc
clear
%可以在此改变参数p的值
%p=4;
%p=10;
p=20;
snr=10;w0=1.2*pi;p2=1;
c=p2/(10^(snr/10));
v=randn(1000,1)*c^0.5;
%2.1
for n=1:1:1000
x(n)=exp(j*pi*(n-1))+exp(j*w0-j*0.7*pi)+v(n);
end
x;
%2.2
r1=0;
r=zeros(501,1);
for l=2:1:501
    for k=(l+1):1:1000
        r1=x(k)*conj(x(k-l))+r1;
    end
    r(l)=r1/1000;
    r1=0;
end
r12=0;
for k1=1:1:1000
    r12=r12+x(k1)*conj(x(k1));
end
r(1)=r12/1000;
r;
%2.3
w=0:pi/1000:2*pi;
pt1=0.0000000;
for k1=2:1:501
  pt1=pt1+real(r(k1)*exp(-i*w*k1));
end
pt=r(1)+2*pt1;
subplot(3,1,1);
plot(w,pt);
title('DSP作业第二次图  姓名:齐家红 学号：200520203020');
%2.4
p1=4;
aspa=zeros(p,p);
for m=1:1:p1
    cl(m)=r(m);
end
R=toeplitz(cl,cl);
B=eye(p1);
aspa=R-c*B;
[vi,ti]=eig(aspa);
syms w
for m1=1:1:p1
    a(m1)=exp(-w*(m1-1)*i);
end
v1=vi(:,p1-1:p1)
d1=(abs(a*vi(:,p1-1)))^2;
d2=(abs(a*vi(:,p1)))^2;
Pmusic=1/(d1+d2);
Pmusic;
subplot(3,1,2);
ezplot(Pmusic,[0,2*pi]);
%2.5;
r11=zeros(1,p);
for j=1:1:p
   r11(j)=r(j+1);    
end
r2=[r(1),r11];
R1=[r(1),r11'.'];
R2=[r2];
R5=toeplitz(R1,R2);
syms a1;
for l=1:p+1
    a1(l)=exp(-w*(l-1)*i);
end
a1=a1.';
Pmvdr=1/(a1'*inv(R5)*a1);
subplot(3,1,3)
ezplot(Pmvdr,[0,2*pi]);