现代数字信号处理理论计算法仿真作业

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

3.17

首先定义函数 u_length(N)，函数的作用是产生长度为N的观测样本u(n),程序如下：

function u = length(N); A1 = sqrt(10^(30/10)); A2 = sqrt(10^(30/10)); A3 = sqrt(10^(27/10)); n = 1:N;

s_signal = A1*exp(i*2*pi*0.15*n)+A2*exp(i*2*pi*0.17*n)+A2*exp(i*2*pi*0.26*n); v_noise = sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N);; u = s_signal+v_noise;

(1)N=32，分别用给予FFT的自相关函数和式3.1.2估计自相关函数r(m),程序如下：

function r11 N=32;

u = u_length(N); %长度为32的随机信号 for m = 1:N

r11(m) = u(m:N)*u(1:N-m+1)'/N; end

r12 = conj(r11(N:-1:2)); %自相关函数的对称性

r1 = [r12,r11]; %r1是按式3.1.2估计自相关函数r(m) u2N = [u,zeros(1,N)]; U_w = fft(u2N,2*N); S_BT = (abs(U_w).^2)/N; r20=ifft(S_BT);

r2=[r20(N+2:2*N) r20(1:N)]; %r2是FFT的自相关函数快速算法得出的自相关函数 hold on

k = -N+1:N-1;

plot(k,abs(r1),'r'); plot(k,abs(r2),'g*'); grid on 仿真图如下：

（2）N=256，分别用BT法和周期图发估计u(n)的功率谱 %BT法估计功率谱密度，自相关采用M=64点

N = 256; M=64;

u = u_length(N);

u_xcorr = xcorr(u,M,'biased'); m = -M:M; k = -100:100;

h = u_xcorr*exp(-j*pi/100).^(m'*k); Sw_BT = 10*log10(abs(h)./max(abs(h))); plot(k/200,Sw_BT);

title('BT法估计功率谱密度')； grid on

BT法估计功率谱密度仿真结果如下图：

%周期图法估计功率谱密度

Uw_N1=abs(fft(u,N)); Uw_N=fftshift(Uw_N1); w=1:N;

plot(w/N-0.5,10*log10(abs(Uw_N)./max(abs(Uw_N)))); title('周期图法估计结果')

周期图法估计功率谱密度仿真结果如下图：

（3）L—D迭代算法求解AR模型参数及功率谱估计，程序如下：

N=256; p=16;

u = u_length(N)

r=xcorr(u,p,'biased'); r=r(p+1:2*p+1); a=-r(2)/r(1);

c=r(1)-r(2)*r(2)'/r(1); for m=2:p

l=conj(r(m:-1:2)'); k=-(r(m+1)+a*l)/c; for n=1:m-1

b(n)=a(n)+k*conj(a(m-n)); end

a=[b k];

c=c*(1-k*k'); end

h=freqz(sqrt(c),[1 a],201,'whole'); h=fftshift(h); b=abs(h).^2; n=-100:100;

plot(n/200,10*log10(b/max(b))); grid

title('阶数为16的AR模型功率谱密度估计曲线') 仿真结果如下图：

3.20

clear; N=1000; M=8;K=2;

v=sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N); n=1:N;

fi=2*pi*rand(1,2);

u(n)=exp(j*0.5*pi*n+j*fi(1))+exp(j*1.2*pi*n+j*fi(2))+v(n); r=xcorr(u,M,'biased'); for n=1:M for m=1:M

R(m,n)=r(M+1+n-m); end end

[G,D]=eigs(R,M-K,'sm'); for n=-100:100;

for m=1:M

a(m)=exp(j*(m-1)*pi*n/100); end

p(n+101)=1/(a*(G*G')*a');

end

p=abs(p); n=-1:0.01:1;

figure;plot( n,10*log10(p/max(p)) );title(??MUSIC??·¨?á1??ˉ) Gr=G*G';

a=zeros(2*M-1,1); for m=1:M

a(m:m+M-1)=a(m:m+M-1)+Gr(M:-1:1,m); end

x=roots(a);

err=abs(abs(x)-1); w=angle(x)/pi; e=sort(err); for n=1:2*K

t(n)=find(err==e(n)); end

figure;bar(w(t)); grid

root—music 估计仿真结果：

MUSIC算法频率估计结果：

本文来源：https://www.bwwdw.com/article/57mh.html