信号处理仿真题作业.doc

重庆邮电大学信号处理基础仿真作业 学号：S120101057 姓名：贾雪婷3.17在计算机上用如下方法产生随机信号的观测样本：首先产生一段零均值、方差为的复高斯白噪声序列；然后在上叠加三个复正弦信号，它们的归一化频率分别是f1=0.15，f2=0.17和f3=0.26。调整和正弦信号的幅度，使在f1、f2和f3处得信噪比分别为30dB、30dB和27dB。（1）令信号观测样本长度N=32，试用3.1.1节讨论的基于FFT的自相关函数快速计算方法估计出自相关函数，并与教材式（3.1.2）估计出的自相关函数做比较。产生零均值、方差为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法和周期图法估计的功率谱，这里设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)图 4 BT法所得功率谱（3）令信号观测样本长度N=256，试用Levinson-Durbin迭代算法求解AR模型的系数并估计的功率谱，模型的阶数取为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);for m=2:p k(m)= -(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)/sigma(m-1); a(m,m)=k(m); for i=1:m-1 a(m, i)=a(m-1, i) + k(m) * conj(a(m-1, m-i); end sigma(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)图 5 16阶AR模型的功率谱估计3.18设随机过程为其中，是零均值、方差为1的白噪声，、是相互独立并在上服从均匀分布的随机相位。使用MVDR方法进行信号频率估计的方针实验，画出频率估计谱线，并给出正弦信号频率的估计值。（要求：信号样本数取1000，估计的自相关矩阵为8阶）解：clear all;%产生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);for m=1:N for n=m:N r(m)=r(m)+un(n)*(un(n-(m-1)'end r(m)=r(m)/N;endM=8;%设定横向滤波器的阶数%利用自相关矩阵和自相关函数的关系，构建自相关矩阵R=zeros(M,M);for i=1:M for j=1:M if i<=j R(i,j)=r(1+j-i);%得到了M阶的矩阵R else R(i,j)=(r(1+i-j)' end end endN3=2048;%设定画图时描点的数目。d=1/(N3-1);%求画图用的横坐标的间隔。h=zeros(1,N3);for i=1:N3 h(i)=-0.5+(i-1)*d;endy=zeros(1,N3);for j=1:N3 w=h(j)*2*pi; y(j)=10*log10(abs(getPMVDR(w,M,R);endplot(h,y);title('MVDR谱估计');xlabel('omega/2pi');ylabel('归一化功率谱/dB');axis(-0.5 0.5 -16 1); 图 6 MVDR谱估计3.19复正弦加白噪声随机过程同题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;for k = 1 : N-M xs(:,k)=un(k+M-1:-1:k).' endR=xs*xs'/(N-M);% 自相关矩阵的特征值分解U,E=svd(R); % U是特征矢量组成的矩阵，E是对角阵，对角元素是由大到小排列的特征值ev=diag(E); %提取对角元素上的特征值； % 根据AIC准则进行信号源个数的估计 for k=1: M dec = 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准则 endAmin, K=min(AIC);N1=K-1; % 信号源个数估计% 根据MDL准则进行信号源个数的估计 for k=1:M dec = 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准则 endAmin, K=min(MDL); %寻找使DML准则最小的索引N2=K-1; %信号源个数的估计（2）根据（1）中信源个数画出相应的MUSIC频率估计谱线。%计算MUSIC谱En=U(:,N1+1:M); %噪声子空间的向量组成的矩阵NF=2048;for n=1:NF Aq=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);图 7 MUSIC谱估计3.20复正弦加白噪声随机过程同题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;for k = 1 : N-M xs(:,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);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);输出结果为：最接近单位圆的两个根是：-0.000102126882707680 + 1.00226420141169i-0.000101665975924671 + 0.997740903255518i上述根的归一化频率为：0.2500162172789630.250016217278959（2）与MUSIC算法的估计结果比较。En=G; NF=2048;for n=1:NF Aq=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);图 8 MUSIC谱估计3.21复正弦加白噪声随机过程同题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+signal2+noise;M=8;for k=1:N-M xs(:,k)=un(k+M-1:-1:k).'endRxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-M-1); %计算自相关矩阵RxxRxy=xs(:,1:end-1)*xs(:,2:end)'/(N-M-1); %计算互相关矩阵Rxy%相关矩阵的特征值分解，寻找最小特征值%U,E=svd(Rxx); %矩阵U是特征矢量组成的矩阵，E是对角阵，对角元素是由大到小排列的特征值ev=diag(E); %提取对角元素上的特征值；emin=ev(end); %获取最小特征值；%构造矩阵对Cxx,Cxy%Z=zeros(M-1,1),eye(M-1);0,zeros(1,M-1); %构造矩阵Z;Cxx=Rxx-emin*eye(M); %计算Cxx;Cxy=Rxy-emin*Z; %计算Cxy;%矩阵对Cxx,Cxy的广义特征值分解%U,E=eig(Cxx,Cxy); %广义特征值分解z=diag(E); %提取对角矩阵中的特征值ph=angle(z)/(2*pi); %求所有特征值的相位对应的归一化频率err=abs(abs(z)-1); %与单位元的距离输出结果为：单次ESPRIT算法中最接近单位圆的两个特征值为:0.588516696017445 - 0.809433346373638i-0.00139435925779186 + 0.999262092748743i上述特征值的相位对应的归一化频率为：-0.1499448116003380.25022208290087311