现代数字信号处理仿真作业.

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,:));。

目录仿真一：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)。

数字信号处理matlab仿真数字信号处理作业设计报告一、目的1.增进对Matlab的认识，加深对数字信号处理理论方面的理解。

2.掌握数字信号处理中IIR和FIR滤波器的设计。

3.了解和掌握用Matlab实现IIR和FIR滤波器的设计方法、课程，为以后的设计打下良好基础。

二、设计内容1.IIR（无限脉冲响应）模拟滤波器设计（1）设计题目：椭圆型模拟带通IIR滤波器技术指标：通带下截止频率fpl=2kHz，上截止频率fph=5kHz，通带内最大衰减ap=1dB；阻带下截止频率fsl=1.5kHz上截止频率fsh=5.5kHz,阻带最小衰减as=40dB. 设计原理：①确定模拟带通滤波器的技术指标，并对边界频率做归一化处理；②确定归一化低通技术要求；③设计归一化低通G(p)；④将低通G(p)转换成带通H(s)。

Matlab原程序如下：clear all;fp=[2000,5000];ap=1;fs=[1500,5500];as=40;wp=2*pi*fp;ws=2*pi*fs;%归一化的截止频率[N,wn]=ellipord(wp,ws,ap,as, &#39;s&#39;);%求椭圆形滤波器的最小阶数和归一化截止频率[B,A]=ellip(N,ap,as,wn, &#39;s&#39;);%求传递函数的分子分母系数[H,w]=freqs(B,A);%频率响应函数f=0:8000;1w=2*pi*f;H=freqs(B,A,w);%求系统在指定频率点w上的频响Hplot(f,20*log10(abs(H)));%绘图显示axis([0 7000 -80 0])仿真波形图如下：（2）设计题目：巴特沃斯低通模拟滤波器技术指标：通带截止频率fp=5kHz，通带内最大衰减ap=2dB；阻带截止频率fs=12kHz，阻带最小衰减as=30dB。

设计原理：①确定模拟带通滤波器的技术指标，并对边界频率做归一化处理；②确定归一化低通技术要求并求出归一化低通原型系统函数Ga(p);③将Ga(p)去归一化。

姓名：潘晓丹学号：班级：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 算法估计的准确。

现代信号处理仿作业（谐波恢复）现代信号处理仿作业(谐波恢复)————————————————————————————————作者：————————————————————————————————日期：现代信号处理仿真作业一(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 -----++++=K (1)的根决定。

由此可得到p 个实正弦波所组成的谐波过程可以用以下的差分方程进行描述：21()()0pi i x n a x n 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均值的高斯白噪声。

数字信号处理上机实验仿真报告学 院 电子工程学院班 级学 号 姓 名实验一： 信号、 系统及系统响应一、实验目的(1) 熟悉连续信号经理想采样前后的频谱变化关系， 加深对时域采样定理的理解。

(2) 熟悉时域离散系统的时域特性。

(3) 利用卷积方法观察分析系统的时域特性。

(4) 掌握序列傅里叶变换的计算机实现方法， 利用序列的傅里叶变换对连续信号、 离散信号及系统响应进行频域分析。

二、实验原理采样是连续信号数字处理的第一个关键环节。

对一个连续信号()a x t 进行理想采样的过程可用(1.1)式表示。

()()()ˆa a xt x t p t =⋅ (1.1) 其中()t xa ˆ为()a x t 的理想采样，()p t 为周期冲激脉冲，即 ()()n p t t nT δ∞=-∞=-∑ (1.2) ()t xa ˆ的傅里叶变换()j a X Ω为 ()()s 1ˆj j j a a m X ΩX ΩkΩT ∞=-∞=-∑ (1.3) 将(1.2)式代入(1.1)式并进行傅里叶变换，()j e ΩnT a n x nT ∞-=-∞=∑ (1.4)式中的()a x nT 就是采样后得到的序列()x n ， 即()()a x n x nT =()x n 的傅里叶变换为()()j j e en n X x n ωω∞-=-∞=∑ (1.5)比较(1.5)和(1.4)可知()()j ˆj e a ΩT X ΩX ωω== (1.6)为了在数字计算机上观察分析各种序列的频域特性，通常对()j e X ω在[]0,2π上进行M 点采样来观察分析。

对长度为N 的有限长序列()x n ，有()()1j j 0e ek k N n n X x n ωω--==∑ (1.7) 其中2π,0,1,,1k k k M Mω==⋅⋅⋅- 一个时域离散线性时不变系统的输入/输出关系为()()()()()m y n x n h n x m h n m ∞=-∞=*=-∑ (1.8) 上述卷积运算也可以转到频域实现()()()j j j e e e Y X H ωωω= (1.9)三、 实验内容及步骤(1) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。

数字信号处理课程仿真大作业一、课程设计题目：基于MATLAB 的利用FFT进行频谱分析二、课程设计目的：1、加深对数字信号处理学习过的基本概念、基本理论和基本方法的理解和掌握；2、学会用MATLAB对信号进行分析和处理，进一步将知识融会贯通；3、加深对FFT原理的理解，学会应用FFT分析信号频谱；4、学会撰写课程设计报告，并应用数字信号处理的基本理论分析结果。

三、课程设计内容：1、自编程序得到一个方波信号（f=50Hz,幅值为1，-1，各半个周期），对其一个周期分别采样256点和1024点，利用基于Matlab 语言所编FFT程序做谐波分析，并与理论分析结果对照。

2、对三角波信号（可以由方波信号求导得到）重复作业一的各项要求。

3、对一、二信号叠加一个白噪声信号（均值为零，方差为0.2）所构成的随机信号用FFT进行频谱分析。

4、对以上结果进行讨论。

5、给出源程序，包含信号如何产生、采样的实现、FFT函数的调用（或自编）、绘图等，给出计算机分析结果的图形截图及理论分析的草图。

四、详细程序及仿真波形：理论分析：用FFT分析周期函数的频谱结构，选择不同的截取长度Ts观察用FFT进行频谱分析时存在的截断效应(频谱泄露和谱间干扰)。

利用matlab的FFT对模拟信号进行谱分析时，只能以有限大的采样频率Fs对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作FFT变换，得到模拟信号的近似频谱。

其误差主要来自以下因素：(1)截断效应（频谱泄露和谱间干扰）截断效应使谱分辨率（能分辨开的两根谱线间的最小间距）降低，并产生谱间干扰。

(2)频谱混叠失真使折叠频率（Fs/2）附近的频谱产生较大失真。

理论和实践都已证明，加大截取长度Ts可以提高分辨率；另外选择合适的窗函数可降低谱间干扰；而频谱混叠失真要通过提高采样频率Fs或预滤波（在采样之前滤除折叠频率以外的频率成分）来改善。

用FFT 进行频谱谐波理论理论分析：在信号处理中，信号的频谱分析是重要的应用领域之一。

实验一 数字信号处理的Matlab 仿真一、实验目的1、掌握连续信号及其MATLAB 实现方法；2、掌握离散信号及其MA TLAB 实现方法3、掌握离散信号的基本运算方法，以及MATLAB 实现4、了解离散傅里叶变换的MA TLAB 实现5、了解IIR 数字滤波器设计6、了解FIR 数字滤波器设计1二、实验设备计算机，Matlab 软件三、实验内容（一）、 连续信号及其MATLAB 实现1、 单位冲击信号()0,0()1,0t t t dt εεδδε-⎧=≠⎪⎨=∀>⎪⎩⎰例1.1：t=1/A=50时，单位脉冲序列的MA TLAB 实现程序如下：clear all;t1=-0.5:0.001:0;A=50;A1=1/A;n1=length(t1);u1=zeros(1,n1);t2=0:0.001:A1;t0=0;u2=A*stepfun(t2,t0);t3=A1:0.001:1;n3=length(t3);u3=zeros(1,n3);t=[t1 t2 t3];u=[u1 u2 u3];plot(t,u)axis([-0.5 1 0 A+2])2、 任意函数()()()f t f t d τδττ+∞-∞=-⎰例1.2：用MA TLAB 画出如下表达式的脉冲序列()0.4(2)0.8(1) 1.2() 1.5(1) 1.0(2)0.7(3)f n n n n n n n δδδδδδ=-+-+++++++clear all;t=-2:1:3;N=length(t);x=zeros(1,N);x(1)=0.4;x(2)=0.8x(3)=1.2;x(4)=1.5;x(5)=1.0;x(6)=0.7;stem(t,x);axis([-2.2 3.2 0 1.7])3、 单位阶跃函数1,0()0,0t u t t ⎧≥⎪=⎨<⎪⎩例1.3：用MA TLAB 实现单位阶跃函数clear all;t=-0.5:0.001:1;t0=0;u=stepfun(t,t0);plot(t,u)axis([-0.5 1 -0.2 1.2])4、 斜坡函数0()()g t B t t =-例1.4：用MA TLAB 实现g(t)=3(t-1)clear all;t=0:0.01:3;B=3;t0=1;u=stepfun(t,t0);n=length(t);for i=1:nu(i)=B*u(i)*(t(i)-t0);endplot(t,u)axis([-0.2 3.1 -0.2 6.2])5、 实指数函数()at f t Ae =例1.5：用MA TLAB 实现0.5()3t f t e =clear all;t=0:0.001:3;A=3;a=0.5;u=A*exp(a*t);plot(t,u)axis([-0.2 3.1 -0.2 14])6、 正弦函数2()cos()t f t A T πϕ=+ 例1.6：用MA TLAB 实现正弦函数f(t)=3cos(10πt+1)clear all;t=-0.5:0.001:1;A=3;f=5;fai=1;u=A*sin(2*pi*f*t+fai);plot(t,u)axis([-0.5 1 -3.2 3.2])（二）、离散信号及其MATLAB 实现1、 单位冲激序列1,0()0,0n n n δ⎧=⎪=⎨≠⎪⎩例2.1：用MA TLAB 产生64点的单位冲激序列clear all;N=64;x=zeros(1,N);x(1)=1;xn=0:N-1;stem(xn,x)axis([-1 65 0 1.1])2、 任意序列()()()m f n f m n m δ∞=-∞=-∑例2.2：用MA TLAB 画出如下表达式的脉冲序列()8.0() 3.4(1) 1.8() 5.6(3) 2.9(4)0.7(5)f n n n n n n n δδδδδδ=+-++-+-+-clear all;N=8;x=zeros(1,N);x(1)=8.0;x(2)=3.4x(3)=1.8;x(4)=5.6;x(5)=2.9;x(6)=0.7;xn=0:N-1;stem(xn,x)axis([-1 8 0 8.2])3、 单位阶跃序列1,0()0,0n u n n ⎧≥⎪=⎨<⎪⎩例2.3：用MA TLAB 实现单位阶跃函数clear all;N=32;x=ones(1,N);xn=0:N-1;stem(xn,x)axis([-1 32 0 1.1])4、 斜坡序列0()()g n B n n =-例2.4：用MA TLAB 实现g(n)=3(n-4)点数为32的斜坡序列clear all;N=32;k=4B=3;t0=1;x=[zeros(1,k) ones(1,N-k)];for i=1:Nx(i)=B*x(i)*(i-k);endxn=0:N-1;stem(xn,x)axis([-1 32 0 90])5、 正弦序列()sin(2)x n A fn πϕ=+例2.5：用MA TLAB 实现幅度A=3，频率f=100，初始相位Φ=1.2，点数为32的正弦信号 clear all;N=32;A=3;f=100;fai=1.2;xn=0:N-1;x=A*sin(2*pi*f*(xn/N)+fai);stem(xn,x)axis([-1 32 -3.2 3.2])6、 实指数序列()n x n Aa =例2.6：用MA TLAB 实现0.7()3x n e =，点数为32的实指数序列clear all;N=32;A=3;a=0.7;xn=0:N-1;x=A*a.^xn;stem(xn,x)7、 复指数序列()(),a j n x n Ae n ω+=∀例2.7：用MA TLAB 实现幅度A=3，a=0.7，角频率ω=314，点数为32的实指数序列 clear all;N=32;A=3;a=0.7;w=314;xn=0:N-1;x=A*exp((a+j*w)*xn);stem(xn,x)8、 随机序列利用MATLAB 产生两种随机信号：rand(1,N)在区间上产生N 点均匀分布的随机序列randn(1,N)产生均值为0，方差为1的高斯随机序列，即白噪声序列例2.8：用MA TLAB 产生点数为32的均匀分布的随机序列与高斯随机序列clear all;N=32;x_rand=rand(1,N);x_randn=randn(1,N);xn=0:N-1;figure(1);stem(xn,x_rand)figure(2);stem(xn,x_randn)（三）、离散信号的基本运算1、 信号的延迟给定离散信号x(n)，若信号y(n)定义为：y(n)=x(n-k)，那么y(n)是信号x(n)在时间轴上右移k 个抽样周期得到的新序列。

第三章仿真作业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'); 仿真结果（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);w/2piw/2piww=[-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。

现代数字信号处理Matlab仿真报告组员：提交时间：一、算法介绍1、LMS算法：LMS算法采用平方误差最小的原则代替最小均方误差最小的原则，信号基本关系：式中，W(n)：n 时刻自适应滤波器的权矢量，，X( n) ：n 时刻自适应滤波器的输入，由最近N 个信号采样值构成，；N：自适应滤波器的阶数；d ( n) ：期望的输出值；e ( n) ：自适应滤波器的输出误差调节信号；μ：收敛因子。

2、RLS算法：RLS算法是用二乘方的时间平均的最小化准则取代最小均方准则，并按时间进行迭代计算。

其基本原理如下所示：称为遗忘因子，它是小于等于1的正数。

第n次迭代的权值。

均方误差。

按照如下准则：越旧的数据对的影响越小。

对滤波器系数求偏导数，并令结果等于零知整理得到标准方程定义标准方程可以化简成形式：经求解可以得到迭代形式定义：，则可知T的迭代方程为系数的迭代方程为其中增益和误差的定义分别为由上边分析可知，RLS算法递推的步骤如下：(1)在时刻n，已经知道和也已经存储在录波器的实验部件中(2)计算，并得到滤波器的输出相应和误差即：(3)进入第次迭代二、模型分析1、AR(2)模型因a1=1.4、a2=-0.7 ，即得：w(n)可用Matlab中高斯白噪声生成函数wgn生成其中在程序中用A描述，随时间变化2、LMS算法算出y(n);；；LMS函数源代码：function [A,en]=LMS(xn,Wn,M,u,i)en = zeros(i,1);A = zeros (i,M);for j = M+1:iyn(j)= A(j-1,:)*xn(j-1:-1:j-M)';en(j-1)=xn(j)-yn(j);x=xn(j-1:-1:j-M);A(j,:)=A(j-1,:)+2*u*en(j-1)*x;end3、RLS算法对赋一个比较大的初值，程序中=100*eye(M,M)；求出、；求出最新时刻T(n)；；RLS函数源代码：function [A,en]=RLS(xn,Wn,M,u,i)en = zeros(i,1);A = zeros(i-1,M);Tn=100*eye(M,M);kn =zeros(M,1);for j=M:i-1en(j,1)=xn(j+1)-A(j-1,:)*xn(j:-1:j-M+1)';kn=Tn*xn(j:-1:j-M+1)'/(u+xn(j:-1:j-M+1)*Tn*xn(j:-1:j-M+1)'); Tn=(Tn-Tn*xn(j:-1:j-M+1)'*xn(j:-1:j-M+1)*Tn/(u+xn(j:-1:j-M+1)*Tn*xn(j:-1:j-M+1)'))/u;A(j,:)=A(j-1,:)+kn(:,1)'.*en(j,1);End3、仿真结果分析1、、收敛曲线u=0.0001=0.982、LMS算法与RLS算法的比较由上图可以看出RLS算法的收敛速度比LMS算法的收敛速度快，但是RLS算法在收敛处的波动比LMS算法大。

