现代数字信号处理仿真作业
现代数字信号处理及其应用 第六章仿真题
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算法大。
信号处理仿真题作业
信号处理基础仿真作业学号：S*********姓名：***3.17在计算机上用如下方法产生随机信号()u n的观测样本：首先产生一段零均值、方差为2σ的复高斯白噪声序列()v n；然后在()v n上叠加三个复正弦信号，它们的归一化频率分别是f1=0.15，f2=0.17和f3=0.26。
调整2σ和正弦信号的幅度，使在f1、f2和f3处得信噪比分别为30dB、30dB和27dB。
（1）令信号观测样本长度N=32，试用3.1.1节讨论的基于FFT的自相关函数快速计算方法估计出自相关函数^()0mr，并与教材式（3.1.2）估计出的自相关函数^()mr做比较。
产生零均值、方差为1的复高斯白噪声序列y >> y=randn(1,32);>> y=y-mean(y);>> y=y/std(y);>> a=0;>> b=sqrt(2);>> y=a+b*y产生三个复正弦信号并产生观察样本：>> N=32;>> 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+y基于FFT的自相关函数快速计算方法：N=32;>> Uk=fft(un, 2*N);Sk=(1/N)* abs(Uk).^2;r0=ifft(Sk);r1=[r0(N+2:2*N),r0(1:N)];>> figure(1);>> stem(real(r1));>> figure(2);>> stem(imag(r1))输出结果为：图 1 基于FFT的自相关函数快速计算实部：虚部：教材中式（3.1.2）估计自相关函数>> r=xcorr(un, N-1,'biased');>> figure(1);>> stem(real(r))>> figure(2);>> stem(imag(r))输出结果为：图 2 教材式（3.1.2）估计的自相关函数实部：虚部：（2）令信号观测样本长度N=256，试用BT法和周期图法估计()u n的功率谱，这里设BT法中所用自相关函数的单边长度M=64。
数字信号处理仿真作业
作业1（1）求2σvA = [1 -a1 -a2; 0 -(1+a2) 0; 0 -a1 -1];b = [r0 a1*r0 a2*r0]';x = inv(A)*b; %x = [sigma_v r1 r2]';sigma_v=x(1); r1=x(2); r2=x(3);只要更改a1和a2的值，得到a1 a2 2vσP1P2r(1) r(2) -0.195 0.95 0.0965 0.0975+0.9698i 0.0975-0.9698i 0.1000 -0.9305 -1.595 0.95 0.0323 0.7975+0.5604i 0.7975-0.5604i 0.8179 0.3546 -1.9114 0.95 0.0038 0.9557+0.1914i 0.9557-0.1914i 0.9802 0.9236 （2）产生)u序列(ndata_len = 1024;trials = 100;n = 1:data_len;a1 = -1.9114;a2 = 0.95;sigma_v_2 =0.003822;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)un=zeros(1,1024)for h=1:1024un(h)=u(h);end（3）)(k r 的估计 p = 128;r0 = xcorr(u,p,'unbiased')'; r = r0(p+1:2 * p+1)' stem(1:257,r0)1002003004005006007008009001000-1-0.8-0.6-0.4-0.200.20.40.60.81r(m)r(0)=0.99579（4）S AR (ω)、S B T (ω)、S PER (ω) S AR (ω)图和程序如下： NF=1024;Spr=fftshift((1/NF)*abs(fft(un,NF)).^2)' Spr=Spr/max(Spr); Spr_unit=10*log10(Spr); plot(1:1024,Spr_unit);S BT(ω)图和程序如下：M=128;r=xcorr(un,M,'biased'); NF=1024;BT=fftshift(fft(r,NF)); BT=BT/max(BT);BT_unit=10*log10(BT); plot(1:1024,BT_unit);S PER(ω)图和程序如下：NF=1024;Par=sigma_v_2./fftshift(abs(fft([1,a1,a2],NF)).^2); Par=Par/max(Par); Par_unit=10*log10(Par); plot(1:1024,Par_unit);作业2（1）产生)(n u 序列 close all ; clear; clc;N=1024; Pv=-10;%-10dB snr=10;a1=10^(snr/20); n=1:N;vn=randn(1,N);un=exp(j*pi*(1.4.*n-1))+exp(j*pi*(1.8.*n-0.9))+a1*vn; figure(); plot(n,un);xlabel('n');ylabel('un')title('产生序列un')（2）)r的估计(kfor m=0:500rm_e(m+1)=un(m+1:N)*un(1:N-m)'/(N-m); %r(m)估计endfigure();plot(1:500,rm_e(1:500));xlabel('m');ylabel('r(m)')title('自相关图');（3）估计四阶样本自相关矩阵M=4;for k=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=xs*xs'/(N-M);（4）S BT(ω)L=256;r=xcorr(un,L,'biased');NF=2048;BT=fftshift(fft(r,NF));BT=BT/max(BT);BT_unit=10*log10(BT);plot(1:2048,BT_unit);S MVDR(ω)NF=2048;for n=1:NFAq=exp(-j*2*pi*(n-1)/NF*(0:M-1)');Pmvdr(n)=1/(Aq'*inv(R)*Aq);endPmvdr=Pmvdr/max(Pmvdr);Pmvdr_unit=10*log10(Pmvdr);plot(1:2048,Pmvdr_unit);S AR(ω)首先估计出自相关函数值r(0),r(1)...r(16);利用Levinson-Durbin迭代算法实现AR模型的系数的求出；最后计算16阶AR模型的功率谱。
数字信号处理实验作业
实验5 抽样定理一、实验目的：1、了解用MA TLAB 语言进行时域、频域抽样及信号重建的方法。
2、进一步加深对时域、频域抽样定理的基本原理的理解。
3、观察信号抽样与恢复的图形，掌握采样频率的确定方法和内插公式的编程方法。
二、实验原理：1、时域抽样与信号的重建 （1）对连续信号进行采样例5-1 已知一个连续时间信号sin sin(),1Hz 3ππ=0001f(t)=(2f t)+6f t f ，取最高有限带宽频率f m =5f 0，分别显示原连续时间信号波形和F s >2f m 、F s =2f m 、F s <2f m 三情况下抽样信号的波形。
程序清单如下：%分别取Fs=fm ，Fs=2fm ，Fs=3fm 来研究问题 dt=0.1; f0=1; T0=1/f0; m=5*f0; Tm=1/fm; t=-2:dt:2;f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); subplot(4,1,1); plot(t,f);axis([min(t),max(t),1.1*min(f),1.1*max(f)]); title('原连续信号和抽样信号'); for i=1:3;fs=i*fm;Ts=1/fs; n=-2:Ts:2;f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); subplot(4,1,i+1);stem(n,f,'filled');axis([min(n),max(n),1.1*min(f),1.1*max(f)]); end程序运行结果如图5-1所示：原连续信号和抽样信号图5-1（2）连续信号和抽样信号的频谱由理论分析可知，信号的频谱图可以很直观地反映出抽样信号能否恢复原模拟信号。
因此，我们对上述三种情况下的时域信号求幅度谱，来进一步分析和验证时域抽样定理。
例5-2编程求解例5-1中连续信号及其三种抽样频率（F s>2f m、F s=2f m、F s<2f m）下的抽样信号的幅度谱。
现代信号处理作业
if (fpe(m)-fpe(m+1))<1e-6 %用FPE准则定阶
break
end
m=m+1;
end
end
disp(['根据模型求出的自相关序列为：',num2str(R)])
disp(['题目给出的AR模型参数为：',num2str(a1)])
R1 = [R(3),R(2),R(1);R(4),R(3),R(2);R(5),R(4),R(3)];
r = [-R(4);-R(5);-R(6)];
A=linsolve(R1,r);%估计模型参数。
A
% 2)估计ARMA模型中的MA参数。
AA=[1,A(1),A(2),A(3)];
v =filter(AA,1,x);
结果与分析：
根据模型求出的自相关序列为：4.9377 4.3287 4.1964 3.8654 3.6481 3.4027 3.1919 2.986 2.797
题目给出的AR模型参数为：-0.58333 -0.375 0.041667
估计出的AR模型参数为：-0.56844 -0.36884 0.039549
（3）由于已经估计出ARMA模型的所有参数，根据ARMA模型与功率谱等价的关系可以算出功率谱，并用plot函数画出图形。
在这里，噪声方差为一，各阶参数都已经估计出，进而画出图形。
结果与分析：
A =
-0.2630
-0.1870
0.0537
B =
0.4696
0.0317
由结果可知，求解结果与已知系数相近，但存在一定误差，当观测的数据个数不断增大时，结果误差会不断减小，一直到不再用明显变化。
数字信号处理大作业(matlab仿真)
《动态数据处理》综合作业2019.1.20目录1 差分方程求解 (3)1.1 基本原理 (3)1.2 测试程序 (4)1.3 运行结果 (4)2 线性卷积 (5)2.1 基本原理 (5)2.2 测试程序 (7)2.3 运行结果 (8)3 DFT与IDFT程序设计 (9)3.1 基本原理 (9)3.2 测试程序 (11)3.3 运行结果 (13)4 基于Hilbert变换的信号包络提取 (14)4.1 基本原理 (14)4.2 测试程序 (15)4.3 运行结果 (17)5 基于窗函数法的FIR滤波器设计 (17)5.1 基本原理 (17)5.2 测试程序 (19)5.3 运行结果 (21)6 相关函数估计 (23)6.1 基本原理 (23)6.2 测试程序 (24)6.3 运行结果 (25)7 基于分段平均的功率谱估计 (25)7.1 基本原理 (25)7.2 测试程序 (27)7.3 运行结果 (28)8 源程序 (29)8.1 差分方程 (29)8.2 线性卷积 (29)8.3 离散傅里叶变换 (30)8.4 基于窗函数的FIR滤波器 (31)8.5 自相关函数 (34)8.6 welch法功率谱估计 (34)1 差分方程求解1.1 基本原理一个线性移不变离散系统可以用一个常系数线性差分方程来描述：其中，、为方程的系数（k=1，…，N；r=0，…，M）。
若给定输入信号及系统的初始条件，则可求出系统的输出，即方程的解。
利用Z变换求解差分方程，在初始条件为零的情况下：1.2 测试程序已知离散系统的差分方程为，输入信号为，求系统响应。
测试程序如下：clc;clear alla=[1,-3/4,1/8];b=[1,1/3,0]; %输入差分方程系数向量，不足补零对齐n=0:0.1:15;xn=(3/4).^n; %输入激励信号yn=filterRe(b,a,xn); %调用filterRe函数plot(n,xn,'-',n,yn,'k-');xlabel('n');ylabel('x(n),y(n)');title('差分方程');1.3 运行结果运行后得到如下信号图，其中下方为输入信号，上方为输出信号，即系统响应。
现代数字信号处理仿真作业
现代数字信号处理仿真作业1.仿真题仿真结果及图形：图 1 基于FFT的自相关函数计算图 2图 3 周期图法和BT法估计信号的功率谱图 4 利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为：a1=--;a2=;a3=3i;a4=7;a5=68i;a6=76i;a7=92i;a8=20ia9=20i;a10=23i;a11=710i;a13=83i ;a14=24i;a15=21i;a16=3i.仿真程序（3_17）:clear allclc%% 产生噪声序列N=32; %基于FFT的样本长度%N=256; %周期图法，BT法，AR模型功率谱估计的长度vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);%%产生复正弦信号f=[ ]; %归一化频率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;%% 利用的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: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.仿真题仿真结果及图形：单次Root-MUSIC算法中最接近单位圆的两个根为：对应的归一化频率为：相同信号的MUSIC谱估计结果如下图 5 对信号进行MUSIC谱估计的结果仿真程序（3_20）:clear allclc%% 信号样本和高斯白噪声的产生N=1000;vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal=[exp(1i**pi*(0:N-1)+1i*2*pi*rand); %复正弦信号exp(-1i**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*+(n-1)/(NF-1))*(0:M-1)');Pmusic(n)=1/(Aq'*En*En'*Aq);endkk=+(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.仿真题仿真结果及图形：单次ESPRIT算法中最接近单位元的两个特征值为：对应的归一化频率为：仿真程序（3_21）:clear allclc%% 信号样本和高斯白噪声的产生N=1000;vn=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal=[exp(1i**pi*(0:N-1)+1i*2*pi*rand); %复正弦信号exp(-1i**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.仿真题仿真结果及图形：步长为时失调参数为m1=；步长为时失调参数为m2=。
现代数字信号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 的点数。
现代数字信号处理理论计算法仿真作业
现代数字信号处理仿真作业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 r11N=32;u = u_length(N); %长度为32的随机信号for m = 1:Nr11(m) = u(m:N)*u(1:N-m+1)'/N;endr12 = 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 onk = -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 onBT法估计功率谱密度仿真结果如下图：%周期图法估计功率谱密度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:pl=conj(r(m:-1:2)');k=-(r(m+1)+a*l)/c;for n=1:m-1b(n)=a(n)+k*conj(a(m-n));enda=[b k];c=c*(1-k*k');endh=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)));gridtitle('阶数为16的AR模型功率谱密度估计曲线')仿真结果如下图：3.20clear;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:Mfor m=1:MR(m,n)=r(M+1+n-m);endend[G,D]=eigs(R,M-K,'sm');for n=-100:100;for m=1:Ma(m)=exp(j*(m-1)*pi*n/100);endp(n+101)=1/(a*(G*G')*a');endp=abs(p);n=-1:0.01:1;figure;plot( n,10*log10(p/max(p)) );title(¡®MUSICËã·¨½á¹û¡¯) Gr=G*G';a=zeros(2*M-1,1);for m=1:Ma(m:m+M-1)=a(m:m+M-1)+Gr(M:-1:1,m);endx=roots(a);err=abs(abs(x)-1);w=angle(x)/pi;e=sort(err);for n=1:2*Kt(n)=find(err==e(n));endfigure;bar(w(t));gridroot—music 估计仿真结果：MUSIC算法频率估计结果：3.21.ESPRIT算法频率估计，仿真程序如下：N=1000;M=8;v=sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N);n=1:N;a=2*pi*rand(1,2);u(n)=exp(j*0.5*pi*n+j*a(1))+exp(j*1.2*pi*n+j*a(2))+v(n); n=1:N-1;d(n)=u(n+1);%ESPRIT估计r=xcorr(u,M,'biased');for n=1:Mfor m=1:MR(m,n)=r(M+1+n-m);Rd(m,n)=r(M+n-m);endendla=eig(R);k=min(la);Z=zeros(M,M);for n=1:M-1;Z(n,n+1)=1;endC=R-k*eye(M);Cd=Rd-k*Z;ew=eig(C,Cd); err=abs(abs(ew)-1);w=angle(ew)/pi; e=sort(err);n1=find(err==e(1));n2=find(err==e(2));ws=[w(n1),w(n2)];stem(ws);4.18.仿真程序源代码与仿真结果如下：%lms预测算法源程序clear%channel system ordersysorder = 2 ;a=[1 -0.975 0.95];mu=0.05; %预测步长% Number of system pointsN=512;for m=1:100;v = randn(1,N)*0.0731; % Input to the filterx = filter(1,a,v); % 产生的需要信号w = zeros ( sysorder , N-sysorder ) ;for n = sysorder+1 : Nu = x(n-1:-1:n-sysorder)' ;y(n)= w(:,n-sysorder)' * u;e(m,n) = x(n) - y(n) ;w(:,n-sysorder+1)= w(:,n-sysorder)+ mu * u * e(m,n) ; endende=sum(e.^2)/100;plot(e);title('预测误差的变化曲线');figurehold onplot (w(1,:),'r');plot (w(2,:),'b');title('权向量的估计')hold off下面的两幅图为步长μ分别取0.05和0.005时，权向量的变化结果：可以看出，步长的不同影响了收敛的速度，但步长变大时使得失调参数增大。
数字信号处理仿真实验
数字信号处理上机实验仿真报告学院电子工程学院班级学号姓名实验一： 信号、 系统及系统响应一、实验目的(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 ˆj e d Ωt a a n X Ωx t t nT t δ∞∞--∞=-∞⎡⎤=-⎢⎥⎣⎦∑⎰ ()()j e d Ωt a n x t t nT t δ∞∞--∞=-∞=-∑⎰()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) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。
现代信号处理大作业
姓名：潘晓丹学号：班级：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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 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。