实用标准文档
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-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); end
Par=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.001503153449793i 0.587376604261220 - 0.810845628739986i 对应的归一化频率为： 0.250183714447964 -0.150223406926494
相同信号的MUSIC谱估计结果如下
图 5 对3.20信号进行MUSIC谱估计的结果
仿真程序（3_20）:
clear all clc
%% 信号样本和高斯白噪声的产生 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-M
xs(:,k)=un(k+M-1:-1:k).'; end
R=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:M
co(m:m+M-1)=co(m:m+M-1)+Gr(M:-1:1,m); end z=roots(co); ph=angle(z)/(2*pi); err=abs(abs(z)-1); %% 计算MUSIC谱 En=U(:,2+1:M); NF=2048; for n=1:NF
Aq=exp(-1i*2*pi*(-0.5+(n-1)/(NF-1))*(0:M-1)'); Pmusic(n)=1/(Aq'*En*En'*Aq); end
kk=-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算法中最接近单位元的两个特征值为： 0.001826505974929 + 1.000690248438859i 0.586994191014025 - 0.809491260856630i 对应的归一化频率为： 0.249709503383161 -0.150146235268272
仿真程序（3_21）:
clear all clc
%% 信号样本和高斯白噪声的产生 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-M
xs(:,k)=un(k+M-1:-1:k).'; end
Rxx=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);
文案大全
%复正弦信号