现代数字信号处理仿真作业
现代数字信号处理及其应用 第六章仿真题
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,:));。
现代数字信号处理仿真作业
现代数字信号处理仿真作业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均值的高斯白噪声。
现代数字信号处理作业2014
现代数字信号处理作业1、请用MATLAB编程举例阐述一维信号压缩与重建（P349 例8.7.1)。
2、请用MATLAB编程举例阐述图像压缩与重建（P350 例8.7.2)。
3、请举例用MATLAB的Spectrum.m进行功率谱估计(P530 13.5)。
4、请用MATLAB编程举例介绍现代普估计各种算法的优缺点（至少3种，P584）。
5、请用MATLAB编程举例阐述自适应谱线增强（P631)。
6、请用MATLAB编程举例阐述自适应滤波在系统辨识中的运用(P625 16.5)。
7、请用MATLAB编程举例阐述自适应噪声抵消（P628 16.5.2)。
8、请用MATLAB编程举例阐述自适应预测（P630 16.5.3)。
9、请用MATLAB编程举例阐述自适应均衡（P632 6.5.4 自适应均衡）。
10、请用MATLAB编程举例阐述卡尔曼滤波的运用。
11有限冲激响应（FIR）数字滤波器设计的设计，请用MATLAB编程举例阐述12图像DCT变换及图像压缩，，请用MATLAB编程举例阐述13自适应抵消的应用，请用MATLAB编程举例阐述14经典功率谱分析应用，请用MATLAB编程举例阐述15卡尔曼滤波应用的一个例子，请用MATLAB编程举例阐述16自适应滤波器应用的一个例子，请用MATLAB编程举例阐述17 谢锦华-基于DCT的一维信号的压缩与重建的一个例子，请用MATLAB编程举例阐述18 自适应均衡应用的一个例子，请用MATLAB编程举例阐述19自适应预测应用的一个例子，请用MATLAB编程举例阐述20 自适应滤波器应用的一个例子，请用MATLAB编程举例阐述21 图像加密与解密的一个例子，请用MATLAB编程举例阐述。
现代信号处理作业
现代信号处理作业1.总结学过的滤波器设计方法，用matlab仿真例子分析不同设计方法的滤波器的性能及适应场合。
答复:1.1模拟低通滤波器的设计方法1.1.1butterworth滤波器设计步骤：⑴.确定阶次n① 已知ωc、ωS和as，求出butterworthdf的阶数n1由：a??10lgh(j?)??10lg2nsas1?(?s/?c)a/10lg(10s?1)求出n：n?2lg(?s/?c)②已知ωc、ωs和ω=ωp(3db)的衰减ap求butterworthdf阶数nPlg(10p?1)n?2lg(?/?)pca/102得到n：③已知ωp、ωs和ω=ωp的衰减ap和as求butterworthdf阶数n1a？？10lgh（j？）？？10lg由Pap1提供？（？P/？C）2n然后：（？P/？C）2n？10ap/102？1，（？s/？c）2n？10as/10？一求出n:lg[(10ap/10?1)(10as/10?1)]n?2lg(?p/?s)⑵.用阶次n确定ha(s)根据公式：|哈（j？）|2s=ha（s）ha（？s）？=J1，分母？0，2n1？（s/j？c）12k？1j[？]22nsk？（？1）12n（j？c）？？ce，k？左半平面上1,2,2nha（s）ha（？s）的极点就是ha（s）的极点，所以hs(s)?n?c?(s?s)kk?1nsk??ce12k?1j[?]?22n，k?1,2,,n1.1.2切比雪夫低通滤波器设计步骤：⑴.确定技术指标?p?p?s?s正常化：？Pp/？P1.ss/？P⑵. 计算过滤器阶数n和？：ch?1(k1?1)100.1??1?1n?其中k1??10.1?ch?s10?1sp0.1?2?1p??10⑶.求出归一化系统函数其中极点由下式求出：圆周率？？嘘？sin[ha（p）-（2k？1）-（2k？1）]？jch？cos[]2n1？？2n？1.（p？p）ii？1便士？s/？pnha（s）？Ha（P）或Ha（P）可以通过直接从N和s中查找表来获得?s?⑷.去归一化：ha(s)?h（ap）=hap?2.数字低通滤波器的设计步骤：（1）确定数字低通滤波器的技术指标：通带截止频率、阻带截止频率、阻带最小衰减系数？s（2）将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。
数字信号处理课程仿真大作业
数字信号处理课程仿真大作业一、课程设计题目：基于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 进行频谱谐波理论理论分析：在信号处理中，信号的频谱分析是重要的应用领域之一。
现代数字信号处理ADSP仿真报告
现代数字信号处理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算法大。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
现代数字信号处理仿真作业1.仿真题3.17仿真结果及图形：图 1 基于FFT的自相关函数计算图 3 周期图法和BT 法估计信号的功率谱图 2 基于式3.1.2的自相关函数的计算图 4 利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为：a1=-0.402637623107952-0.919787323662670i;a2=-0.013530139693503+0.024214641171318i;a3=-0.074241889634714-0.088834852915013i;a4=0.027881022353997-0.040734794506749i;a5=0.042128517350786+0.068932699075038i;a6=-0.0042799971761507 + 0.028686095385146i;a7=-0.048427890183189 - 0.019713457742372i;a8=0.0028768633718672 - 0.047990801912420ia9=0.023971346213842+ 0.046436389191530i;a10=0.026025963987732 + 0.046882756497113i;a11= -0.033929397784767 - 0.0053437929619510i;a12=0.0082735406293574 - 0.016133618316269i;a13=0.031893903622978 - 0.013709547028453i ;a14=0.0099274520678052 + 0.022233240051564i;a15=-0.0064643069578642 + 0.014130696335881i;a16=-0.061704614407581- 0.077423818476583i.仿真程序（3_17）:clear allclc%% 产生噪声序列N=32; %基于FFT的样本长度%N=256; %周期图法，BT法，AR模型功率谱估计的长度vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);%%产生复正弦信号f=[0.15 0.17 0.26]; %归一化频率SNR=[30 30 27]; %信噪比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)];%% 利用3.1.2估计Rr2=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算法中最接近单位圆的两个根为：-0.001156047541561 + 1.001503153449793i0.587376604261220 - 0.810845628739986i对应的归一化频率为：0.250183714447964-0.150223406926494相同信号的MUSIC谱估计结果如下图 5 对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')仿真结果及图形：单次ESPRIT算法中最接近单位元的两个特征值为：0.001826505974929 + 1.000690248438859i0.586994191014025 - 0.809491260856630i对应的归一化频率为：0.249709503383161-0.150146235268272仿真程序（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);仿真结果及图形：步长为0.05时失调参数为m1=0.0493；步长为0.005时失调参数为m2=0.0047。
图 6 步长为0.05时权向量的收敛曲线图 7 步长为0.005时权向量的收敛曲线图 8 步长分别为0.05和0.005时100次独立实验的学习曲线仿真程序（4_18）:clear allclc%% 产生100组独立样本序列data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;v=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,v,Zi); %产生100组独立信号%% LMS迭代mu1=0.05;mu2=0.005;w1=zeros(2,data_len,trials); w2=w1;for m=1:100;temp=zeros(data_len,1);e1(:,:,m)=temp;e2(:,:,m)=temp;d1(:,:,m)=temp;d2(:,:,m)=temp;for n=3:data_len-1w1(:,n+1,m)=w1(:,n,m)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n,1,m)); w2(:,n+1,m)=w2(:,n,m)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n,1,m)); d1(n+1,1,m)=w1(:,n+1,m)'*u(n:-1:n-1,:,m);d2(n+1,1,m)=w2(:,n+1,m)'*u(n:-1:n-1,:,m);e1(n+1,1,m)=u(n+1,:,m)-d1(n+1,1,m);e2(n+1,1,m)=u(n+1,:,m)-d2(n+1,1,m);endendt=1:data_len;w1_mean=zeros(2,data_len);w2_mean=w1_mean;e1_mean=zeros(data_len,1);e2_mean=e1_mean;for m=1:100w1_mean=w1_mean+w1(:,:,m);w2_mean=w2_mean+w2(:,:,m);e1_mean=e1_mean+e1(:,:,m).^2;e2_mean=e2_mean+e2(:,:,m).^2;endw1_mean=w1_mean/100; %100次独立实验权向量的均值w2_mean=w2_mean/100;e1_100trials_ave=e1_mean/100; %100次独立实验的学习曲线均值e2_100trials_ave=e2_mean/100;figure(1)plot(t,w1(1,:,90),t,w1(2,:,90),t,w1_mean(1,:),t,w1_mean(2,:))xlabel('迭代次数');ylabel('权向量')title('步长=0.05')figure(2)plot(t,w2(1,:,90),t,w2(2,:,90),t,w2_mean(1,:),t,w2_mean(2,:)) xlabel('迭代次数');ylabel('权向量')title('步长=0.005')%% 计算剩余误差和失调参数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);endsJmin=sum(Jmin)/trials;Jex1=e1_100trials_ave-sJmin; %剩余均方误差mu1Jex2=e2_100trials_ave-sJmin; %剩余均方误差mu2sum_eig_100trials=sum(sum_eig)/100;Jexfin1=mu1*sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials)); Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials)); M1=Jexfin1/sJmin; %失调参数m1M2=Jexfin2/sJmin; %失调参数m2figure(3)plot(t,e1_100trials_ave,'*',t,e2_100trials_ave)xlabel('迭代次数');ylabel('均方误差')legend('u1=0.05','u2=0.005')axis([0,600,0,1])5.仿真题5.10仿真结果及图形：(1) M=2时，210.99,0.93627a σ=-= ，求解Yule-Walker 方程211(0)(1)(1)(0)0r r a r r σ⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可得到自相关矩阵247.048746.578346.578347.0487R ⎡⎤=⎢⎥⎣⎦相应的计算程序为r2=inv([1,-0.99;-0.99,1])*[0.93627;0];R2=[r2(1),r2(2);r2(2),r2(1)]; % M=2(2) M=3时，2120.99,0,0.93627a a σ=-== ，求解Yule-Walker 方程212(0)(1)(2)1(1)(0)(1)0(2)(1)(0)0r r r r r r a r r r a σ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可得到自相关矩阵为347.048746.578346.112546.578347.048746.578346.112546.578347.0487R ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦相应的计算程序为r3=inv([1,-0.99,0;-0.99,1,0;0,-0.99,1])*[0.93627;0;0];R3=[r3(1),r3(2),r3(3);r3(2),r3(1),r3(2);r3(3),r3(2),r3(1)]; % M=3(3) 计算特征值扩展%% 特征值分解eig_value_1=eig(R2);eig_value_2=eig(R3);%% 特征值扩展eig_spread_1=max(eig_value_1)/min(eig_value_1);eig_spread_2=max(eig_value_2)/min(eig_value_2);M=2时特征值扩展是199.0000；M=3时特征值扩展是444.2790。