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

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

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）:

clear all

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

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算法中最接近单位圆的两个根为：

+

-

对应的归一化频率为：

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

图对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算法中最接近单位元的两个特征值为：

+

-

对应的归一化频率为：

仿真程序（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);

4.仿真题4.18

仿真结果及图形：

步长为0.05时失调参数为m1=0.0493；

步长为0.005时失调参数为m2=0.0047。

图步长为0.05时权向量的收敛曲线

图步长为0.005时权向量的收敛曲线

图步长分别为0.05和0.005时100次独立实验的学习曲线仿真程序（4_18）:

clear all

clc

%%产生100组独立样本序列

data_len=512;

trials=100;

n=1:data_len;

a1=-0.975;

a2=0.95;

sigma_v_2=0.0731;

v=sqrt(sigma_v_2)*randn(data_len,1,trials);

u0=[00];

num=1;

den=[1,a1,a2];

Zi=filtic(num,den,u0);

u=filter(num,den,v,Zi);%产生100组独立信号

%%LMS迭代

mu1=0.05;mu2=0.005;

w1=zeros(2,data_len,trials);w2=w1;

for m=1:100;

temp=zeros(data_len,1);

e1(:,:,m)=temp;e2(:,:,m)=temp;d1(:,:,m)=temp;d2(:,:,m)=temp;

for n=3:data_len-1

w1(:,n+1,m)=w1(:,n,m)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n,1,m));

w2(:,n+1,m)=w2(:,n,m)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n,1,m));

d1(n+1,1,m)=w1(:,n+1,m)'*u(n:-1:n-1,:,m);

d2(n+1,1,m)=w2(:,n+1,m)'*u(n:-1:n-1,:,m);

e1(n+1,1,m)=u(n+1,:,m)-d1(n+1,1,m);

e2(n+1,1,m)=u(n+1,:,m)-d2(n+1,1,m);

end

end

t=1:data_len;

w1_mean=zeros(2,data_len);w2_mean=w1_mean;e1_mean=zeros(data_len,1);e2_mean=e 1_mean;

for m=1:100

w1_mean=w1_mean+w1(:,:,m);

w2_mean=w2_mean+w2(:,:,m);

e1_mean=e1_mean+e1(:,:,m).^2;

e2_mean=e2_mean+e2(:,:,m).^2;

end

w1_mean=w1_mean/100;%100次独立实验权向量的均值

w2_mean=w2_mean/100;

e1_100trials_ave=e1_mean/100;%100次独立实验的学习曲线均值

e2_100trials_ave=e2_mean/100;

figure(1)

plot(t,w1(1,:,90),t,w1(2,:,90),t,w1_mean(1,:),t,w1_mean(2,:))

xlabel('迭代次数');ylabel('权向量')

title('步长=0.05')

figure(2)

plot(t,w2(1,:,90),t,w2(2,:,90),t,w2_mean(1,:),t,w2_mean(2,:))

xlabel('迭代次数');ylabel('权向量')

title('步长=0.005')

%%计算剩余误差和失调参数

wopt=zeros(2,trials);

Jmin=zeros(1,trials);

sum_eig=zeros(trials,1);

for m=1:trials

rm=xcorr(u(:,:,m),'biased');

R=[rm(512),rm(513);rm(511),rm(512)];

p=[rm(511);rm(510)];

wopt(:,m)=R\p;

[v,d]=eig(R);

Jmin(m)=rm(512)-p'*wopt(:,m);

sum_eig(m)=d(1,1)+d(2,2);

end

sJmin=sum(Jmin)/trials;

Jex1=e1_100trials_ave-sJmin;%剩余均方误差mu1

Jex2=e2_100trials_ave-sJmin;%剩余均方误差mu2

sum_eig_100trials=sum(sum_eig)/100;

Jexfin1=mu1*sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials));

Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials));

M1=Jexfin1/sJmin;%失调参数m1

M2=Jexfin2/sJmin;%失调参数m2

figure(3)

plot(t,e1_100trials_ave,'*',t,e2_100trials_ave)

xlabel('迭代次数');ylabel('均方误差')

legend('u1=0.05','u2=0.005')

axis([0,600,0,1])

5.仿真题5.10

仿真结果及图形：

(1)M=2时，210.99,0.93627a σ=-=，求解Yule-Walker 方程

可得到自相关矩阵

相应的计算程序为

r2=inv([1,-0.99;-0.99,1])*[0.93627;0];

R2=[r2(1),r2(2);r2(2),r2(1)];%M=2

(2)M=3时，2120.99,0,0.93627a a σ=-==，求解Yule-Walker 方程

可得到自相关矩阵为

相应的计算程序为

r3=inv([1,-0.99,0;-0.99,1,0;0,-0.99,1])*[0.93627;0;0];

R3=[r3(1),r3(2),r3(3);r3(2),r3(1),r3(2);r3(3),r3(2),r3(1)];%M=3

(3)计算特征值扩展

%%特征值分解

eig_value_1=eig(R2);

eig_value_2=eig(R3);

%%特征值扩展

eig_spread_1=max(eig_value_1)/min(eig_value_1);

eig_spread_2=max(eig_value_2)/min(eig_value_2);

M=2时特征值扩展是199.0000；

M=3时特征值扩展是444.2790。

(3)根据LMS 算法均方误差收敛特性，M=2时步长因子应在区间（0,0.0213）中，M=3时，步长因子应在区间（0,0.0142）之间，因此题中的步长因子不合理。故在仿真中，M=2时采用步长因子0.001，M=3时采用步长因子0.0006.

图500次独立实验M=2步长为0.001时权向量收敛曲线

图500次独立实验M=3步长为0.0006时权向量收敛曲线

图500次独立实验M=2步长为0.001时的学习曲线

图500次独立实验M=3步长为0.0006时的学习曲线仿真程序（5_10_4）:

clear all

clc

%%产生系统输入白噪声

L=10000;

sigma_v1_2=0.93627;

for m=1:500

v(:,m)=sqrt(sigma_v1_2)*randn(L,1);

end

%%生成500组独立的AR模型信号

a1=-0.99;

for m=1:500

u(1,1,m)=v(1,m);

for k=2:L

u(k,1,m)=-a1*u(k-1,1,m)+v(k,m);

end

end

%%LMS迭代算法

M=2;

%M=3;

mu=0.001;

%mu=0.0006;

w=zeros(L,M,500);

for m=1:500

e(1,m)=u(1,m);

uu=zeros(1,M);

w(2,:,m)=w(1,:,m)+mu*e(1,m)*uu;

uu=[u(1,m)uu(1:M-1)];

