现代数字信号处理仿真作业

现代数字信号处理仿真作业

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-0i a9=2+0i; a10=2+3i; a11=7-10i; a12=4-9i; a13=8-3i ;

a14=2+4i; a15=2+1i; a16=3i.

仿真程序（3_17）: clearall clc

%%产生噪声序列

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);

form=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); fori=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算法中最接近单位圆的两个根为： + -

对应的归一化频率为：

相同信号的MUSIC谱估计结果如下

图对3.20信号进行MUSIC谱估计的结果

仿真程序（3_20）: clearall