现代数字信号处理仿真报告3.16计算机模拟实验。

设x(n) = x1(n) + x2(n)，x1(n)是窄带信号，定义为x1(n) = sin(0.05πn +φ)，φ是在[0, 2π)区间上均匀分布的随机相位。

x2(n)是宽带信号，它由一个零均值、方差为1的白噪声信号e(n)激励一个线性滤波器而产生，其差分方程为x2(n) = e(n) + 2e(n-1) + e(n-2)。

(1)计算x1(n)和x2(n)各自的自相关函数，并画出其函数图形。

据此选择合适的延时，以实现谱线增强。

(2)产生一个x(n)序列。

选择合适的μ值。

让x(n)通过谱线增强器。

画出输出信号y(n)和误差信号e(n)的波形，并分别与x1(n)和x2(n)进行比较。

解：一．仿真思路：首先，根据题目要求生成各个信号，包括窄带信号x1(n)、白噪声e(n)、宽带噪声信号x2(n)以及输入信号x(n)。

然后，根据生成的信号，分别计算x1(n)和x2(n)的自相关函数。

其次，观察前述两个自相关数的图形，找出二者的自相关长度，并选择一个介于两自相关长度值之间的一个值，作为Δ值，将输入信号x(n)延迟Δ，即x(n-Δ)。

