现代数字信号处理仿真作业
现代数字信号matlab处理仿真题
3.17（1）相关函数仿真代码：A1=getAk(SNR1);A2=getAk(SNR2);A3=getAk(SNR3); %求得信号的幅度；noise=randn(1,N) + j*randn(1,N); % 构建高斯白噪声；s1=getSk(A1,f1,N);s2=getSk(A2,f2,N);s3=getSk(A3,f3,N);; %产生3个复正弦信号vn=s1+s2+s3+noise;vk=fft(vn,2*N); %对v(n)补N个零，然后做2N点FFTswk=((abs(vk)).^2)/N; %计算功率谱估计S(ω)r0=ifft(swk); %对S(k)做ifft得到r=[r0(N+2 : 2*N) , r0(1 : N)]; %根据教程3.1.8式可得r1=xcorr(vn, N-1,'biased'); %直接计算自相关函数%%%%%%%%%%%%%%%%%%%%%%%%%取序列实部，虚部%%%%%% real_r=real(r);imag_r=imag(r);real_r1=real(r1);imag_r1 = imag(r1);subplot(2,2,1);stem(real_r);xlabel('基于FFT的自相关函数的实部');ylabel('实部');subplot(2,2,2);stem(imag_r);xlabel('基于FFT的自相关函数的虚部');ylabel('虚部');subplot(2,2,3);stem(real_r1);ylabel('实部');xlabel('估计的自相关函数的实部');subplot(2,2,4);stem(imag_r1);xlabel('估计的自相关函数的虚实部');ylabel('虚部');function AK=getAk(SNR) %求得幅度%%%%%%%%%%%%%%%%%%%由SNR=10log(A^2/2*σ^2) %%%%%%% AK=((10^(SNR/10))*2)^0.5;function Sk=getSk(Ak,fk,N)Sk=Ak * exp(j * 2 * pi * fk *(0:N-1));仿真波形：（2）BT 法和周期法估计 仿真程序： clear all; clc;%设定N 值可以改变抽样信号的点数，设定M 值可以设定加窗的大小，设定N3可以补零，确定实际求fft 的点数。
现代数字信号处理及应用仿真题答案
仿真作业姓名：***学号：S*********4.17程序clc;clear;for i=1:500sigma_v1=0.27;b(1)=-0.8458;b(2)=0.9458;a(1)=-(b(1)+b(2));a(2)=b(1)*b(2);datlen=500;rand('state',sum(100*clock));s=sqrt(sigma_v1)*randn(datlen,1);x=filter(1,[1,a],s);%%sigma_v2=0.1;u=x+sqrt(sigma_v2)*randn(datlen,1);d=filter(1,[1,-b(1)],s);%%w0=[1;0];w=w0;M=length(w0);N=length(u);mu=0.005;for n=M:Nui=u(n:-1:n-M+1);y(n)=w'*ui;e(n)=d(n)-y(n);w=w+mu.*conj(e(n)).*ui;w1(n)=w(1);w2(n)=w(2);ee(:,i)=mean(e.^2,2);endendep=mean(ee');plot（ep）;xlabel('迭代次数');ylabel('MSE');title('学习曲线'); plot(w1);hold;plot(w2);仿真结果：步长0.015仿真结果00.10.20.30.40.50.60.7迭代次数M S E 学习曲线步长0.025仿真结果步长0.005仿真结果4.18 程序data_len = 512; %样本序列的长度trials = 100; %随机试验的次数A=zeros(data_len,2);EA=zeros(data_len,1);B=zeros(data_len,2);EB=zeros(data_len,1);for m = 1: trialsa1 = -0.975;a2 = 0.95;sigma_v_2 =0.0731;v = sqrt(sigma_v_2) * randn(data_len, 1, trials);%产生v(n)u0 = [0 0];num = 1;den = [1 a1 a2];Zi = filtic(num, den, u0); %滤波器的初始条件u = filter(num, den, v, Zi); %产生样本序列u(n)%（2）用LMS滤波器来估计w1和w2mu1 = 0.05;mu2 = 0.005;w1 = zeros(2, data_len);w2 = zeros(2, data_len);e1 = zeros(data_len, 1);e2 = zeros(data_len, 1);d1 = zeros(data_len, 1);d2 = zeros(data_len, 1);%LMS迭代过程for n =3 :data_len - 1w1( :, n+1) = w1( :, n) + mu1 * u(n-1 : -1: n-2, : , m) * conj(e1(n));w2( :, n+1) = w2( :, n) + mu2 * u(n-1 : -1: n-2, : , m) * conj(e2(n));d1(n+1) = w1( : , n+1)' * u(n: -1: n-1, :, m);d2(n+1) = w2( : , n+1)' * u(n: -1: n-1, :, m);e1(n+1) = u(n+1, : ,m) - d1(n+1);e2(n+1) = u(n+1, : ,m) - d2(n+1);endA = A + conj(w1)';EA = EA +e1.^2;B = B + conj(w2)';EB = EB + e2.^2;end%剩余均方误差和失调参数wopt=zeros(2,trials);Jmin=zeros(1,trials);sum_eig=zeros(trials,1);for m=1:trials;rm=xcorr(u(:,:,m),'biased');R=[rm(512),rm(513);rm(511),rm(512)];p=[rm(511);rm(510)];wopt(:,m)=R\p;[v,d]=eig(R);Jmin(m)=rm(512)-p'*wopt(:,m);sum_eig(m)=d(1,1)+d(2,2);endsJmin=sum(Jmin)/trials;e1_100trials_ave=sum(e1)/trials;e2_100trials_ave=sum(e2)/trials;Jex1=e1_100trials_ave-sJmin;Jex2=e2_100trials_ave-sJmin;sum_eig_100trials=sum(sum_eig)/100;Jexfin=mu1*sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials)); Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials)); M1=Jexfin/sJminM2=Jexfin2/sJminfigure(1);plot(A/trials);hold on;plot(conj(w1)');xlabel('迭代次数');ylabel('权向量');title('步长为0.05权向量收敛曲线');figure(2);plot(B/trials);hold on;plot(conj(w2)');xlabel('迭代次数');ylabel('权向量');title('步长为0.005权向量收敛曲线');figure(3);plot(EA/trials,'*');hold on;plot(EB/trials,'-');xlabel('迭代次数');ylabel('均方误差');title('步长分别为0.05和0.005学习曲线'); 仿真结果失调参数M1= 0.0545 M2= 0.00524.19程序clear all%产生观测信号和期望信号trials = 100; %随机试验的次数data_len = 1000; %样本数目n =1 : data_len;A1 = zeros(data_len, 2);EA1 = zeros(data_len, 1);for i = 1: trialssigma_v_2 = 0.5;phi = 2 * pi * rand(1, 1); %随机相位signal = sin(pi/2 * n' +phi); %信号s(n)u = signal + sqrt (sigma_v_2) * randn(data_len, 1); %观测信号u(n)d = 2 * cos(pi/2 * n' +phi); %期望响应信号d(n)%LMS迭代算法mu = 0.015;M = 2;w = zeros(M,data_len);e = zeros(data_len,1);y = zeros(data_len,1);for m = 2: data_len-1w(:, m + 1) = w(: , m) + mu * u(m: -1: m - 1) * conj(e(m));y(m + 1) = w(: , m + 1)' * u(m + 1:-1: m);e(m + 1) = d(m + 1) - y(m + 1);endA1 = A1 + conj(w)';EA1 = EA1 +e.^2;endfigure(1);plot(e);xlabel('迭代次数');ylabel('均方误差');title('单次实验学习曲线');figure(2);plot(EA1/trials);xlabel('迭代次数');ylabel('均方误差');title('100次独立试验学习曲线');figure(3);plot(A1/trials);hold on;plot(conj(w)');xlabel('迭代次数');ylabel('权向量');title('权向量收敛曲线'); 仿真结果：5.10(1)247.04846.5783 46.578347.0487R⎡⎤=⎢⎥⎣⎦(2)347.048746.578346.1125 46.578347.048746.5783 46.112546.578647.0487R⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦(3) 特征值分解eig(R2)=diag{0.4704,93.6270}Eig(R3)=diag{0.3148,0.9362,139.8951}特征值扩展：X(R2)=199.0370X(R3)=444.4107（4）程序clear allclc;L=10000;sigma_v1=0.93627;A1 = zeros(L, 2);EA1 = zeros(L, 1);for i=1:100v=sqrt(sigma_v1)*randn(L,1);a1=-0.99;u(1)=v(1);for k=2:Lu(k)=-a1*u(k-1)+v(k);end% u=u(500:end);M=2;w(1,:)=zeros(1,M);e(1)=u(1);mu=0.001;uu=zeros(1,M);w(2,:)=w(1,:)+mu*e(1)*uu;uu=[u(1) uu(1:M-1)];dd=(w(2,:)*uu')';e(2)=u(2)-dd;for k=3:Lw(k,:)=w(k-1,:)+mu*e(k-1)*uu;uu=[u(k-1) uu(1:M-1)];dd=(w(k,:)*uu')';e(k)=u(k)-dd;endA1 = A1 + conj(w);EA1 = EA1 +(e.^2)';endfigure(1);plot(EA1/100);xlabel('迭代次数');ylabel('均方误差');title('迭代500次，步长0.001');figure(2);plot(A1/100);hold on;plot(conj(w));xlabel('迭代次数');ylabel('权向量');title('权向量收敛曲线');5.11clear allclear;clc;for i=1:1500N=1000;M=5;L=2;h=[0.389 1 0.389];sigma=1e-3;vn=sqrt(sigma)*randn(2*M+N,1); H=zeros(2*M+1,2*M+L+1);for k=1:2*M+1H(k,k:1:k+L)=h;ends=randsrc(2*M+L+N,1);S=zeros(2*M+L+1,N);V=zeros(2*M+1,N);for k=1:NS(:,k)=s(2*M+L+k:-1:k);V(:,k)=vn(2*M+k:-1:k);endU=H*S+V;dn=S(M+L+1,:);if (i<=500)mu=0.01;elseif (i>500&&i<=1000)mu=0.025;elsemu=0.05;enda=size(U);M=a(1);N=a(2);err=zeros(N,1);w=zeros(M,N);w((M-1)/2+1,1)=1;err(1)=dn(1)-w(:,1)'*U(:,1);for k=1:N-1w(:,k+1)=w(:,k)+mu*U(:,k)*conj(err(k)); err(k+1)=dn(k+1)-w(:,k+1)'*U(:,k+1);endif (i<=500)ee1(:,i)=mean(abs(err).^2,2);elseif (i>500&&i<=1000)ee2(:,i)=mean(abs(err).^2,2);elseee3(:,i)=mean(abs(err).^2,2);endendep1=mean(ee1');ep2=mean(ee2');ep3=mean(ee3');figure(1);plot(ep1);hold on;plot(ep2);hold on;plot(ep3)xlabel('µü´ú´ÎÊý');ylabel('¾ù·½Îó²î');。
现代数字信号处理及其应用 第六章仿真题
6.12clc;clear;M=15;Lb=10;% hb=[0.407 0.815 0.407];hb=[0.04 -0.05 0.07 -0.21 -0.5 0.72 0.36 0 0.21 0.03 0.07]; Hb=zeros(2*M+1,2*M+Lb+1);for k=1:2*M+1Hb(k,k:1:k+Lb)=hb;endEA1 = zeros(2000, 1);EA2 = zeros(2000, 1);for k=1:100sigma=1e-3;N=2000;s=randsrc(2*M+Lb+N,1);vn=sqrt(sigma)*randn(2*M+N,1);S=zeros(2*M+Lb+1,N);V=zeros(2*M+1,N);for k=1:NS(:,k)=s(2*M+Lb+k:-1:k);V(:,k)=vn(2*M+k:-1:k);endUb=Hb*S+V;errb_LMS=zeros(N,1);wb_LMS=zeros(2*M+1,N);wb_LMS(M+1,1)=1;dn=S(M+Lb+1,:);errb_LMS(1)=dn(1)-wb_LMS(:,1)'*Ub(:,1);mu=0.025;for k=1:N-1wb_LMS(:,k+1)=wb_LMS(:,k)+mu*Ub(:,k)*conj(errb_LMS(k)); errb_LMS(k+1)=dn(k+1)-wb_LMS(:,k+1)'*Ub(:,k+1);endMSEb_LMS=abs(errb_LMS).^2;EA1=EA1+MSEb_LMS;lambda=0.99;delta=0.004;wb_RLS=zeros(2*M+1,N+1);wb_RLS(M+1,1)=1;epsilon=zeros(N,1);P1=eye(2*M+1)/delta;for k=1:NPIn=P1*Ub(:,k);deno=lambda+Ub(:,k)'*PIn;kn=PIn/deno;epsilon(k)=dn(k)-wb_RLS(:,k)'*Ub(:,k);wb_RLS(:,k+1)=wb_RLS(:,k)+kn*conj(epsilon(k)); P1=P1/lambda-kn*Ub(:,k)'*P1/lambda;endMSEb_RLS=abs(epsilon).^2;EA2=EA2+MSEb_RLS;endM=15;Lb=2;hb=[0.407 0.815 0.407];Hb=zeros(2*M+1,2*M+Lb+1);for k=1:2*M+1Hb(k,k:1:k+Lb)=hb;endEA3 = zeros(2000, 1);EA4 = zeros(2000, 1);for k=1:100sigma=1e-3;N=2000;s=randsrc(2*M+Lb+N,1);vn=sqrt(sigma)*randn(2*M+N,1);S=zeros(2*M+Lb+1,N);V=zeros(2*M+1,N);for k=1:NS(:,k)=s(2*M+Lb+k:-1:k);V(:,k)=vn(2*M+k:-1:k);endUb=Hb*S+V;errb_LMS=zeros(N,1);wb_LMS=zeros(2*M+1,N);wb_LMS(M+1,1)=1;dn=S(M+Lb+1,:);errb_LMS(1)=dn(1)-wb_LMS(:,1)'*Ub(:,1);mu=0.025;for k=1:N-1wb_LMS(:,k+1)=wb_LMS(:,k)+mu*Ub(:,k)*conj(errb_LMS(k)); errb_LMS(k+1)=dn(k+1)-wb_LMS(:,k+1)'*Ub(:,k+1);endMSEb_LMS=abs(errb_LMS).^2;EA3=EA3+MSEb_LMS;lambda=0.99;delta=0.004;wb_RLS=zeros(2*M+1,N+1);wb_RLS(M+1,1)=1;epsilon=zeros(N,1);P1=eye(2*M+1)/delta;for k=1:NPIn=P1*Ub(:,k);deno=lambda+Ub(:,k)'*PIn;kn=PIn/deno;epsilon(k)=dn(k)-wb_RLS(:,k)'*Ub(:,k);wb_RLS(:,k+1)=wb_RLS(:,k)+kn*conj(epsilon(k));P1=P1/lambda-kn*Ub(:,k)'*P1/lambda;endMSEb_RLS=abs(epsilon).^2;EA4=EA4+MSEb_RLS;end% figureplot(EA1/100);hold onplot(EA2/100);hold onplot(EA3/100);hold onplot(EA4/100);6.15clc;clear;EA1 = zeros(999, 1);A1 = zeros(2, 1000);for i=1:100a1=0.99;sigma=0.995;N=1000;vn=sqrt(sigma)*randn(N,1);nume=1;deno=[1 a1];u0=zeros(length(deno)-1,1);xic=filtic(nume,deno,u0);un=filter(nume,deno,vn,xic);n0=1;M=2;b=un(n0+1:N);L=length(b);un1=[zeros(M-1,1).',un.'];A=zeros(M,L);for k=1:LA(:,k)=un1(M-1+k:-1:k);enddelta=0.004;lambda=0.98;w=zeros(M,L+1);epsilon=zeros(L,1);P1=eye(M)/delta;for k=1:LPIn=P1*A(:,k);denok=lambda+A(:,k)'*PIn;kn=PIn/denok;epsilon(k)=b(k)-w(:,k)'*A(:,k);w(:,k+1)=w(:,k)+kn*conj(epsilon(k)); P1=P1/lambda-kn*A(:,k)'*P1/lambda; endMSE=abs(epsilon).^2;EA1=EA1+MSE;A1=A1+w;endplot(EA1/100);A2=A1/500;plot(A2(1,:));。
ADSP现代数字信号处理仿真实验报告
目录仿真一：LMS算法和RLS算法 (1)1 自适应滤波的基本原理 (1)1.1 自适应最小均方（LMS）算法 (1)1.2 递归最小二乘方（RLS）算法 (2)2 仿真实验 (4)3 结果分析 (6)仿真二：P阶Levinson-Durbin算法 (8)1 要求： (8)2 算法描述 (8)2.1 产生信号 (8)2.2 L-D算法 (9)2.3 对比信号谱功率和LD算法谱估计 (10)3 结果分析 (11)3.1 AR模型 (11)3.2 MA模型 (12)3.3 总结 (13)仿真一：LMS 算法和RLS 算法1 自适应滤波的基本原理自适应滤波器由参数可调的数字滤波器/自适应处理器和自适应算法两部分组成，如图1所示。
输入信号x(n)通过参数可调数字滤波器后产生的输出信号为y(n)，将其与参考信号d(n)进行比较，得到误差信号e(n)。
误差信号e(n)经过一定的自适应算法后反馈到参数可调数字滤波器，对滤波器进行参数调整（有时还需要利用x(n)），以使得e(n)最终的均方值最小。
这是一种自动控制理论，因此，滤波器在设计时不需要事先知道输入信号和噪声的统计特性，而能够根据输入信号的统计特性变化自动跟踪这种变化，自动调整参数，使滤波器性能达到最佳。
图 1 自适应滤波器框图图1所示自适应滤波器，输入信号为：x(n)和d(n)，两个输出为：y(n)和e(n)。
当误差信号e(n)的均方误差达到最小的时候，可以证明信号y(n)是信号d(n)的最佳估计。
1.1 自适应最小均方（LMS ）算法最陡下降法每次迭代都需要知道性能曲面上某点的梯度值，而梯度值只能根据观测数据进行估计。
LMS 算法是一种有用简单的估计梯度的方法，其最核心的思想是采用平方误差最小代替均方误差最小准则。
信号基本关系：()()()()()()(1)()2()()T y n W n X n e n d n y n W n W n e n X n μ==-+=+式中，W(n) 为 n 时刻自适应滤波器的权矢量，011()[(),(),....()]TN W n w n w n w n -=，下一时刻权矢量 W(n +1) 等于当前权矢量 W (n ) 加上一个修正量，该修正量是误差信号e (n )的加权值，加权系数 2μx(n) 正比于当前的输入信号 x(n)。
现代信号处理大作业
姓名：潘晓丹学号：班级：A1403492作业1LD 算法实现AR 过程估计1.1 AR 模型p 阶AR 模型的差分方程为：)()()(1n w i nx a n x pi i ，其中)(n w 是均值为0的白噪声。
AR 过程的线性预测方法为：先求得观测数据的自相关函数，然后利用Yule-Walker 方程递推求得模型参数，再根据公式求得功率谱的估计。
Yule-Walker 方程可写成矩阵形式：1.2 LD 算法介绍Levinson-Durbin算法可求解上述问题，其一般步骤为：1) 计算观测值各自相关系数pjjrxx,,1,0),(；)0(0xxr；i=1；2) 利用以下递推公式运算：3) i=i+1，若i>p，则算法结束；否则，返回(2)。
1.3 matlab编程实现以AR模型：为例，Matlab 程序代码如下：clear; clc;var = 1;noise = var*randn(1,10000);p = 2;coefficient = [1 -0.5 0.5];x = filter(1,coefficient,noise);divide = linspace(-pi,pi,200);for ii = 1:200w = divide(ii);S1(ii) = var/(abs(1+coefficient(2:3)*exp(-j*w*(1:2))'))^2;end[a_p var_p]=Levinson_Durbin(x,p);for ii = 1:200w = divide(ii);Sxx(ii) = var_p/(abs(1+a_p(2:p+1)*exp(-j*w*(1:p))'))^2;endfigure;subplot(2,2,1);plot(divide,S1,'b');grid onxlabel('w');ylabel('功率');title('AR 功率谱');subplot(2,2,2);plot(divide,Sxx,'r-');grid onxlabel('w');ylabel('功率');title('L-D算法估计');subplot(2,2,3);plot(divide,S1,'b');hold onplot(divide,Sxx,'r--');hold offgrid onxlabel('w');ylabel('功率');title('AR功率谱和算法比较');子函数：Levinson_Durbin.mfunction [a_p var_p] = Levinson_Durbin(x,p)N = length(x);for ii=1:NRxx(ii) = x(1:N-ii+1)*(x(ii:N))'/N;enda(1)=1;a(2)=-Rxx(2)/Rxx(1);for k=1:p-1 % Levinson-Durbin algorithmvar(k+1) = Rxx(0+1)+a(1+1:k+1)*Rxx(1+1:k+1)';reflect_coefficient(k+1+1) = -a(0+1:k+1)*(fliplr(Rxx(2:k+1+1)))'/var(k+1);var(k+1+1) = (1-(reflect_coefficient(k+1+1))^2)*var(k+1);a_temp(1) = 1;for kk=1:ka_temp(kk+1) = a(kk+1)+reflect_coefficient(k+1+1)*a(k+1-kk+1);enda_temp(k+1+1) = reflect_coefficient(k+1+1);a = a_temp;enda_p = a; % prediction coeffecientsvar_p = var(p+1); % prediction error power1.4 仿真结果1）p=2时，仿真结果图如下预测系数：误差功率：var_p=1.01942）p=20时，仿真结果图如下预测系数：误差功率：var_p=0.99983）p=50时，仿真结果图如下预测系数：误差功率：var_p=0.99551.5 结果分析由不同阶数（P值）得到的仿真结果可得：当P的阶数较低时，L-D算法估计AR模型对功率谱估计的分辨率较低，有平滑的效果，从P=2的仿真结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；随着阶数增大，采用L-D算法进行估计后，得到的功率谱会产生振荡，从仿真可以看到，当阶数P较高为50时，估计得到的功率谱与原始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。
现代数字信号处理仿真作业
现代数字信号处理仿真作业1.仿真题3.17仿真结果及图形：图基于FFT的自相关函数计算图图周期图法和BT法估计信号的功率谱图利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为：a1=--;a2=;a3=3i;a4=7;a5=68i;a6=7+6i;a7=9-2i;a8=2-0ia9=2+0i;a10=2+3i;a11=7-10i;a12=4-9i;a13=8-3i ;a14=2+4i;a15=2+1i;a16=3i.仿真程序（3_17）:clear allclc%%产生噪声序列N=32;%基于FFT的样本长度%N=256;%周期图法，BT法，AR模型功率谱估计的长度vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);%%产生复正弦信号f=[0.150.170.26];%归一化频率SNR=[303027];%信噪比A=10.^(SNR./20);%幅度signal=[A(1)*exp(1i*2*pi*f(1)*(0:N-1));%复正弦信号A(2)*exp(1i*2*pi*f(2)*(0:N-1));A(3)*exp(1i*2*pi*f(3)*(0:N-1))];%%产生观察样本un=sum(signal)+vn;%%利用3.1.1的FFT估计Uk=fft(un,2*N);Sk=(1/N)*abs(Uk).^2;r0=ifft(Sk);r1=[r0(N+2:2*N),r0(1:N)];%%r2=xcorr(un,N-1,'biased');%画图k=-N+1:N-1;figure(1)subplot(1,2,1)stem(k,real(r1))xlabel('m');ylabel('实部');subplot(1,2,2)stem(k,imag(r1))xlabel('m');ylabel('虚部');figure(2)subplot(1,2,1)stem(k,real(r2))xlabel('m');ylabel('实部');subplot(1,2,2)stem(k,imag(r2))xlabel('m');ylabel('虚部');%%周期图法NF=1024;Spr=fftshift((1/NF)*abs(fft(un,NF)).^2);kk=-0.5+(0:NF-1)*(1/(NF-1));Spr_norm=10*log10(abs(Spr)/max(abs(Spr)));%%BT法M=64;r3=xcorr(un,M,'biased');BT=fftshift(fft(r3,NF));BT_norm=10*log10(abs(BT)/max(abs(BT)));figure(3)subplot(1,2,1)plot(kk,Spr_norm)xlabel('w/2pi');ylabel('归一化功率谱/DB');title('周期图法')subplot(1,2,2)plot(kk,BT_norm)xlabel('w/2pi');ylabel('归一化功率谱/DB');title('BT法')%%LD迭代算法p=16;r0=xcorr(un,p,'biased');r4=r0(p+1:2*p+1);%计算自相关函数a(1,1)=-r4(2)/r4(1);sigma(1)=r4(1)-(abs(r4(2))^2)/r4(1);for m=2:p%LD迭代算法k(m)=-(r4(m+1)+sum(a(m-1,1:m-1).*r4(m:-1:2)))/sigma(m-1);a(m,m)=k(m);for i=1:m-1a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));endsigma(m)=sigma(m-1)*(1-abs(k(m))^2);endPar=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);%p阶AR模型的功率谱Par_norm=10*log10(abs(Par)/max(abs(Par)));figure(4)plot(kk,Par_norm)xlabel('w/2pi');ylabel('归一化功率谱/DB');title('16阶AR模型')2.仿真题3.20仿真结果及图形：单次Root-MUSIC算法中最接近单位圆的两个根为：+-对应的归一化频率为：相同信号的MUSIC谱估计结果如下图对3.20信号进行MUSIC谱估计的结果仿真程序（3_20）:clear allclc%%信号样本和高斯白噪声的产生N=1000;vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal=[exp(1i*0.5*pi*(0:N-1)+1i*2*pi*rand);%复正弦信号exp(-1i*0.3*pi*(0:N-1)+1i*2*pi*rand)];un=sum(signal)+vn;%%计算自相关矩阵M=8;for k=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=xs*xs'/(N-M);%%自相关矩阵的特征值分解[U,E]=svd(R);%%Root-MUSIC算法的实现G=U(:,3:M);Gr=G*G';co=zeros(2*M-1,1);for m=1:Mco(m:m+M-1)=co(m:m+M-1)+Gr(M:-1:1,m);endz=roots(co);ph=angle(z)/(2*pi);err=abs(abs(z)-1);%%计算MUSIC谱En=U(:,2+1:M);NF=2048;for n=1:NFAq=exp(-1i*2*pi*(-0.5+(n-1)/(NF-1))*(0:M-1)');Pmusic(n)=1/(Aq'*En*En'*Aq);endkk=-0.5+(0:NF-1)*(1/(NF-1));Pmusic_norm=10*log10(abs(Pmusic)/max(abs(Pmusic))); plot(kk,Pmusic_norm)xlabel('w/2*pi');ylabel('归一化功率谱/dB')3.仿真题3.21仿真结果及图形：单次ESPRIT算法中最接近单位元的两个特征值为：+-对应的归一化频率为：仿真程序（3_21）:clear allclc%%信号样本和高斯白噪声的产生N=1000;vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal=[exp(1i*0.5*pi*(0:N-1)+1i*2*pi*rand);%复正弦信号exp(-1i*0.3*pi*(0:N-1)+1i*2*pi*rand)];un=sum(signal)+vn;%%自相关矩阵的计算M=8;for k=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endRxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-M-1);Rxy=xs(:,1:end-1)*xs(:,2:end)'/(N-M-1);%%特征值分解[U,E]=svd(Rxx);ev=diag(E);emin=ev(end);Z=[zeros(M-1,1),eye(M-1);0,zeros(1,M-1)];Cxx=Rxx-emin*eye(M);Cxy=Rxy-emin*Z;%%广义特征值分解[U,E]=eig(Cxx,Cxy);z=diag(E);ph=angle(z)/(2*pi);err=abs(abs(z)-1);4.仿真题4.18仿真结果及图形：步长为0.05时失调参数为m1=0.0493；步长为0.005时失调参数为m2=0.0047。
现代信号处理及其应用仿真
仿真作业周子超200820003043 3.17(1)clear allclose allf1=0.15;f2=0.17;f3=0.26;v=randn(1,256);A1=10^(30/20);A2=10^(30/20);A3=10^(30/27);t=0:255;x1=A1*exp(j*2*pi*f1*t);x2=A2*exp(j*2*pi*f2*t);x3=A3*exp(j*2*pi*f3*t);u=v+x1+x2+x3;u=u(:,1:32);u1=[u,zeros(1,32)];uk=fft(u1);for i=1:64uuk(i)=(1/32)*uk(i)*conj(uk(i));endr0=ifft(uuk);rr1=[r0(1,34:64),r0(1,1:32)];r=0;r10=0;r1=zeros(1,63);u=[u,zeros(1,32)];for n=1:32r=(1/32)*u(n)*conj(u(n));r10=r+r10;endfor m=1:31for n=m+1:32r=(1/32)*u(n)*conj(u(n-m));r1(1,m)=r1(1,m)+r;endendfor m=32:62for n=1:31r=(1/32)*u(n)*conj(u(n+m-31));r1(1,m)=r1(1,m)+r;endendrr2=[r1(1,62:-1:32),r10,r1(1,1:31)];通过上面的计算可以发现基于FFT的自相关函数快速算法估计的自相关函数和式(3.1.2)估计出的自相关函数相等。
(2)clear allclose allf1=0.15;f2=0.17;f3=0.26;v=randn(1,256);A1=10^(30/20);A2=10^(30/20);A3=10^(27/20);t=0:255;x1=A1*exp(j*2*pi*f1*t);x2=A2*exp(j*2*pi*f2*t);x3=A3*exp(j*2*pi*f3*t);u=v+x1+x2+x3;U=fft(u);for i=1:256SPER(i)=abs(U(i)*conj(U(i)))/256;endSPER=10*log(SPER);n=0:length(u)-1;k=n/(length(u)-1);figureplot(k,SPER)xlabel('归一化频率')ylabel('功率谱(dB)')title('周期图法')u=u(:,1:256);r=0;r10=0;r1=zeros(1,511);u=[u,zeros(1,256)];for n=1:256r=(1/256)*u(n)*conj(u(n));r10=r+r10;endfor m=1:255for n=m+1:256r=(1/256)*u(n)*conj(u(n-m));r1(1,m)=r1(1,m)+r;endendfor m=256:510for n=1:255r=(1/256)*u(n)*conj(u(n+m-255));r1(1,m)=r1(1,m)+r;endendrr=[r1(1,510:-1:256),r10,r1(1,1:255)];rr=rr(1,193:319);k=-63:63;w=pi*k/63+pi;Sbtw=rr*(exp(-j)).^(k'*w);SBTW=abs(Sbtw);P=10*log(SBTW);figureplot(w/(2*pi),P)xlabel('归一化频率')ylabel('功率谱(dB)')title('BT 法')00.10.20.30.40.50.60.70.80.91-40-2020406080100120140归一化频率功率谱（d B ）周期图法00.10.20.30.40.50.60.70.80.9160708090100110120归一化频率功率谱（d B ）BT 法(3)clear allclose allf1=0.15;f2=0.17;f3=0.26;v=randn(1,256);A1=10^(30/20);A2=10^(30/20);A3=10^(27/20);t=0:255;x1=A1*exp(j*2*pi*f1*t);x2=A2*exp(j*2*pi*f2*t);x3=A3*exp(j*2*pi*f3*t);u=v+x1+x2+x3;r=0;r10=0;r1=zeros(1,511);u=[u,zeros(1,256)];for n=1:256r=(1/256)*u(n)*conj(u(n));r10=r+r10;endfor m=1:255for n=m+1:256r=(1/256)*u(n)*conj(u(n-m));r1(1,m)=r1(1,m)+r;endendr00=0;r11=0;a=zeros(16,16);delta=zeros(16);r=[r10,r1(1,1:16)];a(1,1)=-r(2)/r(1);delta(1)=r(1)-r(2)*conj(r(2))/r(1);for m=2:16r11=0;for i=1:m-1r00=a(i,m-1)*r(m+1-i);r11=r11+r00;endk=-(r(m+1)+r11)/delta(m-1);a(m,m)=k;for i=1:m-1a(i,m)=a(i,m-1)+k*conj(a(m-i,m-1));enddelta(m)=delta(m-1)*(1-abs(k)^2);endap=a(:,16).';delta=delta(16);ap=[1,ap,zeros(1,1007)];i=1;for w=0:2*pi/1023:2*pie=exp(-j*w.*(0:1023));AP=ap*e.';SARW(i)=delta/((1+AP)*conj(1+AP)); i=i+1;endP=abs(10*log10(SARW));w=0:1/1023:1;figureplot(w,P)xlabel('归一化频率')ylabel('功率谱(dB)')title('Levinson-Durbin迭代算法')00.10.20.30.40.50.60.70.80.91归一化频率功率谱（d B ）Lev inson-Durbin 迭代算法3.20clear allclose allM=8;N=1000;n=1:N;x=exp(j*2*pi*0.25*n)+exp(j*2*pi*0.6*n)+randn(size(n));Rxx=zeros(M,M);for n=M:NX=zeros(1,M);X=x(n:-1:n-M+1);Rxx=Rxx+X'*X;endRxx=(1/(N-M+1))*Rxx;[V D]=eig(Rxx,'nobalance');weights=diag(D);weights=weights(1:6);for k=1:6;noise_eigenvects(:,k)=V(:,k);end% root-musicP=0;for i=1:length(weights),P=P+conv(noise_eigenvects(:,i),conj(flipud(noise_eigenvects(:,i)))); endroots_P=roots(P);roots_P1=roots_P(abs(roots_P)<1);[not_used,indx]=sort(abs(abs(roots_P1)-1));sorted_roots=roots_P1(indx);if(angle(sorted_roots(1))>0)f1rootmusic=angle(sorted_roots(1))/(2*pi);elsef1rootmusic=(angle(sorted_roots(1))+2*pi)/(2*pi);endif(angle(sorted_roots(2))>0)f2rootmusic=angle(sorted_roots(2))/(2*pi);elsef2rootmusic=(angle(sorted_roots(2))+2*pi)/(2*pi);endfrootmusic=[f1rootmusic f2rootmusic];% musici1=1;for w=0:2*pi/1023:2*piaw=[1 exp(j*w) exp(j*2*w) exp(j*3*w) exp(j*4*w) exp(j*5*w) exp(j*6*w) exp(j*7*w)].';Pmusic(i1)=1/abs(aw'*noise_eigenvects*noise_eigenvects'*aw);Pmusic(i1)=10*log10(Pmusic(i1));i1=i1+1;endw=0:1/1023:1;plot(w,Pmusic)xlabel('归一化频率')ylabel('功率谱(dB)')title('Music算法')通过上面的计算，Root-Music算法估计的频率为f1=0.25，f2=0.5999，而Music 算法的估计结果为f1=0.248，f2=0.63，所以Root-Music算法估计的频率较Music 算法估计的准确。
现代信号处理仿真作业一(318谐波恢复)
现代信号处理仿真作业一(3.18谐波恢复)班级：自研42 姓名：李琳琳 学号：2004211068一． 谐波恢复的基本理论和方法：1． Pisarenko 谐波分解理论谐波过程可用差分方程描述，首先利用Pisarenko 谐波分解理论推导谐波过程所对应的差分方程。
对单个正弦波()sin(2)x n fn πθ=+利用三角函数恒等式，有：()2cos(2)(1)(2)0x n f x n x n π--+-=对上式作z 变换，得： 12[12cos(2)]()0f z z X z π---+=得到特征多项式：1212cos(2)0f z z π---+=。
由此可见，正弦波的频率可以由相应特征方程的一对共轭根来决定：|arctan[Im()/Re()]|/2i i i f z z π=将单个正弦波推广到多个正弦波的情形，得：如果p 个实的正弦波信号没有重复频率的话，则这p 个频率应该由特征多项式1(21)212110p p p a z a z z -----++++= (1)的根决定。
由此可得到p 个实正弦波所组成的谐波过程可以用以下的差分方程进行描述：21()()0pi i x n a xn i =+-=∑这是一个无激励的AR 过程。
2. 谐波恢复的ARMA 建模法在无激励的AR 模型差分方程21()()0pi i x n a x n i =+-=∑两边同乘()x n k -，并取数学期望，则有：21()()0,px i x i R k a R k i k =+-=∀∑ (2)正弦波过程一般是在加性白噪声中被观测的，设加性白噪声为()w n ，即观测过程为： 1()()()sin(2)()pi i i i y n x n w n A f n w n πθ==+=++∑，其中()w n 为0均值的高斯白噪声。
由于谐波过程和加性白噪声统计独立，所以有：2()()()()()y x w x w R k R k R k R k k σδ=+=+ (3)将(3)式代入(2)式得： 21()()0,2py i y i R k a R k i k p =+-=>∑..................... . (4)这就是加性白噪声中观测到的p 个正弦信号所组成的谐波信号的ARMA 过程的所服从的法方程。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
现代数字信号处理仿真作业第三章仿真作业3.17（1）代码clear;N=32;m=[-N+1:N-1];noise=(randn(1,N)+j*randn(1,N))/sqrt(2); f1=0.15;f2=0.17;f3=0.26;SNR1=30;SNR2=30;SNR3=27;A1=10^(SNR1/20);A2=10^(SNR2/20);A3=10^(SNR3/20);signal1=A1*exp(j*2*pi*f1*(0:N-1));signal2=A2*exp(j*2*pi*f2*(0:N-1));signal3=A3*exp(j*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;uk=fft(un,2*N);sk=(1/N) *abs(uk).^2;r0=ifft(sk);r1=[r0(N+2:2*N),r0(1:N)];r=xcorr(un,N-1,'biased');figuresubplot(2,2,1)stem(m,real(r1));xlabel('m');ylabel('FFT估计r1实部');subplot(2,2,2)stem(m,imag(r1));xlabel('m');ylabel('FFT估计r1虚部');subplot(2,2,3)stem(m,real(r));xlabel('m');ylabel('平均估计r实部');subplot(2,2,4)stem(m,imag(r));xlabel('m');ylabel('平均估计r虚部');仿真结果（2）代码 clear; N=256;noise=(randn(1,N)+j*randn(1,N))/sqrt(2); f1=0.15; f2=0.17; f3=0.26; SNR1=30; SNR2=30; SNR3=27;A1=10^(SNR1/20); A2=10^(SNR2/20); A3=10^(SNR3/20);signal1=A1*exp(j*2*pi*f1*(0:N-1)); signal2=A2*exp(j*2*pi*f2*(0:N-1)); signal3=A3*exp(j*2*pi*f3*(0:N-1)); un=signal1+signal2+signal3+noise;-40-2002040-200020004000mF F T 估计r 1实部-40-2002040-2000-1000010002000mF F T 估计r 1虚部-40-2002040-200020004000m平均估计r 实部-40-2002040-2000-1000010002000m平均估计r 虚部spr=fftshift((1/NF)*abs(fft(un,NF)).^2);f1=(0:length(spr)-1)*(1/(length(spr)-1))-0.5;M=64;r=xcorr(un,M,'biased');bt=fftshift(abs(fft(r,NF)));f2=(0:length(bt)-1)*(1/(length(bt)-1))-0.5;figuresubplot(1,2,1)plot(f1,10*log10(spr/max(spr)));xlabel('w/2pi');仿真结果w/2pi w/2pi（3）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果-4-0.5-0.4-0.3-0.2-0.100.10.20.30.40.53.20（1）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);nw=128;ww=[-128:128]/128*pi; e=exp(-j*ww'*[0:p-1])%k 行m 列 ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果距离单位圆最近的两个解为-0.2363-0.9717i 和0.3747+0.9271i ，对应的归一化频率为0.1889和-0.2880 （2）代码 clear; N=1000;fai1=rand(1,1)*2*pi; fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;cx=xcorr(un,p,'biased'); rxx=cx(p+1:2*p)'; R=toeplitz(rxx); [u,s]=eig(R); nw=128;ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k 行m 列 ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw)); 仿真结果-2-10123456-33.21代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;for k=1:N-pxs(:,k)=un(k+p-1:-1:k)';endrxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-p-1);rxy=xs(:,1:end-1)*xs(:,2:end)'/(N-p-1);[u,e]=svd(rxx);ev=diag(e);emin=ev(end);z=[zeros(p-1,1),eye(p-1);0,zeros(1,p-1)];cxx=rxx-emin*eye(p);cxy=rxy-emin*z;[U,E]=eig(cxx,cxy);Z=diag(E);ph=angle(Z)/(2*pi);err=abs(abs(Z)-1);仿真结果最接近单位圆的两个解分别为0.5867+0.8097i和0.0349-0.9984i，对应的归一化频率为0.1502和-0.2444。
第四章仿真题作业4-18（1）产生随机序列代码clear;data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;vn=sqrt(sigma_v_2)*randn(data_len,1,trials);u0=[0 0];num=1;den=[1 a1 a2];zi=filtic(num,den,u0);u=filter(num,den,vn,zi);(2)单次实验估计最优权向量代码clear;data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;vn=sqrt(sigma_v_2)*randn(data_len,1,trials);u0=[0 0];num=1;m=1;den=[1 a1 a2];zi=filtic(num,den,u0);u=filter(num,den,vn,zi);mu1=0.05;mu2=0.005;w1=zeros(2,data_len);w2=zeros(2,data_len);e1=zeros(data_len,1);e2=zeros(data_len,1);d1=zeros(data_len,1);d2=zeros(data_len,1);%for m=1:trialsfor n=3:data_len-1w1(:,n+1)=w1(:,n)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n));w2(:,n+1)=w2(:,n)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n)); d1(n+1)=w1(:,n+1)'*u(n:-1:n-1,:,m); d2(n+1)=w2(:,n+1)'*u(n:-1:n-1,:,m); e1(n+1)=u(n+1,:,m)-d1(n+1); e2(n+1)=u(n+1,:,m)-d2(n+1); endfigure(1) hold on;plot(w1(1,:));plot(w1(2,:));plot(w2(1,:),’g ’); plot(w2(2,:),’g ’) 仿真结果（3）100次实验计算剩余均方误差 代码clear;data_len=512; trials=100; n=1:data_len; a1=-0.975; a2=0.95;sigma_v_2=0.0731;vn=sqrt(sigma_v_2)*randn(data_len,1,trials);0100200300400500600-1.5-1-0.50.511.5u0=[0 0];num=1;%m=1;den=[1 a1 a2];zi=filtic(num,den,u0);u=filter(num,den,vn,zi);mu1=0.05;mu2=0.005;w1=zeros(2,data_len);w2=zeros(2,data_len);e1=zeros(data_len,1);e2=zeros(data_len,1);d1=zeros(data_len,1);d2=zeros(data_len,1);jn1=zeros(1,trials);jn2=zeros(1,trials);mse1=zeros(1,trials);mse2=zeros(1,trials);for m=1:trialsfor n=3:data_len-1w1(:,n+1)=w1(:,n)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n)); w2(:,n+1)=w2(:,n)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n)); d1(n+1)=w1(:,n+1)'*u(n:-1:n-1,:,m);d2(n+1)=w2(:,n+1)'*u(n:-1:n-1,:,m);e1(n+1)=u(n+1,:,m)-d1(n+1);e2(n+1)=u(n+1,:,m)-d2(n+1);endjn1(m)=e1(256);jn2(m)=e2(256)mse1(m)=abs(sum(jn1(1:m))/m);mse2(m)=abs(sum(jn2(1:m))/m);endfigure(1)hold on;plot (mse1,'g');plot(mse2,'b');wopt=zeros(2,trials);jmin=zeros(1,trials);sum_eig=zeros(trials,1);for m=1:trialsrm=xcorr(u(:,:,m),'biased');r=[rm(512),rm(513);rm(511),rm(512)];p=[rm(511);rm(510)]wopt(:,m)=r\p;[v,d]=eig(r);jmin(m)=rm(512)-p'*wopt(:,m);sum_eig(m)=d(1,1)+d(2,2); endfigure(2)plot(jmin,'r');sjmin=sum(jmin)/trials;e1_100trials_ave=sum(jn1)/trials; e2_100trials_ave=sum(jn2)/trials; jex1=e1_100trials_ave-sjmin; jex2=e2_100trials_ave-sjmin;sum_eig_100trials=sum(sum_eig)/100;jexfin=mu1*sjmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials)); jexfin2=mu2*sjmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials)); m1=jexfin/sjmin; m2=jexfin2/sjmin; 仿真结果LMS 算法学习曲线剩余均方误差jex1=-0.0569,jex2=-0.0959 失调参数：m1=0.0523,m2=0.005010203040506070809010000.050.10.150.20.250.30.350.40.45迭代次数M S Eu=0.05u=0.005第五章仿真题1、当M=2时构造yule_walk方程=得到故自相关矩阵2、当M=3时构造Yule_walk方程=得到故自相关矩阵为3、r1=[47.0487 46.578346.5783 47.0487]eig_value=eig(r);eig_spred=max(eig_value)/min(eig_value);r2=[47.0487 46.5783 46.112546.5783 47.0487 46.578346.1125 46.5783 47.0487]eig_value=eig(r2);eig_spred=max(eig_value)/min(eig_value);特征值扩展结果=199.0370,=444.4107 4、代码clear;data_len=10000;trials=100;n=1:data_len;a1=-0.99;M=2;wsum=zeros(data_len,M);esum=zeros(1,data_len);sigma_v_2=0.93627;for m=1:trialsvn=sqrt(sigma_v_2)*randn(data_len,1);u(1)=vn(1);for k=2:data_lenu(k)=-a1*u(k-1)+vn(k)endw(1,:)=zeros(1:M);e(1)=u(1);mu=0.001;uu=zeros(1,M);w(2,:)=w(1,:)+mu*e(1)*uu;uu=[u(1) uu(1:M-1)];dd=(w(2,:)*uu')';e(2)=u(3)-dd;for k=3:data_lenw(k,:)=w(k-1,:)+mu*e(k-1)*uu;uu=[u(k-1) uu(1:M-1)];dd=(w(k,:)*uu')';e(k)=u(k)-dd;endwsum=wsum+w;esum=esum+abs(e);endfigure(1)hold on;plot(wsum/100);plot(w);figure(2)plot(esum/100);仿真结果w迭代次数n迭代次数nJ （n ）第六章仿真作业6-13 代码clear;m=4;n=1000;f1=[0.1 0.25 0.27];SNR=[30 30 27];sigma=1;am=sqrt(sigma*10.^(SNR/10));t=linspace(0,1,n);pha=2*pi*rand(size(f1));vn=sqrt(sigma/2)*randn(size(t))+j*sqrt(si gma/2)*randn(size(t));un=vn;for k=1:length(f1)s=am(k)*exp(j*2*pi*n*f1(k).*t+j*pha(k)); un=un+s;end %产生观测信号A=zeros(m,n-m+1);for i=1:n-m+1A(:,i)=un(m+i-1:-1:i);end[U,S,V]=svd(A');invpha=V*inv(S'*S)*V';p=1024;f2=linspace(-0.5,0.5,p);omega=2*pi*f2;a=zeros(m,p);for k=1:pfor i=1:ma(i,k)=exp(-j*omega(k)*(i-1));endendpmvdr=zeros(size(omega));for k=1:ppmvdr(k)=1/(a(:,k)'*invpha*a(:,k)); endpmvdr=abs(pmvdr/max(abs(pmvdr)));plot(omega/(2*pi),10*log10(pmvdr));仿真结果6-15a1=0.99; sigma=0.995; N=1000; trials=500; n0=1; M=2; L=N-n0;vn=sqrt(sigma)*randn(N,1,trials); nume=1; deno=[1 a1];w/(2*pi)归一化的功率谱密度/d Bu0=zeros(length(deno)-1,1);w=zeros(M,L+1,trials);epsilon=zeros(L,1,trials);W_mean=zeros(M,L+1);MSE=zeros(L,1);for m=1:trialsxic=filtic(nume,deno,u0);un=filter(nume,deno,vn(:,1,m),xic); b=un(n0+1:N);un1=[zeros(M-1,1).',un.'].';A=zeros(M,L);for k=1:LA(:,k )= un1(M-1+k:-1:k);enddelta=0.004;lambda=0.98;P1=eye(M)/delta;for k=1:LPIn=P1*A(:,k);denok=lambda+A(:,k)'*PIn;kn=PIn/denok;epsilon(k,1,m)=b(k)-w(:,k,m)'*A(:,k);w(:,k+1,m)=w(:,k,m)+kn*conj(epsilon(k,1,m ));P1=P1/lambda-kn*A(:,k)'*P1/lambda;endMSE=MSE+abs(epsilon(:,1,m)).^2;W_mean=W_mean+w(:,:,m);endfigure(1);plot(1:L,1/trials*MSE);xlabel('迭代次数');ylabel('MSE');figure(2);plot(1:L+1,1/trials*W_mean); xlabel('迭代次数');ylabel('权值');第七章仿真作业7-14代码clear all;N=3000;gv=0.0332;trials=100;v=randn(N,1,trials)*sqrt(gv);a1=-2;a2=1.4;a3=-0.4;a4=0.0384;u1=zeros(N,1,trials);for m=1:trialsfor i=1:N-4u1(i+4,1,m)=-a1*u1(i+3,1,m)-a2*u1(i+2,1,m )-a3*u1(i+1,1,m)...-a4*u1(i,1,m)+v(i+4,1,m);endendN2=2000;Jmin=0.005;W_esti=zeros(4,N2,trials);W=zeros(4,N2+1);e=zeros(N2,1);for m=1:trialsP_esti=eye(4);for n=5:N2P_pre=P_esti;U(:,n,m)=[u1(n-1,1,m);u1(n-2,1,m);u1(n-3, 1,m);u1(n-4,1,m)];A=(U(:,n,m))'*P_pre*U(:,n,m)+Jmin; K=P_pre*U(:,n,m)/A;alpha(n,m)=u1(n,1,m)-(U(:,n,m))'*W_esti(: ,n,m);W_esti(:,n+1,m)=W_esti(:,n,m)+K*alpha(n,m );P_esti=P_pre-K*(U(:,n,m))'*P_pre;endW=W+W_esti(:,:,m);e=e+abs(alpha(:,m)).^2;endW=1/trials*W;e=1/trials*e;h=1:N2+1;figure(1);plot(h,W);xlabel('迭代次数'); ylabel('权值'); gtext('w1');gtext('w2');gtext('w3');gtext('w4');figure(2);plot(1:N2,e);xlabel('迭代次数'); ylabel('MSE');第八章仿真题clc;clear;N=100;M=10;K=2;theta=[-10;40]*pi/180;SNR=[10;20];Am=[sqrt(10.^(SNR/10))];S=Am*ones(1,N);S(2,:)=S(2,:).*exp(j*2*pi*rand(1,N)); for a=1:Mfor b=1:KA(a,b)=exp(-j*(a-1)*pi*sin(theta(b))); endendV=zeros(M,N);for m=1:Mv=wgn(1,N,0,'complex');v=v-mean(v);v=v/std(v);V(m,:)=v;endX=A*S+V;R=zeros(M,M);for i=1:NR=R+X(:,i)*X(:,i)';endR=R/N;%%%%%%%%%%%%MUSIC算法%%%%%%%%%%%%%% [VR,D]=eig(R);[B,IX]=sort(diag(D));G=VR(:,IX(M-K:-1:1));P=[];for n=-pi/2:pi/180:pi/2a=exp(-j*[0:M-1]'*pi*sin(n));P=[P,1/(a'*G*G'*a)];endt=-90:1:90;figure(1);plot(t,10*log10(abs(P)/max(abs(P))))title('MUSIC算法')ylabel('归一化功率谱/dB')xlabel('空间角度')% %%% %%%%%%%%%%%%RootMUSIC算法%%%%%%%%%%%%%%% syms z% pz=z.^([0:M-1]');% pz1=(z^(-1)).^([0:M-1]);% fz=z^(M-1)*pz1*G*G'*pz;% a=sym2poly(fz);% r=roots(a);% r1=abs(r);% for i=1:2*K% [Y,I(i)]=min(abs(r1-1));% r1(I(i))=inf;% end% for i=1:2*K%theta_esti1(i)=asin(-angle(r(I(i)))/pi)*1 80/pi;% end% %%% % %%%%%%%%%%%%%%%ESPRIT算法%%%%%%%%%%%%%%%%% S=VR(:,IX(M:-1:M-K+1));% S1=S(1:M-1,:);% S2=S(2:M,:);% fai=S1\S2;%% [U_fai,V_fai]=eig(fai);% for i=1:K%theta_esti2(i)=asin(-angle(V_fai(i,i))/pi )*180/pi;% end% %%% %%%%%%%%%%%%%%%%%%MVDR算法%%%%%%%%%%%%%%%%% P=[];% for n=-pi/2:pi/180:pi/2% a=exp(-j*[0:M-1]'*pi*sin(n));% P=[P,1/(a'*inv(R)*a)];% end% P=abs(P/max(abs(P)));% P=10*log10(P);% figure(2);% plot(t,P);% title('MVDR算法')% ylabel('归一化功率谱/dB')% xlabel('空间角度')% %%% %%%%%%%%%%%%%%%%%%F-SAPES算法%%%%%%%%%%%%%%%%% P=6;% L=M+1-P;% Rf=zeros(L,L);% for i=1:P% Rf=Rf+X(i:i+L-1,:)*X(i:i+L-1,:)'/N; % end% Rf=Rf/P;% n1=0:P-1;% n2=0:L-1;% cc=[1 zeros(1,L-1)];% for n3=-90:.5:90% fy=exp(j*pi*sin(n3/180*pi));% tt=[(fy.^(n1')).' zeros(1,M-P)];% Tfy=toeplitz(cc,tt);% GfTheta=1./(P^2)*Tfy*R*Tfy';% Qf=Rf-Rf-GfTheta;% aTheta=fy.^(-n2');%Wof=((inv(Qf))*aTheta)./(aTheta'*inv(Qf)* aTheta);%sigma2Theta(((n3+90)/.5+1))=Wof'*GfTheta* Wof;% end%sigma2Theta=abs(sigma2Theta/max(abs(sigma 2Theta)));% sigma2Theta=10*log10(sigma2Theta);% figure(3);% t1=-90:.5:90;% plot(t1,sigma2Theta);% title('F-SAPES算法')% ylabel('归一化功率谱/dB')% xlabel('空间角度')-35-30-25-20-15-10-5MVDR 算法归一化功率谱/d B 空间角度-100-80-60-40-20020406080100-30-25-20-15-10-5F-SAPES 算法归一化功率谱/d B 空间角度。