现代数字信号处理理论计算法仿真作业.docVIP

现代数字信号处理仿真作业 3.17 首先定义函数 u_length(N)，函数的作用是产生长度为N的观测样本u(n),程序如下： function u = length(N); A1 = sqrt(10^(30/10)); A2 = sqrt(10^(30/10)); A3 = sqrt(10^(27/10)); n = 1:N; s_signal = A1*exp(i*2*pi*0.15*n)+A2*exp(i*2*pi*0.17*n)+A2*exp(i*2*pi*0.26*n); v_noise = sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N);; u = s_signal+v_noise; (1)N=32，分别用给予FFT的自相关函数和式3.1.2估计自相关函数r(m),程序如下： function r11 N=32; u = u_length(N); %长度为32的随机信号 for m = 1:N r11(m) = u(m:N)*u(1:N-m+1)/N; end r12 = conj(r11(N:-1:2)); %自相关函数的对称性 r1 = [r12,r11]; %r1是按式3.1.2估计自相关函数r(m) u2N = [u,zeros(1,N)]; U_w = fft(u2N,2*N); S_BT = (abs(U_w).^2)/N; r20=ifft(S_BT); r2=[r20(N+2:2*N) r20(1:N)]; %r2是FFT的自相关函数快速算法得出的自相关函数 hold on k = -N+1:N-1; plot(k,abs(r1),r); plot(k,abs(r2),g*); grid on 仿真图如下： （2）N=256，分别用BT法和周期图发估计u(n)的功率谱 %BT法估计功率谱密度，自相关采用M=64点 N = 256; M=64; u = u_length(N); u_xcorr = xcorr(u,M,biased); m = -M:M; k = -100:100; h = u_xcorr*exp(-j*pi/100).^(m*k); Sw_BT = 10*log10(abs(h)./max(abs(h))); plot(k/200,Sw_BT); title(BT法估计功率谱密度)； grid on BT法估计功率谱密度仿真结果如下图： %周期图法估计功率谱密度 Uw_N1=abs(fft(u,N)); Uw_N=fftshift(Uw_N1); w=1:N; plot(w/N-0.5,10*log10(abs(Uw_N)./max(abs(Uw_N)))); title(周期图法估计结果) 周期图法估计功率谱密度仿真结果如下图： （3）L—D迭代算法求解AR模型参数及功率谱估计，程序如下： N=256; p=16; u = u_length(N) r=xcorr(u,p,biased); r=r(p+1:2*p+1); a=-r(2)/r(1); c=r(1)-r(2)*r(2)/r(1); for m=2:p l=conj(r(m:-1:2)); k=-(r(m+1)+a*l)/c; for n=1:m-1 b(n)=a(n)+k*conj(a(m-n)); end a=[b k]; c=c*(1-k*k); end h=freqz(sqrt(c),[1 a],201,whole); h=fftshift(h); b=abs(h).^2; n=-100:100; plot(n/200,10*log10(b/max(b))); grid title(阶数为16的AR模型功率谱密度估计曲线) 仿真结果如下图： 3.20 clear; N=1000; M=8;K=2; v=sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N); n=1:N; fi=2*pi*rand(1,2); u(n)=exp(j*0.5*pi*n+j*fi(1))+exp(j*1.2*pi*n+j*fi(2))+v(n); r=xcorr(u,M,biased); for n=1:M for m=1:M R(m,n)=r(M+1+n-m); end end [G,D]=eigs(R,M-K,sm); for n=-100:100; for m=1:M a(m)=exp(j*(m-1)*pi*n/100); end p(n+101)=1/(a*(G*G)*a); end p=ab