最后，按照课本上图3.45所示的自适应预测原理完成自适应滤波，其中，自适应滤波器采用LMS算法，算法内容见课本3.7节。

二．数学原理1．自适应的最小均方(LMS)算法LMS算法的核心思想是用平方误差代替均方误差，并用这个平方误差来近似二次性能曲面的梯度。

这样可以推导出LMS算法的权系数基本关系式w(n+1) = w(n) - μ▽(n) = w(n) + 2μe(n)x(n) 滤波器每输入一个数据x(n)，按照当前权系数w(n)滤波输出数据y(n)后，再根据期望信号d(n)和输出信号y(n)计算出相应的误差e(n)，接着按照权系数更新关系式更新权矢量w(n+1)。

下一时刻，再利用新的权矢量w(n+1)对数据进行滤波。

如果μ值选取的合适，那么经过一定的更新步数后，系统的权矢量w将会收敛到一个固定的w*附近，也就是，这个时候，自适应滤波器已达到了稳定滤波阶段。

现代数字信号处理实验报告题目：基于LMS和RLS的自适应滤波器应用仿真基于LMS 和RLS 的自适应滤波器应用仿真一、自适应滤波器自适应滤波器是指利用前一时刻的结果，自动调节当前时刻的滤波器参数，以适应信号和噪声未知或随机变化的特性，得到有效的输出，主要由参数可调的数字滤波器和自适应算法两部分组成，如图1.1所示图1.1 自适应滤波器原理图x(n)称为输入信号，y(n)称为输出信号，d(n)称为期望信号或者训练信号，e(n)为误差僖号，其中，e(n)=d(n)-y(n)，自适应滤波器的系数(权值)根据误差信号e(n)，通过一定的自适应算法不断的进行更新，以达到使滤波器实际输出y(n)与期望响应d(n)之间的均方误差最小。

