《现代数字信号matlab处理仿真题》由会员分享，可在线阅读，更多相关《现代数字信号matlab处理仿真题（28页珍藏版）》请在金锄头文库上搜索。
1、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、 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);yla
3、bel(实部);xlabel(估计的自相关函数的实部);subplot(2,2,4);stem(imag_r1);xlabel(估计的自相关函数的虚实部);ylabel(虚部);function AK=getAk(SNR) %求得幅度%由SNR=10log(A2/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可以补零，确定
4、实际求fft的点数。N=256; %样本观察长度M=64; %加窗长度f1=0.15;f2=0.17;f3=0.26; %定义3个归一化正弦频率SNR1=30;SNR2=30;SNR3=27; %定义三个正弦信号信噪比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个复正弦信号un=s1+s2+s3+noise;%下面采用周期法估
5、计的频谱%N2=2*N;u2=zeros(1,N2);u2(1:N)=un;%对信号序列进行补零Uw=fft(u2);%求DFT变换Sw1=zeros(1,N2);for w=1:N2 Sw1(w)=(abs(Uw(w)2)/N;%计算功率谱endr0=zeros(1,N2);r0=ifft(Sw1);r=zeros(1,N2);r(1:N)=r0(N+1:N2);r(N+1:N2)=r0(1:N);M=64;%加窗处理for n=1:2*Mr1(n)=r(N2-2*M)/2+n);end%加窗之后的序列为r1N3=256;%实际进行fft变换的点数。r2=zeros(1,N3);r2(1:2
6、*M)=r1;%将加窗之后的r1进行了补零得到的是r2.Sw2=fft(r2);Sw2=log10(abs(Sw2);%取模取对数temp=zeros(1,N3/2);temp=Sw2(1:N3/2);Sw2(1:N3/2)=Sw2(N3/2+1):N3);Sw2(N3/2+1):N3)=temp;%将0-2*pi的序列折换成-pi到pi的序列。d=1/(N3-1);x=zeros(1,N3);for i=1:N3 x(i)=-0.5+(i-1)*d;endSw1=log10(abs(Sw1);%取模取对数temp=zeros(1,N2/2);temp=Sw1(1:N2/2);Sw1(1:N2
7、/2)=Sw1(N2/2+1):N2);Sw1(N2/2+1):N2)=temp;%将0-2*pi的序列折换成-pi到pi的序列。l=1/(N2-1);y=zeros(1,N2);for k=1:N2 y(k)=-0.5+(k-1)*l;endsubplot(1,2,1);plot(x,Sw2);xlabel(BT法);ylabel(归一化的功率谱/dB);subplot(1,2,2);plot(y,Sw1);xlabel(周期法);ylabel(归一化的功率谱/dB);仿真波形：由仿真结果可以看出：BT法和周期法在相同的阶数下，周期法比BT法分辩率更高，但是波动比BT法更大。（3）Levin
8、son-Durbin迭代法仿真代码：clear all;clc;%设定N值可以改变抽样信号的点数，设定M值可以设定加窗的大小，设定N3可以补零，确定实际求fft的点数。N=256; %样本观察长度p=16; %设定AR模型参数f1=0.15;f2=0.17;f3=0.26; %定义3个归一化正弦频率SNR1=30;SNR2=30;SNR3=27; %定义三个正弦信号信噪比A1=getAk(SNR1);A2=getAk(SNR2);A3=getAk(SNR3); %求得信号的幅度；noise=randn(1,N) + j *randn(1,N); % 构建高斯白噪声；s1=getSk(A1,f1
9、,N);s2=getSk(A2,f2,N);s3=getSk(A3,f3,N); %产生3个复正弦信号un=s1+s2+s3+noise;r0=xcorr(un, p,biased); %直接计算自相关函数rp=r0(p+1:2*p+1);%Lervinson_Durbin算法 %a(1,1)=-rp(2)/rp(1);e(1)=rp(1)-abs(rp(2)2/rp(1); %计算1b=0;for m=2:p k(m)=-(rp(m+1)+sum(a(m-1,1:m-1).*rp(m:-1:2)/e(m-1); a(m,m)=k(m); for d=1:m-1 a(m,d)=a(m-1,d)
10、+k(m)*conj(a(m-1,m-d); end e(m)=e(m-1)*(1-abs(k(m)2);end%计算16阶AR功率谱%NF=1024; %FFT点数Sw=e(p)./(abs(fft(1,a(p, :),NF).2);%此下作用相当于函数fftshifttemp=zeros(1,NF/2);temp=Sw(1:NF/2);Sw(1:NF/2)=Sw(NF/2+1):NF);Sw(NF/2+1):NF)=temp;%将0-2*pi的序列折换成-pi到pi的序列。SAR=zeros(1,NF);for w=1:NF SAR(w)=log10(abs(Sw(w)2);endd=1/
11、(NF-1);x=zeros(1,NF);for i=1:NF x(i)=-0.5+(i-1)*d;endplot(x,SAR);xlabel(w/2);ylabel(归一化的功率谱/dB);仿真波形：3.20（1）Root-MUSIC仿真代码：clear all;clc;N=1000; %样本个数f1=0.5; %信号频率 f2=-0.3;s1=getSk(f1,N); %产生第一个信号s2=getSk(f2,N); %产生第二个信号noise=(randn(1,N) + j*randn(1,N)/sqrt(2); % 构建高斯白噪声；un=s1+s2+noise; %产生带噪声的信号K=2
12、; %信号源个数M=8; %自相关矩阵阶数R=getR(un,N,M); %根据教材3-5-18求自相关矩阵%root_MUSIC进行频率估计%U,E=svd(R); %矩阵U是特征值向量组成的矩阵，E是对角矩阵，对角矩阵的排列是从大到小e=diag(E);G=U(:,3:M);GG=G*G;co=zeros(2*M-1,1); %由教材可知是关于变量z的2（M-1）次方程，共有2（M-1）个根for m=1:M co(m:m+M-1)=co(m:m+M-1)+GG(M:-1:1,m);endz=roots(co); %求多项式的跟erro=abs(abs(z)-1);%除掉大于1的增根（方程
13、式在此条件下变形的）ph=angle(z)/(2*pi);%挑选出最接近单位圆的N个根disp(z)disp(erro)disp(ph)function aw=getaw(w,M)aw=zeros(M,1);for j=1:M aw(j)=exp(-w*(j-1)*i);endfunction R=getR(x,N,M)L=N-M+1;tempx=zeros(M,1);R=zeros(M,M);for n=M:N for j=1:M tempx(j)=x(n-(j-1); end R=R+tempx*tempx;endR=R/L; function Sk=getSk(fk,N) Sk=exp(j * 2 * pi * fk *(0:N-1)+j*2*pi*rand); %产生0到2上均分布的随机相位信号仿真结果：w1= -0.3004 w2= 0.4996（2）MUSIC算法仿真程序：clear all;clc;N=1000; %样本个数f1=0.5; %信号频率 f2=-0.3;s1=getSk(f1,N); %产生第一个信号s2=getSk(f2,N); %产生第二个信号noise=(randn(1,N) + j*randn(1,N)/sqrt(2); % 构建高斯白噪声；un=s1+s2+noise;
