《现代数字信号处理》仿真题(34章).docVIP

3.17 (1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %FFT法与直接法求自相关函数的比较 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc; f1=0.15; f2=0.17; f3=0.26; n=1:32; d=1; %白噪声的方差取1 N=length(n); A1=sqrt((10^3)*2); A2=A1; A3=sqrt((10^2.7)*2); s=A1*cos(2*pi*f1*n+unifrnd(0,2*pi))+A2*cos(2*pi*f2*n+unifrnd(0,2*pi))+A3*cos(2*pi*f3*n+unifrnd(0,2*pi))+normrnd(0,1,1,N); %直接法求自相关函数 r1=zeros(1,N); k=1; for k1=0:N-1 for k2=k1+1:N r1(k)=r1(k)+s(k2)*(s(k2-k1)); end r1(k)=r1(k)/N; k=k+1; end %FFT法求自相关函数 s=[s,zeros(1,N)]; r2=fft(s,2*N); r2=abs(r2).^2/N; r3=ifft(r2); 比较结果：两种方法求出的自相关函数相同。 (2） %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %用BT法和周期图法分别进行功率谱估计 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %用FFT法求出自相关函数 s1=[s,zeros(1,N)]; r2=fft(s1,2*N); r3=abs(r2).^2/N; r=ifft(r3); r4=rot90(rot90(r)); %用BT法进行功率谱估计 M=64; f=-0.5:0.002:0.5; S1=r(1)*ones(1,length(f)); for k=2:M S1=S1+r(k)*exp(-j*2*pi*f*(k-1))+r4(k-1)*exp(j*2*pi*f*(k-1)); end subplot(2,1,1),plot(f,S1,red),grid; title(BT法功率谱估计得到的曲线); %用周期图法进行功率谱估计 U=zeros(1,length(f)); for k=1:N U=U+s(k)*exp(-j*2*pi*f*k); end S2=(abs(U).^2)/N; subplot(2,1,2),plot(f,S2,black),grid; title(周期图法功率谱估计得到的曲线); 仿真结果如下： (3) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %用AR模型法进行功率谱估计 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求自相关函数 r1=zeros(1,N); r2=fft(s,2*N); r3=abs(r2).^2/N; r=ifft(r3); a=zeros(1,p); a1=zeros(1,p); %用Levinson-Drubin迭代算法求解AR模型系数 for m=1:p if m==1 a(1)=-r(2)/r(1); D=r(1)-abs(r(2))^2/r(1); else k=0; for i=1:m-1 k=k+a(i)*r(m+1-i); end k=-(k+r(m+1))/D; a1=a; a(m)=k; for i=m-1:-1:1 a(i)=a1(i)+k*a1(m-i); end D=D*(1-abs(k)^2); end end f=-0.5:0.002:0.5; S=zeros(1,length(f)); S1=zeros(1,length(f)); for i=1:p S1=S1+a(i)*exp(-j*2*pi*f*i); end S1=abs(S1+1).^2; S=D./S1; plot(f,S),grid; title(16阶AR模型法功率谱估计曲线);