1.1 自适应滤波器原理图中参数可调的数字滤波器和自适应算法组成自适应滤波器。

自适应滤波算法是滤波器系数权值更新的控制算法，根据输入信号与期望信号以及它们之间的误差信号，自适应滤波算法依据算法准则对滤波器的系数权值进行更新，使其能够使滤波器的输出趋向于期望信号。

记数字滤波器脉冲响应为：h(k)=[h 0(k) h 1(k) … h n-1(k)]T输入采样信号为：x(k)=[x(k) x(k-1) … x(k -n-1)]误差信号为：)()()(^k y k y k e -=()()()()T e k y k h k x k =-优化过程就是最小化性能指标J(k),它是误差的平方和：21()[()()()]kT i J k y i h k x i ==-∑求使J(k)最小的系数向量h(k),即使J(k)对h(k)的导数为零，也就是0)()(=k dh k dJ 。

把J(k)的表达式代入，得：12[()()()]()0kT i y i h k x i x i =-=∑和11()()()()()kkTTT i i xi y i h k x i x i ===∑∑由此得出滤波器系数的最优向量：11()()()()()kTT i k Ti xi y i h k x i xi ===∑∑这个表达式由输入信号自相关矩阵()xx c x 和输入信号与参考信号的相关矩阵()yx c k 组成，如下所示，维数都为（n,n ）:1()()()kT xx i c k x i x i ==∑1()()()kT yx i c k x i y i ==∑系数最优向量也可以写成如下形式：1()()()T opt yx xx h k c k c k -=自相关和互相关矩阵的递归表达式如下：()(1)()()T xx xx c k c k x k x k =-+()(1)()()T yx yx c k c k y k x k =-+把()yx c k 的递归表达式代入系数向量表达式，得：1()()()T yx xx h k c k c k -=即1()[(1)()()]()T T yx xx h k c k x k y k c k -=-+考虑到(1)(1)(1)T yx xx c k h k c k -=--可以记1()()[(1)(1)()()]xx xx h k c x c k h k y k x k -=--+用前面得到的表达式求出(1)xx c k -，并代入上式：1()(){[()()()](1)()()}T xx xx h k c x c k x k x k h k y k x k -=--+或 1()(1)()[()()()()(1)]T xx h k h k c x y k x k x k x k h k -=-+-- 则滤波器系数的递归关系式可以记作1()(1)()[()()()()(1)]T xx h k h k c x y k x k x k x k h k -=-+--其中()()()(1)T e k y k x k h k =--e(k)表示先验误差。

MUSIC算法仿真例题3.6代码Main.mclear all;clc;N=256;M=4;K=3;SNR1=30;SNR2=30;SNR3=27;%设定信噪比；A1=getAk(SNR1);A2=getAk(SNR2);A3=getAk(SNR3);%求得信号的幅度；f1=0.1;f2=0.25;f3=0.27;%设定频率值s1=zeros(1,N);s2=zeros(1,N);s3=zeros(1,N);%初始化sks1=getsk(A1,f1,N);%求sks2=getsk(A2,f2,N);s3=getsk(A3,f3,N);v=randn(1,N);% 构建高斯白噪声；x=s1+s2+s3+v;%构建了时域上的信号N3=2048;%设定画图时描点的数目。

R=getR(x,N,M);[G,D]=getG(R,M,K);d=1/(N3-1);%求画图用的横坐标的间隔。

h=zeros(1,N3);for i=1:N3h(i)=-0.5+(i-1)*d;endy=zeros(1,N3);for j=1:N3w=h(j)*2*pi;aw=getaw(w,M);MinValue=min(y(j));y(j)=10*log10(abs(1/((aw')*G*(G')*aw)));endplot(h,y);getAk.mfunction AK=getAk(SNR) %定义信号幅度AK=((10^(SNR/10))*2)^0.5;Getaw.m %求awfunction aw=getaw(w,M)aw=zeros(M,1);for j=1:Maw(j)=exp(-w*(j-1)*i);endgetR.m %求自相关矩阵Rfunction R=getR(x,N,M)L=N-M+1;tempx=zeros(M,1);R=zeros(M,M);for n=M:Nfor j=1:Mtempx(j)=x(n-(j-1));endR=R+tempx*tempx';endR=R/L;getG.m %求G矩阵function [G,D]=getG(R,M,K)[V,D]=eig(R);G=zeros(M,M-K);z=zeros(M,1);for j=1:Mz(j)=D(j,j);%将特征值放入了z里面end[z,y]=sort(z);%对z进行了排序目的是，找到最小的M-K个特征值多对应的特征向量。