信号处理仿真题作业

信号处理基础仿真作业学号：S120101057姓名：贾雪婷3.17在计算机上用如下方法产生随机信号()un的观测样本：首先产生一段零均值、方差为2的复高斯白噪声序列()vn；然后在()vn上叠加三个复正弦信号，它们的归一化频率分别是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的复高斯白噪声序列yy=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法和周期图法估计()un的功率谱，这里设BT法中所用自相关函数的单边长度M=64。周期图法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));xn=signal1+signal2+signal3+noiseNF=1024;Spr=fftshift((1/NF)*abs(fft(xn,NF)).^2);Spr=Spr/max(Spr);t=[-0.5:1/NF:0.5-1/NF];plot(t,10*log10(Spr));图3周期图法所得功率谱BT法M=64;r=xcorr(xn,M,'biased');NF=1024;BT=fftshift(fft(r,NF));BT=BT/max(BT);plot((-511:512)/1024,10*log10(BT))图4BT法所得功率谱（3）令信号观测样本长度N=256，试用Levinson-Durbin迭代算法求解AR模型的系数并估计()un的功率谱，模型的阶数取为p=16。p=16;r0=xcorr(xn,p,'biased');r=r0(p+1:2*p+1);a(1,1)=-r(2)/r(1);sigma(1)=r(1)-(abs(r(2))^2)/r(1);form=2:pk(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1);a(m,m)=k(m);fori=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);endNF=1024;Par=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);Par=Par/max(Par);plot((-511:512)/1024,10*log10(Par))图516阶AR模型的功率谱估计3.18设随机过程()un为12()exp(0.5)exp(0.3)()unjnjjnjvn其中，()vn是零均值、方差为1的白噪声，1、2是相互独立并在[0,2]上服从均匀分布的随机相位。使用MVDR方法进行信号频率估计的方针实验，画出频率估计谱线，并给出正弦信号频率的估计值。（要求：信号样本数取1000，估计的自相关矩阵为8阶）解：clearall;%产生0均值，方差为1的复高斯白噪声序列v(n)N=1000;FS=1;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);%产生带噪声的信号样本u(n)signal1=exp(j*0.5*pi*(0:N-1)/FS+j*2*pi*rand);signal2=exp(-j*0.3*pi*(0:N-1)/FS+j*2*pi*rand);un=signal1+signal2+noise;%利用un来估计自相关函数r=zeros(1,N);form=1:Nforn=m:Nr(m)=r(m)+un(n)*(un(n-(m-1)))';endr(m)=r(m)/N;endM=8;%设定横向滤波器的阶数%利用自相关矩阵和自相关函数的关系，构建自相关矩阵R=zeros(M,M);fori=1:Mforj=1:Mifi=jR(i,j)=r(1+j-i);%得到了M阶的矩阵RelseR(i,j)=(r(1+i-j))';endendendN3=2048;%设定画图时描点的数目。d=1/(N3-1);%求画图用的横坐标的间隔。h=zeros(1,N3);fori=1:N3h(i)=-0.5+(i-1)*d;endy=zeros(1,N3);forj=1:N3w=h(j)*2*pi;y(j)=10*log10(abs(getPMVDR(w,M,R)));endplot(h,y);title('MVDR谱估计');xlabel('\omega/2\pi');ylabel('归一化功率谱/dB');axis([-0.50.5-161]);图6MVDR谱估计3.19复正弦加白噪声随机过程()un同题3.18中所给。试用MUSIC方法进行信号频率估计的仿真实验。要求：信号样本数取1000，估计的自相关矩阵为8阶）（1）采用AIC和MDL准则估计信号源个数；N=1000;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);signal1=exp(j*0.5*pi*(0:N-1)+j*2*pi*rand);signal2=exp(-j*0.3*pi*(0:N-1)+j*2*pi*rand);un=signal1+signal2+noise;M=8;fork=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=xs*xs'/(N-M);%自相关矩阵的特征值分解[U,E]=svd(R);%U是特征矢量组成的矩阵，E是对角阵，对角元素是由大到小排列的特征值ev=diag(E);%提取对角元素上的特征值；%根据AIC准则进行信号源个数的估计fork=1:Mdec=prod(ev(k:M).^(1/(M-k+1)));%计算第一项中对数的自变量的分子nec=mean(ev(k:M));%计算第一项中对数的自变量的分母lnv=(dec/nec)^((M-k+1)*N);%计算第一项中对数的自变量AIC(k)=-2*log(lnv)+2*(k-1)*(2*M-k+1);%计算AIC准则end[Amin,K]=min(AIC);N1=K-1;%信号源个数估计%根据MDL准则进行信号源个数的估计fork=1:Mdec=prod(ev(k:M).^(1/(M-k+1)));nec=mean(ev(k:M));lnv=(dec/nec)^((M-k+1)*N);MDL(k)=-log(lnv)+(k-1)/2*(2*M-k+1)*log(N);%计算DML准则end[Amin,K]=min(MDL);%寻找使DML准则最小的索引N2=K-1;%信号源个数的估计（2）根据（1）中信源个数画出相应的MUSIC频率估计谱线。%计算MUSIC谱En=U(:,N1+1:M);%噪声子空间的向量组成的矩阵NF=2048;forn=1:NFAq=exp(-j*2*pi*(n-1)/NF*(0:M-1)');Pmusic(n)=1/(Aq'*En*En'*Aq);%music谱Pmusic(n)=10*log10(Pmusic(n));endf=-0.5:1/NF:0.5-1/NF;plot(f,Pmusic);图7MUSIC谱估计3.20复正弦加白噪声随机过程()un同题3.18中所给。试使用Root-MUSIC方法进行信号频率估计的仿真实验。（要求：信号样本数取1000，估计的自相关矩阵为8阶。）（1）计算正弦信号的频率估计值。N=1000;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);signal1=exp(j*0.5*pi*(0:N-1)+j*2*pi*rand);signal2=exp(-j*0.3*pi*(0:N-1)+j*2*pi*rand);un=signal1+signal2+noise;M=8;fork=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=xs*xs'/(N-M);[U,E]=svd(R);ev=diag(E);G=U(:,3:M);Gr=G*G';co=zeros(2*M-1,1);form=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);输出结果为：最接近单位圆的两个根是：-0.000102126882707680+1.00226420141169i-0.000101665975924671+0.997740903255518i上述根的归一化频率为：0.2500162172789630.250016217278959（2）与MUSIC算法的估计结果比较。En=G;NF=2048;forn=1:NFAq=exp(-j*2*pi*(n-1)/NF*(0:M-1)');Pmusic(n)=1/(Aq'*En*En'*Aq);Pmusic(n)=10*log10(Pmusic(n));endf=-0.5:1/NF:0.5-1/NF;plot(f,Pmusic);图8MUSIC谱估计3.21复正弦加白噪声随机过程()un同题3.18中所给。试使用ESPRIT算法进行信号频率估计的仿真实验，给出正弦信号频率的估计值（要求：信号样本数取1000，估计的自相关矩阵为8阶。）解：N=1000;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);signal1=exp(j*0.5*pi*(0:N-1)+j*2*pi*rand);signal2=exp(-j*0.3*pi*(0:N-1)+j*2*pi*rand);un=signal1+sign