dd=(w(2,:,m)*uu')';

e(2,m)=u(3,m)-dd;

for k=3:L

w(k,:,m)=w(k-1,:,m)+mu*e(k-1,m)*uu;

uu=[u(k-1,1,m)uu(1:M-1)];

dd=(w(k,:,m)*uu')';

e(k,m)=u(k,m)-dd;

end

end

%%M=2

e_mean=zeros(10000,1);

w_mean=zeros(10000,2);

for m=1:500

w_mean=w_mean+w(:,:,m);

e_mean=e_mean+e(:,m).^2;

end

w_mean=w_mean/500;

e_mean=e_mean/500;

t=1:L;

figure(1)

plot(t,w(:,1,100),t,w(:,2,100),t,w_mean(:,1),t,w_mean(:,2))

xlabel('迭代次数n');ylabel('抽头权值');

title('M=2,步长0.001的权向量收敛曲线')

figure(2)

plot(t,e_mean)

xlabel('迭代次数n');ylabel('MSE');title('M=2,步长0.001的学习曲线')

%%M=3

e_mean=zeros(10000,1);

w_mean=zeros(10000,3);

for m=1:500

w_mean=w_mean+w(:,:,m);

e_mean=e_mean+e(:,m).^2;

end

w_mean=w_mean/500;

e_mean=e_mean/500;

t=1:L;

figure(1)

plot(t,w(:,1,100),t,w(:,2,100),t,w(:,3,100),t,w_mean(:,1),t,w_mean(:,2),t,w_m ean(:,3))

xlabel('迭代次数n');ylabel('抽头权值');

title('M=3,步长0.0006的权向量收敛曲线')

figure(2)

plot(t,e_mean)

xlabel('迭代次数n');ylabel('MSE');title('M=2,步长0.0006的学习曲线')

6.仿真题6.13

仿真结果及图形：

滤波器抽头个数为4时

图M=4时的MVDR谱

图M=4时基于奇异值分解的MVDR谱

从上面两图可以看出，M=4时并没有将3个频点分辨出来，增加滤波器阶数可以解决此问题，因此当M=20时仿真结果如下两图所示：

图M=20时的MVDR谱

图M=20时基于奇异值分解的MVDR谱

仿真程序（6_13）:

clear all

clc

%%产生观测信号

M=4;

%M=20;

N=1000;

f=[0.10.250.27];

SNR=[303027];

sigma=1;

Am=sqrt(sigma*10.^(SNR/10));

t=linspace(0,1,N);

phi=2*pi*rand(size(f));

vn=sqrt(sigma/2)*randn(size(t))+1i*sqrt(sigma/2)*randn(size(t)); Un=vn;

for k=1:length(f)

s=Am(k)*exp(1i*2*pi*N*f(k).*t+1i*phi(k));

Un=Un+s;

end

Un=Un.';

%%构建矩阵

A=zeros(M,N-M+1);

for n=1:N-M+1

A(:,n)=Un(M+n-1:-1:n);

end

[U,S,V]=svd(A');

invphi=V*inv(S'*S)*V';

%%构建向量

P=1024;

f=linspace(-0.5,0.5,P);

omega=2*pi*f;

a=zeros(M,P);

for k=1:P

for m=1:M

a(m,k)=exp(-1i*omega(k)*(m-1));

end

end

%%计算MVDR谱

Pmvdr=zeros(size(omega));

for k=1:P

Pmvdr(k)=1/(a(:,k)'*invphi*a(:,k));

end

Pmvdr=abs(Pmvdr/max(abs(Pmvdr)));

Pmvdr=10*log10(Pmvdr);

kk=-0.5+(0:P-1)*(1/(P-1));

figure(1)

plot(kk,Pmvdr)

%%基于习题6.11的奇异值分解的MVDR方法

for k=1:P

Sw=zeros(1,M);

for i=1:M

Sw(i)=(a(:,k)')*V(:,i)/S(i,i);

end

P_svd(k)=1/sum((abs(Sw)).^2);

end

P_svd=abs(P_svd/max(abs(P_svd)));

P_svd=10*log10(P_svd);

xlabel('w/2*pi');ylabel('归一化频谱/dB')

title('M=4的MVDR谱')

%title('M=20的MVDR谱')

figure(2)

plot(kk,P_svd)

xlabel('w/2*pi');ylabel('归一化频谱/dB')

title('M=4的基于SVD的MVDR频谱')

%title('M=20的基于SVD的MVDR频谱')

7.仿真题6.15

仿真结果及图形：

图单次实验估计权值以及500次独立实验的估计权值

图500次独立实验的学习曲线

仿真程序（6_15）:

clear all

clc

%%产生AR模型的输入信号

a1=0.99;

sigma=0.995;

N=1000;

trials=500;

vn=sqrt(sigma)*randn(N,1,trials);

nume=1;

deno=[1a1];

u0=zeros(length(deno)-1,1);

xic=filtic(nume,deno,u0);

un=filter(nume,deno,vn,xic);%产生500组独立的信号

%%产生期望信号和观测数据矩阵

n0=1;

M=2;

b=un(n0+1:N,:,:);

L=length(b);

A=zeros(M,L,trials);

for m=1:trials

un1=[zeros(M-1,1).',un(:,:,m).'].';

for k=1:L

A(:,k,m)=un1(M-1+k:-1:k);

end

end

%%RLS算法求最优权向量

delta=0.004;

lambda=0.98;

w=zeros(M,L+1,trials);

for m=1:trials

epsilon=zeros(L,1);

P1=eye(M)/delta;

for k=1:L

PIn=P1*A(:,k,m);

denok=lambda+A(:,k,m)'*PIn;

kn=PIn/denok;

epsilon(k)=b(k,1,m)-w(:,k,m)'*A(:,k,m);

w(:,k+1,m)=w(:,k,m)+kn*conj(epsilon(k));

P1=P1/lambda-kn*A(:,k,m)'*P1/lambda;

end

MSE(m,:)=abs(epsilon).^2;

end

MSE_mean=zeros(1,L);w_mean=zeros(2,N);

for m=1:trials

w_mean=w_mean+w(:,:,m);

MSE_mean=MSE_mean+MSE(m,:);

end

w_mean=w_mean/trials;%500次独立实验权向量的均值

MSE_mean=MSE_mean/trials;%500次独立实验的MSE

t=1:1000;

figure(1)

plot(t,w(1,:,100),t,w(2,:,100),t,w_mean(1,:),t,w_mean(2,:)) xlabel('迭代次数n');ylabel('权值')

figure(2)

plot(1:L,MSE_mean)

xlabel('迭代次数n');ylabel('MSE')

8.仿真题6.15

仿真结果及图形：

图100次独立实验权值变化曲线

图单次实验权值变化曲线

图100次独立实验MSE取对数的变化曲线仿真程序（7_14）:

clear all

clc

%%产生白噪声

N=2000;

gv=0.0332;

for m=1:100

v(m,:)=randn(1,N)*sqrt(gv);

end

%%产生AR模型序列

a=[1.6-1.460.616-0.1525];

u1=zeros(1,N,100);

for m=1:100%产生100组独立实验信号

for i=1:(N-4)

u1(1,i+4,m)=a(1)*u1(1,i+3,m)+a(2)*u1(1,i+2,m)+a(3)*u1(1,i+1,m)+a(4)*u1(1,i,m) +v(m,i+4);

end

end

%%卡尔曼滤波¨

N2=2000;

Jmin=0.005;

for m=1:100

for i=5:N2

U(:,i,m)=[u1(1,i-1,m);u1(1,i-2,m);u1(1,i-3,m);u1(1,i-4,m)];

end

end

W_esti=zeros(4,N2,100);

for m=1:100

P_esti=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];

for l=1:N2-1

P_pre=P_esti;

A=(U(:,l,m))'*P_pre*U(:,l,m)+Jmin;

K=P_pre*U(:,l,m)/A;

alpha(l)=u1(1,l,m)-(U(:,l,m))'*W_esti(:,l,m);

W_esti(:,l+1,m)=W_esti(:,l,m)+K*alpha(l);

P_esti=P_pre-K*(U(:,l,m))'*P_pre;

end

end

w=zeros(4,N2);e=zeros(1,N2,100);d=zeros(1,N2,100);MSE=zeros(1,N2);

for m=1:100

w=w+W_esti(:,:,m);

for n=5:N2

d(1,n,m)=W_esti(:,n,m)'*u1(1,n-1:-1:n-4,m)';

e(1,n,m)=u1(1,n,m)-d(1,n,m);

end

MSE=MSE+e(:,:,m).^2;

end

w=w/100;%100次独立实验的权向量均值

MSE=MSE/100;%100次独立实验的均方误差

t=1:N2;

figure(1)

plot(t,w(1,:),t,w(2,:),t,w(3,:),t,w(4,:))

legend('w1','w2','w3','w4')

xlabel('迭代次数');ylabel('权值')

figure(2)

plot(t,W_esti(1,:,50),t,W_esti(2,:,50),t,W_esti(3,:,50),t,W_esti(4,:,50)) legend('w1','w2','w3','w4');xlabel('迭代次数');ylabel('权值')

figure(3)

semilogy(t,MSE)

xlabel('迭代次数');ylabel('对数MSE')

9.仿真题8.16

仿真结果及图形：

单次RootMUSIC算法得到的DOA估计为：

-9.993839.9904

单次ESPRIT算法得到的DOA估计为：

-9.971639.7934

图MUSIC算法实现DOA估计

图MVDR算法实现DOA估计仿真程序（8_16）:

clear all

clc

%%产生阵列接收信号

N=100;

M=10;

K=2;

theta=[-10;40]*pi/180;

SNR=[10;20];

Am=[sqrt(10.^(SNR/10))];

S=Am*ones(1,N);

S(2,:)=S(2,:).*exp(1i*2*pi*rand(1,N));

for a=1:M

for b=1:K

A(a,b)=exp(-1i*(a-1)*pi*sin(theta(b)));

end

end

V=zeros(M,N);

for m=1:M

v=wgn(1,N,0,'complex');

v=v-mean(v);

v=v/std(v);

V(m,:)=v;

end

X=A*S+V;

%%估计空间自相关矩阵

R=zeros(M,M);

for i=1:N

R=R+X(:,i)*X(:,i)';

end

R=R/N;

%%MUSIC算法·¨

[VR,D]=eig(R);

[B,IX]=sort(diag(D));

G=VR(:,IX(M-K:-1:1));

Pmusic=[];

for n=-pi/2:pi/180:pi/2

a=exp(-1i*[0:M-1]'*pi*sin(n));

Pmusic=[Pmusic,1/(a'*G*G'*a)];

end

%%MVDR算法

Pmvdr=[];

for n=-pi/2:pi/180:pi/2

a=exp(-1i*[0:M-1]'*pi*sin(n));

Pmvdr=[Pmvdr,1/(a'*inv(R)*a)];

end

%%RootMUSIC算法¨

syms z

pz=z.^([0:M-1]');

pz1=(z^(-1)).^([0:M-1]);

fz=z^(M-1)*pz1*G*G'*pz;

a=sym2poly(fz);

r=roots(a);

r1=abs(r);

for i=1:2*K

[Y,I(i)]=min(abs(r1-1));

r1(I(i))=inf;

end

for i=1:2*K

theta_esti_music(i)=asin(-angle(r(I(i)))/pi)*180/pi;

end

%%ESPRIT算法

S=VR(:,IX(M:-1:M-K+1));

S1=S(1:M-1,:);

S2=S(2:M,:);

fai=S1/S2;

[U_fai,V_fai]=eig(fai);

for i=M-K:M-1

theta_esti_esprit(M-i)=asin(angle(V_fai(i,i))/pi)*180/pi; end

t=-90:90;

Pmusic_norm=10*log10(abs(Pmusic)/max(abs(Pmusic)));

Pmvdr_norm=10*log10(abs(Pmvdr)/max(abs(Pmvdr)));

figure(1)

plot(t,Pmusic_norm)

xlabel('空间角度/(单位度)');ylabel('归一化空间谱/dB')

title('MUSIC')

figure(2)

plot(t,Pmvdr_norm)

xlabel('空间角度/(单位度)');ylabel('归一化空间谱/dB')

title('MVDR')

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

现代数字信号处理仿真作业 1.仿真题3.17 仿真结果及图形： 图 1 基于FFT的自相关函数计算

图 3 周期图法和BT 法估计信号的功率谱 图 2 基于式3.1.2的自相关函数的计算

图 4 利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为： a1=-0.402637623107952-0.919787323662670i; a2=-0.013530139693503+0.024214641171318i; a3=-0.074241889634714-0.088834852915013i; a4=0.027881022353997-0.040734794506749i; a5=0.042128517350786+0.068932699075038i; a6=-0.0042799971761507 + 0.028686095385146i; a7=-0.048427890183189 - 0.019713457742372i; a8=0.0028768633718672 - 0.047990801912420i a9=0.023971346213842+ 0.046436389191530i; a10=0.026025963987732 + 0.046882756497113i; a11= -0.033929397784767 - 0.0053437929619510i; a12=0.0082735406293574 - 0.016133618316269i; a13=0.031893903622978 - 0.013709547028453i ; a14=0.0099274520678052 + 0.022233240051564i; a15=-0.0064643069578642 + 0.014130696335881i; a16=-0.061704614407581- 0.077423818476583i. 仿真程序（3_17）: clear all clc %% 产生噪声序列 N=32; %基于FFT的样本长度

现代数字信号处理及其应用——LMS算法结果及分析

LMS 算法MATLAB 实现结果及其分析 一、LMS ：为课本155页例题 图1.1：LMS 算法学习曲线（初始权向量[]T 00w ?=） 图1.2滤波器权系数迭代更新过程曲线（步长075.0=μ） 图1.3滤波器权系数迭代更新过程曲线（步长025.0=μ）图1.4滤波器权系数迭代更新过程曲线（步长015.0=μ） 分析解释： 在图1.1中，收敛速度最慢的是步长为015.0=μ的曲线，收敛速度最快的是步长075.0=μ的曲线，所以可以看出LMS 算法的收敛速度随着步长参数的减小而相应变慢。图1.2、1.3、1.4分别给出了步长为075.0=μ、025.0=μ、025.0=μ的滤波器权系数迭代更新过程曲线，可以发现其不是平滑的过程，跟最抖下降法不一样，体现了其权向量是一个随机过程向量。

LMS2：为课本155页例题，156页图显示结果 图2.1：LMS 算法学习曲线（初始权向量[]T 00w ?=） 图2.2滤波器权系数迭代更新过程曲线（步长025.0=μ） 图2.3滤波器权系数迭代更新过程曲线（步长025.0=μ）图2.4最陡下降法权值变化曲线（步长025.0=μ） 分析解释： 图2.1给出了步长为025.0=μ的学习曲线，图2.2给出了滤波器权向量的单次迭代结果。图2.3给出了一 次典型实验中所得到的权向量估计()n w ?=，以及500次独立实验得到的平均权向量()}n w ?E{=的估计，即()∑==T t n w T 1 t )(?1n w ?，其中)(?n w t 是第t 次独立实验中第n 次迭代得到的权向量，T 是独立实验次数。可以发现，多次独立实验得到的平均权向量()}n w ?E{=的估计平滑了随机梯度引入的梯度噪声，使得其结果与使用最陡下降法（图2.4）得到的权向量趋于一致，十分接近理论最优权向量[]T 7853.08361.0w 0-=。 LMS3：为课本172页习题答案

数字信号处理实验七小信号放大器特性分析与仿真

实验七小信号放大器特性分析与仿真1，实验目的 使用matlab分析各种小信号放大器的结构、参数及特性，加深对各种小信号放大器的理解和认识 二、实验原理 小信号放大器是电子线路的重要组成部分之一，由于他工作在晶体管的线性区域之内，因此又称为线性放大器。使用MATLAB可仿真小信号放大器的各种参数，如电压增益，输入阻抗，输出阻抗，频率响应等等。 1、晶体三极管的等效电路 常见的晶体三极管等效电路有：低频h参数，共基极T型高频等效电路，混合π型高频等效电路，他们通常用于分析各种小信号晶体管放大器的特性。 共发射极h参数的等效电路如图(a)所示，它适用于对低频放大器进行分析。另外，还存在着一种简化的h参数等效电路，其中忽略晶体管内部的电压反馈系数。共发射极的h参数与各电压电流的关系为。 共基极T型高频等效电路如图(b)所示，适用于共基极高频放大电路进行分析，工作频率可达100MHZ以上。 混合π型高频等效电路如图(c)所示，适用于分析共发射极的高频发达电路。在较宽的频率范围之内，等效电路的参数和工作频率无关。另外还存在着简化的混合π型高频等效电路，其中和处于开路状态。 2、共发射极放大电路 共发射极放大电路是一种使用的最为广泛的放大电路形式，其特点是电压增益和电流增益都比较高。自定义M函数amplifl..m用来仿真共发射极放大电路，使用它可以计算该放大器的的智力参数和交流参数。该

放大器的电路如下图。 MATLAB的特点之一就是适合进行线性代数运算，因此午在分析直流参数或分析交流参数时，都可以采用基尔霍夫定理，然后采用矩阵求逆的方式求出电压和电流的具体数值，进一步便可得到该放大器的各种参数。在分析共发射极放大的交流参数时，采用的晶体管模型是低频H 参数等效电路。一般来说，每个晶体管都可以用三个节点来表示，他们分别是基极集电极和发射极。在计算交流参数过程中，忽略各电容器的容抗。 3、直接耦合放大器 在两个或三个晶体管之间进行直接耦合的放大器称为直接耦合放大器，他多用作音响系统中的前置放大器，录音机内的磁头放大器。直接耦合放大器的主要特点是工作点稳定，电压增益高，下图是一个典型的直接耦合放大电路，它有三个晶体管构成，第一级为低噪声放大，第二级为高增益放大，第三极为射随器，整个放大器的电压增益由负反馈电路确定。由于采用了串联电压负反馈，同时又使用了射随器，因此该电路具有较高的输入阻抗和较低的输出阻抗。 4、差分放大器

数字信号处理实验作业

实验6 数字滤波器的网络结构 一、实验目的： 1、加深对数字滤波器分类与结构的了解。 2、明确数字滤波器的基本结构及其相互间的转换方法。 3、掌握用MA TLAB 语言进行数字滤波器结构间相互转换的子函数及程序编写方法。 二、实验原理： 1、数字滤波器的分类 离散LSI 系统对信号的响应过程实际上就是对信号进行滤波的过程。因此，离散LSI 系统又称为数字滤波器。 数字滤波器从滤波功能上可以分为低通、高通、带通、带阻以及全通滤波器；根据单位脉冲响应的特性，又可以分为有限长单位脉冲响应滤波器（FIR ）和无限长单位脉冲响应滤波器（IIR ）。 一个离散LSI 系统可以用系统函数来表示： M -m -1-2-m m m=0 012m N -1-2-k -k 12k k k=1 b z b +b z +b z ++b z Y(z)b(z)H(z)=== =X(z)a(z) 1+a z +a z ++a z 1+a z ∑∑ 也可以用差分方程来表示： N M k m k=1 m=0 y(n)+a y(n-k)=b x(n-m)∑∑ 以上两个公式中，当a k 至少有一个不为0时，则在有限Z 平面上存在极点，表达的是以一个IIR 数字滤波器；当a k 全都为0时，系统不存在极点，表达的是一个FIR 数字滤波器。FIR 数字滤波器可以看成是IIR 数字滤波器的a k 全都为0时的一个特例。 IIR 数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、直接Ⅲ型、级联型和并联型。 FIR 数字滤波器的基本结构分为横截型（又称直接型或卷积型）、级联型、线性相位型及频率采样型等。本实验对线性相位型及频率采样型不做讨论，见实验10、12。 另外，滤波器的一种新型结构——格型结构也逐步投入应用，有全零点FIR 系统格型结构、全极点IIR 系统格型结构以及全零极点IIR 系统格型结构。 2、IIR 数字滤波器的基本结构与实现 （1）直接型与级联型、并联型的转换 例6-1 已知一个系统的传递函数为 -1-2-3 -1-2-3 8-4z +11z -2z H(z)=1-1.25z +0.75z -0.125z 将其从直接型（其信号流图如图6-1所示）转换为级联型和并联型。

数字信号处理实验

实验六: 用FFT对信号作频谱分析 一、实验目的 1．了解双音多频信号的产生、检测、包括对双音多频信号进行DFT时的参数选择等。 2．初步了解数字信号处理在是集中的使用方法和重要性。 3．掌握matlab的开发环境。 二、实验原理与方法 1、引言 双音多频（Dual Tone Multi Frequency, DTMF）信号是音频电话中的拨号信号，由美国AT&T贝尔公司实验室研制，并用于电话网络中。这种信号制式具有很高的拨号速度，且容易自动监测识别，很快就代替了原有的用脉冲计数方式的拨号制式。这种双音多频信号制式不仅用在电话网络中，还可以用于传输十进制数据的其它通信系统中，用于电子邮件和银行系统中。这些系统中用户可以用电话发送DTMF信号选择语音菜单进行操作。DTMF信号系统是一个典型的小型信号处理系统，它要用数字方法产生模拟信号并进行传输，其中还用到了D/A变换器；在接收端用A/D变换器将其转换成数字信号，并进行数字信号处理与识别。为了系统的检测速度并降低成本，还开发一种特殊的DFT算法，称为戈泽尔(Goertzel)算法，这种算法既可以用硬件（专用芯片）实现，也可以用软件实现。下面首先介绍双音多频信号的产生方法和检测方法，包括戈泽尔算法，最后进行模拟实验。下面先介绍电话中的DTMF信号的组成。在电话中，数字0~9的中每一个都用两个不同的单音频传输，所用的8个频率分成高频带和低频带两组，低频带有四个频率：679Hz,770Hz,852Hz和941Hz；高频带也有四个频率：1209Hz,1336Hz,1477Hz和1633Hz.。每一个数字均由高、低频带中各一个频率构成，例如1用697Hz和1209Hz两个频率，信号用表示。这样8个频率形成16种不同的双频信号。具体号码以及符号对应的频率如表10.6.1所示。表中最后一列在电话中暂时未用。DTMF信号在电话中有两种作用，一个是用拨号信号去控制交换机接通被叫的用户电话机，另一个作用是控制电话机的各种动作，如播

数字信号处理实验作业

实验5 抽样定理 一、实验目的： 1、了解用MA TLAB 语言进行时域、频域抽样及信号重建的方法。 2、进一步加深对时域、频域抽样定理的基本原理的理解。 3、观察信号抽样与恢复的图形，掌握采样频率的确定方法和插公式的编程方法。 二、实验原理： 1、时域抽样与信号的重建 （1）对连续信号进行采样 例5-1 已知一个连续时间信号sin sin(),1Hz 3 ππ=0001f(t)=(2f t)+6f t f ，取最高有限带宽频率f m =5f 0，分别显示原连续时间信号波形和F s >2f m 、F s =2f m 、F s <2f m 三情况下抽样信号的波形。 程序清单如下： %分别取Fs=fm ，Fs=2fm ，Fs=3fm 来研究问题 dt=0.1; f0=1; T0=1/f0; m=5*f0; Tm=1/fm; t=-2:dt:2; f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); subplot(4,1,1); plot(t,f); axis([min(t),max(t),1.1*min(f),1.1*max(f)]); title('原连续信号和抽样信号'); for i=1:3; fs=i*fm;Ts=1/fs; n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); subplot(4,1,i+1);stem(n,f,'filled'); axis([min(n),max(n),1.1*min(f),1.1*max(f)]); end 程序运行结果如图5-1所示：

原连续信号和抽样信号 图5-1 （2）连续信号和抽样信号的频谱 由理论分析可知，信号的频谱图可以很直观地反映出抽样信号能否恢复原模拟信号。因此，我们对上述三种情况下的时域信号求幅度谱，来进一步分析和验证时域抽样定理。 例5-2编程求解例5-1中连续信号及其三种抽样频率（F s>2f m、F s=2f m、F s<2f m）下的抽样信号的幅度谱。 程序清单如下： dt=0.1;f0=1;T0=1/f0;fm=5*f0;Tm=1/fm; t=-2:dt:2;N=length(t); f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); wm=2*pi*fm;k=0:N-1;w1=k*wm/N; F1=f*exp(-j*t'*w1)*dt;subplot(4,1,1);plot(w1/(2*pi),abs(F1)); axis([0,max(4*fm),1.1*min(abs(F1)),1.1*max(abs(F1))]); for i=1:3; if i<=2 c=0;else c=1;end fs=(i+c)*fm;Ts=1/fs; n=-2:Ts:2;N=length(n); f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); wm=2*pi*fs;k=0:N-1; w=k*wm/N;F=f*exp(-j*n'*w)*Ts; subplot(4,1,i+1);plot(w/(2*pi),abs(F)); axis([0,max(4*fm),1.1*min(abs(F)),1.1*max(abs(F))]); end 程序运行结果如图5-2所示。 由图可见，当满足F s≥2f m条件时，抽样信号的频谱没有混叠现象；当不满足F s≥2f m 条件时，抽样信号的频谱发生了混叠，即图5-2的第二行F s＜2f m的频谱图，，在f m=5f0的围，频谱出现了镜像对称的部分。

现代数字信号处理及应用仿真题答案

仿真作业 姓名：李亮 学号：S130101083

4.17程序 clc; clear; for i=1:500 sigma_v1=0.27; b(1)=-0.8458; b(2)=0.9458; a(1)=-(b(1)+b(2)); a(2)=b(1)*b(2); datlen=500; rand('state',sum(100*clock)); s=sqrt(sigma_v1)*randn(datlen,1); x=filter(1,[1,a],s); %% sigma_v2=0.1; u=x+sqrt(sigma_v2)*randn(datlen,1); d=filter(1,[1,-b(1)],s); %% w0=[1;0]; w=w0; M=length(w0); N=length(u); mu=0.005; for n=M:N ui=u(n:-1:n-M+1); y(n)=w'*ui; e(n)=d(n)-y(n); w=w+mu.*conj(e(n)).*ui; w1(n)=w(1); w2(n)=w(2); ee(:,i)=mean(e.^2,2); end end ep=mean(ee'); plot（ep）; xlabel('迭代次数');ylabel('MSE');title('学习曲线'); plot(w1); hold; plot(w2); 仿真结果：

步长0.015仿真结果 0.10.20.30.4 0.50.60.7迭代次数 M S E 学习曲线

步长0.025仿真结果

步长0.005仿真结果 4.18 程序 data_len = 512; %样本序列的长度 trials = 100; %随机试验的次数 A=zeros(data_len,2);EA=zeros(data_len,1); B=zeros(data_len,2);EB=zeros(data_len,1); for m = 1: trials a1 = -0.975; a2 = 0.95; sigma_v_2 =0.0731; v = sqrt(sigma_v_2) * randn(data_len, 1, trials);%产生v(n) u0 = [0 0]; num = 1; den = [1 a1 a2]; Zi = filtic(num, den, u0); %滤波器的初始条件 u = filter(num, den, v, Zi); %产生样本序列u(n) %（2）用LMS滤波器来估计w1和w2 mu1 = 0.05; mu2 = 0.005; w1 = zeros(2, data_len);

数字信号处理实验

数字信号处理实验

实验一 自适应滤波器 一、实验目的 1、掌握功率谱估计方法 2、会用matlab 对功率谱进行仿真 二、实验原理 功率谱估计方法有很多种，一般分成两大类，一类是经典谱估计；另一类是现代谱估计。经典谱估计可以分成两种，一种是BT 法，另一种是周期法；BT 法是先估计自相关函数，然后将相关函数进行傅里叶变换得到功率谱函数。相应公式如下所示： ||1 *0 1 ?()()()(11) ??()(12) N m xx n jwn BT xx m r m x n x n m N P r m e --=∞ -=-∞ =+-=-∑ ∑ 周期图法是采用功率谱的另一种定义，但与BT 法是等价的，相应的功率谱估计如下所示： 21 1? ()()01 (13)N jw jwn xx n P e x n e n N N --== ≤≤--∑ 其计算框图如下所示： 观测数据x(n) FFT 取模的平方 1/N ) (jw xx e ∧ 图1.1周期图法计算用功率谱框图

由于观测数据有限，所以周期图法估计分辨率低，估计误差大。针对经典谱估计的缺点，一般有三种改进方法：平均周期图法、窗函数法和修正的周期图平均法。 三、实验要求 信号是正弦波加正态零均值白噪声，信噪比为10dB，信号频率为2kHZ，取样频率为100kHZ。 四、实验程序与实验结果 （1）用周期图法进行谱估计 A、实验程序： %用周期法进行谱估计 clear all; N1=128;%数据长度 N2=256; N3=512; N4=1024; f=2;%正弦波频率，单位为kHZ fs=100;%抽样频率,单位为kHZ n1=0:N1-1; n2=0:N2-1; n3=0:N3-1; n4=0:N4-1; a=sqrt(20);%由信噪比为10dB计算正弦信号的幅度

西安电子科技大学数字信号处理大作业

数字信号处理大作业 班级：021231 学号： 姓名： 指导老师：吕雁

一写出奈奎斯特采样率和和信号稀疏采样的学习报告和体会 1、采样定理 在进行A/D信号的转换过程中，当采样频率fs.max大于信号中最高频 率fmax的2倍时(fs.max>2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍；采样定 理又称奈奎斯特定理。 (1)在时域 频带为F的连续信号 f(t)可用一系列离散的采样值f(t1),f(t1±Δt)，f(t1±2Δt)，...来表示,只要这些采样点的时间间隔Δt≤1/2F，便可根据各 采样值完全恢复原始信号。 (2)在频域 当时间信号函数f(t)的最高频率分量为fmax时,f(t)的值可由一系列 采样间隔小于或等于1/2fo的采样值来确定,即采样点的重复频率fs ≥2fmax。 2、奈奎斯特采样频率 (1)概述 奈奎斯特采样定理：要使连续信号采样后能够不失真还原，采样频率必须 大于信号最高频率的两倍（即奈奎斯特频率）。 奈奎斯特频率（Nyquist frequency）是离散信号系统采样频率的一半，因哈里·奈奎斯特（Harry Nyquist）或奈奎斯特－香农采样定理得名。采样定理指出，只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽，就可 以真实的还原被测信号。反之，会因为频谱混叠而不能真实还原被测信号。 采样定理指出，只要离散系统的奈奎斯特频率高于采样信号的最高频率或 带宽，就可以避免混叠现象。从理论上说，即使奈奎斯特频率恰好大于信号带宽，也足以通过信号的采样重建原信号。但是，重建信号的过程需要以一个低 通滤波器或者带通滤波器将在奈奎斯特频率之上的高频分量全部滤除，同时还 要保证原信号中频率在奈奎斯特频率以下的分量不发生畸变，而这是不可能实 现的。在实际应用中，为了保证抗混叠滤波器的性能，接近奈奎斯特频率的分 量在采样和信号重建的过程中可能会发生畸变。因此信号带宽通常会略小于奈 奎斯特频率，具体的情况要看所使用的滤波器的性能。需要注意的是，奈奎斯 特频率必须严格大于信号包含的最高频率。如果信号中包含的最高频率恰好为

数字信号处理综合设计实验报告

数字信号处理实验八 调制解调系统的实现 一、实验目的: （1）深刻理解滤波器的设计指标及根据指标进行数字滤波器设计的过程（2）了解滤波器在通信系统中的应用 二、实验步骤: 1.通过SYSTEMVIEW软件设计与仿真工具，设计一个FIR数字带通滤波器，预先给定截止频率和在截止频率上的幅度值，通过软件设计完后，确认滤波器的阶数和系统函数，画出该滤波器的频率响应曲线，进行技术指标的验证。 建立一个两载波幅度调制与解调的通信系统，将该滤波器作为两个载波分别解调的关键部件，验证其带通的频率特性的有效性。系统框图如下： 规划整个系统，确定系统的采样频率、观测时间、细化并设计整个系统，仿真调整并不断改进达到正确调制、正确滤波、正确解调的目的。（参考文件

zhan3.svu） （1）检查滤波器的波特图，看是否达到预定要求； （2）检查幅度调制的波形以及相加后的信号的波形与频谱是否正常； （3）检查解调后的的基带信号是否正常，分析波形变形的原因和解决措施；（4）实验中必须体现带通滤波器的物理意义和在实际中的应用价值。 2.熟悉matlab中的仿真系统； 3.将1.中设计的SYSTEMVIEW（如zhan3.svu）系统移植到matlab中的仿真环境中，使其达到相同的效果； 4.或者不用仿真环境，编写程序实现该系统，并验证调制解调前后的信号是否一致。 实验总共提供三个单元的时间（6节课）给学生，由学生自行学习和自行设计与移植 三、系统设计 本系统是基于matlab的simulink仿真软件设计的基带信号调制与解调的系统，利用matlab自带的数字信号仿真模块构成其原理框图并通过设置载波、带通滤波器以及低通滤波器等把基带信号经过载波调制后再经乘法器、带通滤波器和低通滤波器等电路系统能解调出基带信号。 1、实验原理框图

现代数字信号处理习题

1.设()u n 是离散时间平稳随机过程，证明其功率谱()w 0S ≥。 证明：将()u n 通过冲激响应为()h n 的LTI 离散时间系统，设其频率响应()w H 为 ()001,w -w w 0, w -w w H w ???? 输出随机过程()y n 的功率谱为()()()2y S w H w S w = 输出随机过程()y n 的平均功率为()()()00201 1r 022w w y y w w S w dw S w dw π π π+?-?= =?? 当频率宽度w 0???→时，上式可表示为()()()01 r 00y S w w π =?≥ 由于频率0w 是任意的，所以有()w 0 S ≥ 3、已知：状态方程 )()1,()1()1,()(1n n n n x n n F n x ν-Γ+--=观测方程 )()()()(2n n x n C n z ν+= )()]()([111n Q n n E H =νν )()]()([222n Q n n E H =νν 滤波初值 )]0([)|0(0x E x =ξ } )]]0([)0()]][0([)0({[)0(H x E x x E x E P --= 请简述在此已知条件下卡尔曼滤波算法的递推步骤。 解：步骤1 状态一步预测，即 1 *11)|1(?)1,()|(N n n C n x n n F n x ∈--=--∧ ξξ 步骤2 由观测信号z(n)计算新息过程，即 1*11)|(?)()()|(?)()(M n n C n x n C n z n z n z n ∈-=-=--ξξα 步骤3 一步预测误差自相关矩阵 N N H H C n n n Q n n n n F n P n n F n n P *1)1,()1()1,() 1,()1()1,()1,(∈-Γ--Γ+---=- 步骤4 新息过程自相关矩阵M M H C n Q n C n n P n C n A *2)()()1,()()(∈+-= 步骤5 卡尔曼增益M N H C n A n C n n P n K *1)()()1,()(∈-=- 或 )()()()(1 2n Q n C n P n K H -= 步骤6 状态估计 1*1)()()|(?)|(?N n n C n n K n x n x ∈+=-αξξ 步骤7 状态估计自相关矩阵 N N C n n P n C n K I n P *)1,()]()([)(∈--= 或 )()()()]()()[1,()]()([)(2n K n Q n K n C n K I n n P n C n K I n P H H +---= 步骤8 重复步骤1-7，进行递推滤波计算 4、经典谱估计方法：

数字信号处理实验仿真结果

西安电子科技大学 统计与自适应信号处理实验仿真 题目基于LMS算法的自适应滤波器的 matlab仿真 学院电子工程学院 专业电路与系统 学号 学生姓名 授课教师 撰写日期： 2012年 12 月 20 日

在自适应滤波器中，参数可调的数字滤波器一般为FIR 数字滤波器，IIR 数字滤波器或格型滤波器。图1中，()x n 表示时刻n 的输入信号，()y n 表示时刻n 的输出信号，()d n 表示时刻n 的信号或期望响应信号，()e n 表示时刻n 的误差信号。误差信号为期望响应信号()d n 与输出信号()y n 之差，记为()()()e n d n y n =-。自适应滤波器的系统参数受误差信号控制，并根据()e n 的值而自动调整，使之适合下一时刻(1)n +的输入(1)x n +，以使输出信号(1)y n +更加接近期望信号(1)d n +，并使误差信号(1)e n +进一步减小。当均方误差2[()]E e n 达到最小值时，()y n 最佳地逼近()d n ，系统已经适应了外界环境。 2.2 LMS 算法 1）2[()]E e n 与权值W 的关系 LMS 自适应滤波器通过算法，当2[()]E e n 最小时，滤波器已经调节出适合外部环境的滤波器权值W 。 我们可以先推导2[()]E e n 与加权系数W 的关系式。 写成矩阵形式： 1 0()()[()][][][()]N T T i i i y j W x j X j W W X j -====∑ （1） 误差： ()()()()[][()]T e j d j y j d j W X j =-=- （2） 则： 22 2[()][()[][()]] [()]2[()[()]][][[][()][()][]] T T T T E e n E d j W X j E d j E d j X j W E W X j X j W =-=-+ （3） 令2[][()[()]],[][[()][()]],[()](0)T dd P E d j X j R E X j X j E d j ?===代入式（3），则有： 2211 1 [()][()]2[][][][][] (0)(0)2(0) i m i T T N N N dd i m x x i x d i m i E e j E d j P W W R W WW W ???====-+=+-∑∑∑ （4） 可以从上式看出均方误差2[()]E e n 是加权系数W 的二次函数，它是一个中间上凹的超抛物线曲面，是具有唯一最小值的函数，即2[()]E e n 与W 的关系在几何上是一个“碗形”的多维曲面。为了简单，设W 是唯一的，则2[()]E e n 与W 的关系成为一个抛物线。调节加权系数W 使

现代数字信号处理-第七章-7.17 仿真题

仿真题7.17 现有一个在二维平面内运动的目标，它从（60000m,40000m ）处，以 （-172m/s,246m/s ）的速度出发。在400s 的运动过程中，目标运动速率保持为300m/s ，并在56~105s,182~245s,285~314s 和348~379s 期间分别以1g,-1.5g,3g 和-2.5g(g=9.8m/2s )的转弯速率进行机动，其余时间段则进行匀速运动。系统在两个方向的观测噪声标准差为m y x 100==σσ。采用IMM 算法实现对该目标的跟踪，其中的模型集合由具有不同转弯速率的协同转弯模型构成。定义状态向量由目标在各方向的位置和速度分量构成，即 ()()()[] T y x n v n y n v n x n x )()(= 在协同转弯模型中，状态转移矩阵及状态噪声输入矩阵分别为 ()()()??????? ?????????T T T T -T -T T --T =-ωωωωωωωωωωωωωcos 0sin 0)sin(1) cos(10)sin(0cos 0)cos(10)sin(1)1,(n n F ()?????? ????????T T T T =-Γ2/00002/1,22n n 其中，ω为转弯速率，T 为采样周期。模型集合由7个协同转弯模型组成，转弯速率分别为 s s s s s s s /6.5,/74.3,/87.1,/0,/87.1,/74.3,/6.57654321 ====-=-=-=ωωωωωωω。转速0ω对应模型的系统状态噪声标准差为1.8m/2s ，其余模型的系统状态噪声标准差为2.5m/2s 。模型初始概率为{0.03，0.03，0.03，0.92，0.03，0.03，0.03}，转移概率矩阵为 ?????????? ????????????=0.90.1000000.10.80.10000 00.10.80.1000000.10.80.1000000.10.80.1000000.10.80.1000000.10.9π 请给出： （1）目标的真实运动轨迹。

DSP与数字信号处理作业

1、什么是DSP？简述DSPs的特点？简述DSPs与MCU、FPGA、ARM的区别？学习DSP开发需要哪些知识？学习DSP开发需要构建什么开发环境？(15分) 答：（1）DSP是Digital Signal Processing(数字信号处理的理论和方法）的缩写，同时也是Digital Signal Processor（数字信号处理的可编程微处理器）的缩写。通常流过器件的电压、电流信号都是时间上连续的模拟信号，可以通过A/D器件对连续的模拟信号进行采样，转换成时间上离散的脉冲信号，然后对这些脉冲信号量化、编码，转化成由0和1构成的二进制编码，也就是常说的数字信号。DSP能够对这些数字信号进行变换、滤波等处理，还可以进行各种各样复杂的运算，来实现预期的目标。 （2）DSP既然是特别适合于数学信号处理运算的微处理器，那么根据数字信号处理的要求，DSP芯片一般具有下面所述的主要特点：1）程序空间和数据空间分开，CPU可以同时访问指令和数据； 2）在一个指令周期内可以完成一次乘法和一次加法运算； 3）片内具有快速RAM，通常可以通过独立的数据总线在程序空间和数据空间同时访问； 4）具有低开销和无开销循环及跳转的硬件支持； 5）具有快速的中断处理和硬件I/O支持； 6）可以并行执行多个操作； 7）支持流水线操作，使得取址、译码和执行等操作可以重复执行。（3）DSP采用的是哈佛结构，数据空间和存储空间是分开的，通过

独立的数据总线在数据空间和程序空间同时访问。而MCU采用的是冯·诺依曼结构，数据空间和存储空间共用一个存储器空间，通过一组总线（地址总线和数据总线）连接到CPU）。很显然，在运算处理能力上，MCU不如DSP；但是MCU价格便宜，在对性能要求不是很高的情况下，还是很具有优势的。 ARM是Advanced RISC（精简指令集）Machines的缩写是面向低运算市场的RISC微处理器。ARM具有比较强的事务管理功能，适合用来跑跑界面、操作系统等，其优势主要体现在控制方面，像手持设备90%左右的市场份额均被其占有。而DSP的优势是其强大的数据处理能力和较高的运算速度，例如加密/解密、调制/解调等。 FPGA是Field Programmable Gate Array（现场可编程门阵列）的缩写，它是在PAL、GAL、PLD等可编程器件的基础上进一步发展的产物，是专用集成电路中集成度最高的一种。FPGA采用了逻辑单元阵列LCA（Logical Cell Array）的概念，内部包括了可配置逻辑模块CLB、输入/输出模块IOB、内部连线三个部分。用户可以对FPGA内部的逻辑模块和I/O模块进行重置配置，已实现用户自己的逻辑。它还具有静态可重复编程和动态在系统重构的特性，使得硬件的功能可以像软件一样通过编程来修改。使用FPGA来开发数字电路，可以大大缩短设计时间，减少PCB面积，提高系统的可靠性；同时FPGA可以用VHDL或Verilog HDL来编程，灵活性强。由于FPGA能够进行编程、除错、再编程和重复操作，因此可以充分地进行设计开发和验证。当电路有少量改动时，更能显示出FPGA的优势，其现场编程能力可

现代数字信号处理实验报告

现代数字信号处理实验报告 1、估计随机信号的样本自相关序列。先以白噪声()x n 为例。 (a) 产生零均值单位方差高斯白噪声的1000个样点。 (b)用公式： 999 1?()()()1000x n r k x n x n k ==-∑ 估计()x n 的前100个自相关序列值。与真实的自相关序列()()x r k k δ=相比较，讨论你的估计的精确性。 (c) 将样本数据分成10段，每段100个样点，将所有子段的样本自相关的平均值作为()x n 自相关的估值，即： 999 00 1?()(100)(100) , 0,1,...,991000x m n r k x n m x n k m k ===+-+=∑∑ 与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列()()x r k k δ=吗？ (d)再将1000点的白噪声()x n 通过滤波器1 1 ()10.9H z z -= -产生1000点的y (n )，试重复(b)的工作，估计y (n )的前100个自相关序列值，并与真实的自相关序列()y r k 相比较，讨论你的估计的精确性。 仿真结果： (a)

图1.1零均值单位方差高斯白噪声的1000个样本点 分析图1.1：这1000个样本点是均值近似为0，方差为1的高斯白噪声。(b) 图1.2() x n的前100个自相关序列值 分析上图可知：当k=0时取得峰值，且峰值大小比较接近于1，而当k≠0时估计的自相关值在0附近有小幅度的波动，这与真实自相关序列r (k)=δ(k) x 比较接近，k≠0时估计值非常接近0，说明了估计的结果是比较精确的。

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

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

第三章仿真作业3.17 （1）代码 clear; N=32; m=[-N+1:N-1]; 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)); un=signal1+signal2+signal3+noise; uk=fft(un,2*N); sk=(1/N) *abs(uk).^2; r0=ifft(sk); r1=[r0(N+2:2*N),r0(1:N)]; r=xcorr(un,N-1,'biased'); figure subplot(2,2,1) stem(m,real(r1)); xlabel('m'); ylabel('FFT估计r1实部'); subplot(2,2,2) stem(m,imag(r1)); xlabel('m'); ylabel('FFT估计r1虚部'); subplot(2,2,3) stem(m,real(r)); xlabel('m'); ylabel('平均估计r实部'); subplot(2,2,4) stem(m,imag(r)); xlabel('m'); ylabel('平均估计r虚部'); 仿真结果

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

xlabel('空间角度' % %% % %%%%%%%%%%%%RootMUSIC 算法%%%%%%%%%%%%%% % syms z % pz=z.^([0:M-1]'; % pz1=(z^(-1.^([0:M-1]; % fz=z^(M-1*pz1*G*G'*pz; % a=sym2poly(fz; % r=roots(a; % r1=abs(r; % for i=1:2*K % [Y,I(i]=min(abs(r1-1; % r1(I(i=inf; % end % for i=1:2*K % theta_esti1(i=asin(-angle(r(I(i/pi*180/pi; % end % %% % % %%%%%%%%%%%%%%%ESPRIT 算法%%%%%%%%%%%%%%%% % S=VR(:,IX(M:-1:M-K+1; % S1=S(1:M-1,:; % S2=S(2:M,:; % fai=S1\S2; % % [U_fai,V_fai]=eig(fai; % for i=1:K % theta_esti2(i=asin(-angle(V_fai(i,i/pi*180/pi; % end % %% % %%%%%%%%%%%%%%%%%%MVDR 算法%%%%%%%%%%%%%%%% % P=[]; % for n=-pi/2:pi/180:pi/2 % a=exp(-j*[0:M-1]'*pi*sin(n; % P=[P,1/(a'*inv(R*a]; % end % P=abs(P/max(abs(P; % P=10*log10(P; % figure(2; % plot(t,P; % title('MVDR 算法' % ylabel('归一化功率谱/dB' % xlabel('空间角度' % %% % %%%%%%%%%%%%%%%%%%F-SAPES 算法%%%%%%%%%%%%%%%%

现代数字信号处理

现代数字信号处理Advanced Digital Signal Processing 东南大学信息科学与工程学院 杨绿溪

教科书、参考书 ?杨绿溪, 现代数字信号处理, 科学出版社, 2008年12月。?胡广书，数字信号处理----理论、算法与实现，清华大学出版社，1997(或2003)年。 ?皇甫堪等，现代数字信号处理，电子工业出版社，2004年6月。 ?丁玉美等，数字信号处理-----时域离散随机信号处理，西安电子科技大学出版社，2002年12月。 ?金连文，韦岗，现代数字信号处理简明教程，清华大学出版社，2004年1月。 ?何子述等，现代数字信号处理及其应用，清华大学出版社，2009年5月。 ?S.Haykin, Adaptive Filter Theory, Prentice Hall, 2001.

课程基本内容 1.离散时间信号处理基础(本科内容复习) 2.离散随机信号分析基础 –离散时间随机信号基本概念? –基本的正交变换(与信号正交展开、去相关) –基本的参数估计方法 3.线性预测和格型滤波器(语音编码应用)? 4.随机信号的线性建模? 5.功率谱估计(与频率估计、子空间分析)? 6.最优线性滤波: 维纳滤波与卡尔曼滤波? 7.自适应滤波器(线性系统的学习)?

可能选讲或简介的内容 8.多速率数字信号处理和滤波器组 9. 神经智能信息处理；压缩感知等 10. 盲信号处理 11.空时、阵列与MIMO信号处理 12.信号的时频分析

第一章离散时间信号处理基础??本科课程内容复习?? ?数字信号与数字信号处理(DSP)概述 ?滤波器－－简单的数字信号处理系统 ?信号的变换－z变换、DTFT、DFT和FFT ?特殊的序列（和对应的滤波器） –全通序列、最小相位序列、线性相位、半正定序列

数字信号处理实验报告

1． 复指数信号的离散傅里叶变换。其中 n j e n x )9.0()(3/π=，n =[0,10] 用MATLAB 求这一有限时宽的序列的傅里叶变换。 2． 试用Mablab 求其有限长序列)100()8.0()(1≤≤=n n x n 与 )180()6.0()(2≤≤=n n x n 的圆周卷积， （N=20），并画出其结果图。（待定） 3． 试用MATLAB 的residuez 函数，求出 12181533325644162)(234234-+-+++++=z z z z z z z z z X 的部分分式展开和。 4． 试用MATLAB 命令求解以下离散时间系统的单位取样响应。 （1）)1()()2()1(4)(3-+=-+-+n x n x n y n y n y （2） )()2(10)1(6)(2 5n x n y n y n y =-+-+ 5． 已知某系统的单位取样响应为()()()[]10)87(--=n u n u n h n ，试用MATLAB 求当激励信号为)5()()(--=n u n u n x 时，系统的零状态响应。 6．

7． 8．a=[3 4 1]; 9．>> b=[1 1]; 10．>> n=0:30; 11．>> x=impDT(n); 12．??? Undefined function or variable 'impDT'. 13． 14．>> h=filter(b,a,x); 15．>> x=impDT(n); 16．>> h=filter(b,a,x); 17．>> 18．>> stem(n,h,'fill'),grid on 19．>> xlabel('n'),title('系统单位取样响应h(n)') 20．>> a=[2.5 6 10]; 21．>> b=[1]; 22．>> n=0:30; 23．>> x=impDT(n); 24．>> h=filter(b,a,x); 25．>> h=filter(b,a,x); 26．>> h=filter(b,a,x); 27．>> stem(n,h,'fill'),grid on 28．>> xlabel('n'),title('系统单位取样响应h(n)')