随机信号实验.docxVIP

求x(n)=sin(10πn)+sin(24πn)的功率谱密度：一：直接做傅里叶变换，再去模平方x=sin(10*pi*t)+sin(24*pi*t);Fs=1024; N=1024; Cx=fft(x,N);Y=abs(Cx)Pxx=Y.*conj(Y)/N; %power spectrumf=Fs*(0:(N/2))/N;plot(f,Pyy(1:(N/2+1)))title(功率谱)xlabel(频率/(Hz))二：求自相关函数，在做傅里叶变换N=1024;idx=0:round(N/2-1);k=idx/N;n=1:Nxn=sin(10*pi*n)+sin(24*pi*n);yn=xcorr(xn,unbiased);subplot(211),plot(n,yn),title(自相关函数);Yk=fft(yn,N);Px=abs(Yk);subplot(212),plot(k,10*log10(Px(idx+1)))title(功率谱密度);三．Burg法（1）算法主要思想Burg法不是直接估计AR模型的参数，而是先估计反射系数。使用线性预测的方法来计算不同阶数下的反射系数，其同时使用前向和后向线性预测，使前向和后向预测误差的平均功率相对各阶反射系数最小，将反射系数代入Levinson-Durbin公式即可求解。（2）运算简要框图用L-D公式求解估计反射系数X(n) 输出N=512;n=0:N-1;xn=2*sin(10*pi*n)+ cos(24*pi*n);Rx=zeros(1,N+1);%从课本上的公式来看，Rx（m）中的m属于（0，m），即共有m+1个，故在这里设Rx是一个一行，N+1列的向量figure(1)plot(n,xn);title(（a）两正弦信号加白噪声波形)%下面用书中所讲自相关函数估计中的渐进无偏估计来估计自相关函数for m=1:N+1;%由于在matlab中，下角标不能是0，m属于（0，m），在此只能从1到N+1 sum=0; for n=1:(N+1-m);%同样道理，把书中公式里m换成m-1,N换成N+1，求和下限变为1 sum=sum+xn(n).*xn(n+m-1); end Rx(m)= sum/N;%切记，这里的Rx(1)才是自相关函数在0点的取值。Rx（m）只是一个存储数据的代号，为了跟书中公式一致，才叫RxendP=50;a=zeros(P,P);%a中有两个变量m,i，所以设a是P行P列的向量p=zeros(1,P+1);%由于matlab中没有ρ，故用p来代替表示，ρ的范围是（0，P）共有P+1个%下面计算AR模型参数a(1,1)=-Rx(2)/Rx(1);ef=zeros(P,N);eb=zeros(P,N);ef(1,:)=xn;eb(1,:)=xn;for m=2:P+1; km1=0;km2=0; for n=m:N km1=km1+ef(m-1,n).*eb(m-1,n-1); km2=km2+(ef(m-1,n)).^2+(eb(m-1,n-1)).^2; end a(m,m)=(-2)*km1./km2; for n=m:N ef(m,n)=ef(m-1,n)+a(m,m).*eb(m-1,n-1); eb(m,n)=eb(m-1,n-1)+a(m,m).*ef(m-1,n);end endp(1)=Rx(1);p(2)=Rx(1).*(1-abs(a(1,1).^2));a=a(2:P+1,2:P+1);for m=2:P %由于m=1时的各个值在上面已经给出，故从m=2开始求 for i=1:m-1 a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i); end p(m)=p(m-1)*(1-abs(a(m,m)).^2);endz=[1,a(P,:)];G=sqrt(p(P));[H, w]=freqz(G,z,512);%调用计算数字滤波器频响的函数figure(2)plot(w/(2*pi),10*log10(abs(H).^2));ylabel(10log(PSD))title((b)Burg法估计功率谱密度 )比较：第一种方法得到的功率谱分辨率较高，但方差性较差，功率谱起伏剧烈，容易出现虚假谱峰。第二种方法得到的功率谱方差性较好，功率谱比第一种方法估计的要平滑，分辨率比第一种低。第三种方法得到的功率谱比前两种都平滑，而且分辨率最高。但当阶数越高时，得到的谱的分辨率也越高，但方差也越大，将会产生更多的虚假谱峰。总结：当采样的点数减少的时候，古典法的估计就十分不准确了，而参数法功率谱估计却还可以让人清晰地分辨两个挨得较近的谱峰。优点十分突出。