现代数字信号matlab处理仿真题

3.17

(1)相关函数

仿真代码：

Al=getAk(SNRl);

A2=getAk(SNR2);

A3=getAk(SNR3);%求得信号的幅度；

noise=randn(l,N)+j*randn(l,N);%构建高斯白噪声；

sl=getSk(Al,fl,N);

s2=getSk(A2,f2,N);

S3:getSk(A3,f3,N);;%产生3个复正弦信号

vn=sl+s2+s3+noise;

vk=fft(vn,2*N);%对v(n)补N个零，然后做2N点FFT

swk=((abs(vk)).A2)/N;%计算功率谱估计S(⑹

rO=ifft(swk);%对S(k)做ifft得到

r=[rO(N+2:2*N),rO(l:N)];%根据教程3.1.8式可得

rl=xcorr(vn,N-l,'biased')；%直接计算自相关函数

%%%%%%%%%%%%%%%%%%%%%%%%%取序列实部，虚部％%%%%%

real_r=real(r);

imag_r=imag(r);

real_rl=real(rl);

imag_rl=imag(rl);

subplot(2,2,l);

stem(real_r);

xlabel。基于FFT的自相关函数的实部，)；

ylabelC实部

subplot(2,2,2);

stem(imag_r);

xlabel。基于FFT的自相关函数的虚部')；

ylabelC虚部)

subplot(2,2,3J；

stem(real_rl);

ylabel。实部、

xlabel]估计的自相关函数的实部〕

subplot(2,2,4)；

stem(imag_rl);

xlabel]估计的自相关函数的虚实部〕

ylabel(虚部

functionAK=getAk(SNR)%求得幅度

%%%%%%%%%%%%%%%%%%%由SNR=10log(AA2/2*。A2)%%%%%%%

AK=((10A(SNR/10])*2)A0.5;

functionSk=getSk(Ak,fk,N)

Sk=Ak*exp(j*2*pi*fk

仿真波形:

(2)BT法和周期法估计

仿真程序：

clearall;

clc;

%设定N值可以改变抽样信号的点数，设定M值可以设定加窗的大小，设定N3可以补零，确定实际求fft

的点数。

N=256;%样本观察长度

M=64;%加窗长度

fl=0.15;

f2=0.17;

f3=0.26;%定义3个归一化正弦频率

SNR1=3O;

SNR2=30;

SNR3=27;%定义三个正弦信号信噪比

Al=getAk(SNRl);

A2=getAk(SNR2);

A3=getAk(SNR3);%求得信号的幅度；

noise=randn(l,N)+j*randn(l,N);%构建高斯白噪声；

sl=getSk(Al,fl,N);

s2=getSk(A2,f2,N);

s3=getSk(A3,f3,N);%产生3个复正弦信号

un=sl+s2+s3+noise;

%%%%%%%%%%%%下面采用周期法估计的频谱％%%%%%%%%%%%%

N2=2*N;

u2=zeros(l,N2);

u2(l:N)=un;%对信号序列进行补零

Uw=求DFT变换

Swl=zeros(l,N2);

forw=l:N2

Swl(w)=((abs(Uw(w)))A2)/N;%计算功率谱

end

rO=zeros(l,N2);

rO=ifft(Swl);

r=zeros(l,N2);

r(l:N)=rO(N+l:N2);

r(N+l:N2)=rO(l:N);

M=64;%加窗处理

forn=l:2*M

rl(n)=r((N2-2*M)/2+n);

end%加窗之后的序列为rl

N3=256;%实际进行fft变换的点数。

r2=zeros(l,N3);

r2(l:2*M)=rl;%将加窗之后的rl进行了补零得到的是r2.

Sw2=fft(r2);

Sw2=loglO(abs(Sw2));%取模取对数

temp=zeros(l,N3/2);

temp=Sw2(l:N3/2);

Sw2(l:N3/2)=Sw2((N3/2+l):N3);

Sw2((N3/2+l):N3]=temp;%将0-2*pi的序列折换成-pi到pi的序歹!]。

d=l/(N3-l);

x=zeros(l,N3);

fori=l:N3

x(i)=-0.5+(i-l)*d;

end

Swl=loglO(abs(Swl));%取模取对数

temp=zeros(l,N2/2);

temp=Swl(l:N2/2);

Swl(l:N2/2)=Swl((N2/2+l):N2);

Swl((N2/2+l):N2)=temp;%将0-2*pi的序列折换成-pi到pi的序列。

1=1/(N2-1);

y=zeros(l,N2);

fork=l:N2

y(k)=-0.5+(k-l)*l;

end

subplotfl,2,1);

plot(x,Sw2);

xlabelf'BT法')；

ylabel。归一化的功率谱/dB]

subplotfl,2,2);

plot(y,Swl);

xlabel。周期法

ylabel]归一化的功率谱/dB]

仿真波形：

B

、p

拨

>

m

g

N

®—

2.5-2-

-0.500.5-0.5o0.5

BT法周期法

由仿真结果可以看出：BT法和周期法在相同的阶数下,周期法比BT法分辩率更高，但是波动比BT法更

大。

(3)Levinson-Durbin迭代法

仿真代码：

clearall;

clc;

%设定N值可以改变抽样信号的点数，设定M值可以设定加窗的大小，设定N3可以补零，确定实际求fft

的点数。

N=256;%样本观察长度

p=16;%设定AR模型参数

fl=0.15;

f2=0.17;

f3=0.26;%定义3个归一化正弦频率

SNR1=3O;

SNR2=30;

SNR3=27;%定义三个正弦信号信噪比

Al=getAk(SNRl);

A2=getAk(SNR2);

A3=getAk(SNR3);%求得信号的幅度；

noise=randn(l,N)+j*randn(l,N);%构建高斯白噪声；

sl=getSk(Al,fl,N);

s2=getSk(A2,f2,N);

s3=getSk(A3,f3,N);%产生3个复正弦信号

un=sl+s2+s3+noise;

r0=xcorr(un,p；biased');%直接计算自相关函数

rp=r0(p+l:2*p+l);

%%%%%%%%Lervinson_Durbin算法%%%%%%%%

a(l,l)=-rp(2)/rp(l);

e(l)=rp(l)-abs(rp(2))A2/rp(l);%计算。1

b=0;

form=2:p

ford=l:m-l

a(m,dj=a(m-l,dj+k(mj*conj(a(m-l,m-djj;

end

e(m)=e(m-l)*(l-abs(k(m))A2);

end

%%%%计算16阶AR功率谱％%%%%%%%%%%%%%%%%%%%%%

NF=1024;%%%FFT点数

Sw=e(p)./(abs(fft([l,a(p,:]],NF)].A2);

%%%%%%此下作用相当于函数fftshift

temp=zeros(l,NF/2);

temp=Sw(l:NF/2);

Sw(l:NF/2)=Sw((NF/2+l):NFJ;

Sw((NF/2+l):NF)=temp;%将0-2*pi的序列折换成-pi到pi的序列。

SAR=zeros(l,NF);

forw=l:NF

SAR(w)=loglO((abs(Sw(w))A2));

end

d=l/(NF-l);

x=zeros(l,NF);

fori=l:NF

x(i)=-0.5+(i-l)*d;

end

plot(x,SAR);

xlabel('w/2ir')；

ylabel('归一化的功率谱/dB]

仿真波形：

3.20

(1)Root-MUSIC

仿真代码：

clearall;

cic;

N=1000;%样本个数

fl=0.5;%信号频率

f2=-0.3;

sl=getSk(fl,N);%产生第一个信号

s2=getSk(f2,N);%产生第二个信号

noise=(randn(l,N)+j*randn(l,N))/sqrt(2);%构建高斯白噪声；

un=sl+s2+noise;%产生带噪声的信号

K=2;%信号源个数

M=8;%自相关矩阵阶数

R=getR(un,N,M);%根据教材3-5-18求自相关矩阵

%%%%%%%root_MUSIC进行频率估计％%%%%%%%%%%%%

[U,E]=svd(R);%矩阵U是特征值向量组成的矩阵，E是对角矩阵，对角矩阵的排列是从大到小

e=diag(E);

G=U(:,3:M);

GG=G*G';

co=zeros(2*M-l,l);%由教材可知是关于变量z的2(M-1)次方程，共有2(M-1)个根

form=l:M

co(m:m+M-l)=co(m:m+M-l)+GG(M:-l:l,m);

end

z=roots(co);%求多项式的跟

erro=abs(abs(z)；);%除掉大于1的增根(方程式在此条件下变形的)

ph=angle(z)/(2*pi);

%挑选出最接近单位圆的N个根

disp(z)

disp(erro)

disp(ph)

functionaw=getaw(w,M)

aw=zeros(M,l)；

forj=l:M

aw(j)=exp(-w*(j-l)*i);

end

functionR=getR(x,N,M)

L=N-M+1;

tempx=zeros(M,l);

R=zeros(M,M);

forn=M:N

forj=l:M

tempx(j)=x(n-(j-l));

end

R=R+tempx*tempx,;

end

R=R/L;

functionSk=getSk(fk,N)

Sk=exp(j*2*pi*fk*(0:N-l)+j*2*pi*rand);%产生0到2n上均分布的随机相位信号

仿真结果：

wl=-0.3004w2=0.4996

(2)MUSIC算法

仿真程序：

clearall;

cic;

N=1000;%样本个数

fl=0.5;%信号频率

f2=-0.3;

sl=getSk(fl,N);%产生第一个信号

s2=getSk(f2,N);%产生第二个信号

noise=(randn(l,N)+j*randn(l,N))/sqrt(2);%构建高斯白噪声;

un=sl+s2+noise;%产生带噪声的信号

K=2;%信号源个数

M=8;%自相关矩阵阶数

R=getR(un,N,M);%根据教材3-5-18求自相关矩阵

NF=2048;%定义扫描点数

%%%%%%%MUSIC进行谱峰收索％%%%%%%%%%%%%

[U,E]=svd(R);%矩阵U是特征值向量组成的矩阵，E是对角矩阵，对角矩阵的排列是从大到小

e=diag(E);

G=U(:,3:M);

d=l/(NF-l);%求画图用的横坐标的间隔。

h=zeros(l,NF);

fori=l:NF

h(i)=-0.5+(i-l)*d;

end

y=zeros(l,NF);

forj=l:NF

w=h(j)*2*pi;

aw=getaw(w,M);

y(j)=loglO(abs(l/((aw')*G*(G')*aw)));

end

plot(h,y);

xlabel('w/2n');

ylabelC归一化的功率谱/dB1);

title('MUSIC谱估计,)

仿真波形：

CQP

/

期

神

旨

g

^

—

S

由仿真结果可知：Root-MUSIC算法与MUSIC算法相比，两者误差都很小。

3.21ESPRIT■算法

仿真代码：

clearall;

cic;

formatlong

N=1000;%样本个数

fl=0.5;%信号频率

f2=-0.3;

sl=getSk(fl,N);%产生第一个信号

s2=getSk(f2,N);%产生第二个信号

noise=(randn(l,N)+j*randn(l,N))/sqrt(2);%构建高斯白噪声；

un=sl+s2+noise;%产生带噪声的信号

K=2;%信号源个数

M=8;%自相关矩阵阶数

[X,R]=corrmtx(un,M);

Rxx=R(l:M,l:M);

Rxy=R(l:M,2:M+l);%计算Rxy

%%%%%%%EPRIT进行频率估计％%%%%%%%%%%%%

[U,E]=svd(Rxx);

e=diag(E);

emin=e(end);%获取最小特征值

%%%构造Z矩阵％%%%%%%%%%%%%%%%%%%%%%%%%%%%

Z=zeros(M,M);

fori=l:M-l

Z(i,i+l)=l;

end

Cxx=Rxx-emin*eye(M);

Cxy=Rxy-emin*Z;

[U,E]=eig(Cxx,Cxy);

z=diag(E);

ph=angle(z)/(2*pi);%求所有特征值的相应归一化特征值

erro=abs(abs(z)-l);%与单位圆之间的距离

disp(ph);

disp(erro);

functionR=getR(x,N,M)

L=N-M+1;

tempx=zeros(M,l)；

R=zeros(M,M);

forn=M:N

forj=l:M

tempx(j)=x(n-(j-l));

end

R=R+tempx*tempx,;

end

R=R/L;

functionSk=getSk(fk,N)

Sk=exp(j*2*pi*fk*(0:N-l)+j*2*pi*rand);%产生0至!J2兀上均分布的随机相位信号

仿真结果：

wl=0.496818689408838w2=-0.285563358651823

4.18

仿真代码：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%程序功能：产生512点的样本函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearall;

cic;

N=512;%样本序列长度

M=100;%独立试验次数

sigma_v_2=0.0731;

vn=sqrt(sigma_v_2)*randn(N,l,M);%产生均值为零，方差为0.0731的加性噪声，产生100

次512X1的噪声用于100次独立试验

u0=[00];

a=[l-0.9750.95];

un=filter(l,a,vn);%构建AR模型，其中vn为输入，un为输出产生工00组独

立信号1x512

%display(un);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%程序功能：用LMS算法来估计wl,w2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ul=0.05;%定义两个步长

u2=0.005;

wl=zeros(2,N,M);%初始化权向量.两个列向量对应要估计的wl,w2100个2xN权

向量

w2=zeros(2,N,M);%分别针对ul=0.05,u2=0.005

form=l:M;%100次独立实验

el(:,:,m)=zeros(N,l);%初始化相关参数

e2(:,:,m)=zeros(N,l)；

dl(:,:,m)=zeros(N,l)；

d2(:,:,m)=zeros(N,l)；

forn=3:N-l%信号向量时刻迭代

wl(:,n+l,m)=wl(:,n,m)+ul*un(n-l:-l:n-2,:,m)*conj(el(n,l,m));

w2(:,n+l,m)=w2(:,n,m)+u2*un(n-l:-l:n-2,:,m)*conj(e2(n,l/m));

,

dl(n+l,l,m)=wl(:,n+l/m)*un(n:-l:n-l,:,m);

d2(n+l,l,iri)=w2(:,n+l,m),*un(n:-l:n-l,:,m);

el(n+l,l,m)=un(n+l,:,m)-dl(n+l,l,m);

e2(n+l,l/rn)=un(n+l,:,m)-d2(n+l,l,m);

end

end

wl_ave=zeros(2,N);

w2_ave=zeros(2,N);

el_ave=zeros(N,l)；

e2_ave=zeros(N,l)；

form=l:M

wl_ave=wl_ave+wl(:,:,m);

w2_ave=w2_ave+w2(m);

el_ave=el_ave+el(:,:,m),A2;

e2_ave=e2_ave4-e2(:,:,m),A2;

end

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

w2_ave=w2_ave/M;

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

e2_ave=e2_ave/M;

t=l:N;

figure(l)

plot(t/wl(l,:,100),t,wl(2,:,100),t,wl_ave(l,:),t,wl_ave(2,:))

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

title('步长=0.05')

figure(2)

plot(t,w2(L：,100),t,w2(2,:,100),t,w2_ave(L：),t,w2_ave(2,:))

xlabelC迭代次数);ylabel(权向量,)

title('步长=0.005')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%程序功能：独立实验100次计算剩余误差和失调参数，画出学习曲线

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

wopt=zeros(2,M);

Jmin=zeros(l,M);

sum_eig=zeros(M,l)；

%%%%通过维纳-霍夫方程计算最小均方误差

%直接计算自相关函数

form=l:M

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

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

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

wopt(:,m)=R\p;%单次最佳权值

[v,d]=eig(R);

Jmin(m)=rm(512)-p'*wopt(:,m);%单次最小均方误差

sum_eig(m)=d(l,l)-d(2,2);%单次特征值求和

end

sJmin=sum(Jmin)/M;%100次平均误差

Jexl=el_ave-sJmin;

Jex2=e2_ave-sJmin;%计算剩余均方误差

sum_eig_100=sum(sum_eig)/M;%100次特征值总和

Jexfinl=ul*sJmin*(sum_eig_100/(2-ul*sum_eig_100));

Jexfin2=u2*sJmin*(sum_eig_100/(2-u2*sum_eig_100));

Ml=Jexfinl/sJmin;%失调参数ml

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

figure⑶

plot(t,el_ave,'*大e2_ave,T');

xlabel('迭代次数，)；

ylabel”匀方误差。；

legendCul=0.05'「u2=0.005。

title。学习曲线)

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

display(Ml);

display(M2);

仿真波形：

迪

迭代次数迭代次数

学习曲线

1

专u1=0.05

u2=0.005

Ml=-0.002061056056484

M2=-2.064886318291408e-04

由仿真结果可知：当步长较大时，收敛速度较快，但失调和剩余均方误差较大，从而使

稳态性能变差；而步长较小时，收敛速度虽然较慢，但失调和剩余均方误差减小，从而可以

改善稳态性能；

问题：当计算失调参数时，采用教材式(4.4.28)时，仿真得到的失调参数，较式(4.4.27)

误差较大。分析原因，式(4.4.28)只有在步长非常小时，误差才小。

%Ml=ul*sum_eig_100/(2-ul*sum_eig_100);%由教材式(4.4.28)

%M2=u2*sumeig100/(2-u2*sumeig100);

5.10仿真结果及图形

⑴M=2时，q=-0.99,cr2=0.93627,求解Yule-Walker方程

r(0)r(l)1cr2

r(l)r(0)£aJ-[0_

得到相关矩阵

27.048746.5783一

&=

'46.578347.0487

计算程序为：

%由yuler-walker方程计算M=2时的自相关矩阵

r2=inv([l,-0.99;-0.99,l])*[v_sigmal;0];

R2=[r2(l),r2(2);r2(2),r2(l)];%M1=2

⑵M=3时，%=一0・99,%=°b=0.93627,求解Yule-Walker方程

>(0)r(Dr(2)~-1-CT2

—

rd)r(0)rd)a10

/(2)r(Dr(0)_〃2_0

得到自相关矩阵为

-47.048746.578346.1125'

居=46.578347.048746.5783

46.112546.578347.0487

计算程序为：

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

R3=[r3(l),r3(2),r3(3);r3(2),r3(l)/r3(2);r3(3),r3(2),r3(l)];%M2=3

⑶计算特征值扩展

%%%%%%%计算特征扩展值％%%%%%%%

eigl=eig(R2);%计算特征值

eig2=eig(R3);

eig_spreadl=max(eigl)/min(eigl);

eig_spread2=max(eig2)/min(eig2);

%display(eig_spreadl);

%display(eig_spread2);

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

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

⑷根据LMS算法均方误差收敛特性，M=2时步长因子应在区间(0,0.0106)之间，M=3时，

步长因子应在区间(0,0.00708)之间，因此题中的步长因子不合理，仿真时不收敛。故在仿

真中，M=2时采用步长因子0.001,M=3时采用步长因子0.0006.

仿真程序：％%%%%%%%%用LMS算法实现对u(n)的线性预测估计％%%%%%%%%%%%

clearall;

cic;

Ml=2;

M2=3;%滤波器权系数参数

L^IOOOO;%产生系统输入白噪声点数,用于迭代次数

v_sigmal=0.93627;%v(n)的方差

N=100;%独立试验次数500次

v=sqrt(v_sigmal)*randn(L,l,N);%产生均值为零，方差为0.93627的加性噪声，产生500次

5000x1的噪声用于500次独立试验

a=[l-0.99];

u=filter(l,a,v);%构建AR模型，其中v为输入，u为输出，产生500组独立信号

5000x1

%由yuler-walker方程计算M=2时的自相关矩阵

r2=inv([l,-0.99;-0.99,l])*[v_sigmal;0];

R2=[r2(l),r2(2);r2(2),r2(l)];%M1=2

%display(R2);

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

R3=[r3(l)/r3(2)/r3(3);r3(2),r3(l),r3(2);r3(3),r3(2)/r3(l)];%M2=3

%display(R3);

%%%%%%%计算特征扩展值％%%%%%%%

eigl=eig(R2);%计算特征值

eig2=eig(R3);

eig_spreadl=max(eigl)/min(eigl);

eig_spread2=max(eig2)/min(eig2);

%display(eig_spreadl);

%display(eig_spread2);

%用1乂5算法实现对u(n)的线性预测估计

mul=0.001;

mu2=0.0005;

wl=zeros(Ml,L,N);%初始化权向量.两个列向量对应要估计的wl,w2500个2xN权向

量

w2=zeros(M2,L,N);%分别针对Ml,M2

form=l:N%独立试验500次针对mul,wl两个权系数

el(:,:,m)=zeros(L,l)；%初始化相关参数

dl(:,:,m)=zeros(L,l);

forn=3:L-l

wl(:/n+lzm)=wl(:,n,m)+mul*u(n-l:-l:n-2/:,m)*conj(el(n,l,m));

dl(n+l,l,rn)=wl(:,n+l,m)'*u(n:-l:n-l,:,m);

el(n+l,l,rn)=u(n+l,:,m)-dl(n+l,l,iTi);

end

end

form=l:N%独立试验500次针对mu2,w2三个权系数

%初始化相关参数

e2(:z:,m)=zeros(L,l);

d2(:,:,m)=zeros(L,l);

forn=4:L-l

w2(:/n+l,m)=w2(:,n,m)+mu2*u(n-l:-l:n-3,:,m)*conj(e2(n,l,m));

d2(n+l,l,rn)=w2(:,n+l,m)'*u(n:-l:n-2,:,m);

e2(n+l,l,iri)=u(n+l,:,m)-d2(n+l,l,iTi);

end

end

wl_ave=zeros(M1,L);

w2_ave=zeros(M2,L);

el_ave=zeros(L,l)；

e2_ave=zeros(L,l)；

form=l:N

wl_ave=wl_ave+wl(:,:,m);

w2_ave=w2_ave+w2(m);

el_ave=el_ave+el(:,:,m),A2;

e2_ave=e2_ave+e2(:,:,m),A2;

end

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

w2_ave=w2_ave/N;

el_ave=el_ave/N;%5100次独立实验的学习曲线均值

e2_ave=e2_ave/N;

t=l:L;

figure(l)

plot(t,wl(l,:,100),t,wl(2,:,100),t,wl_ave(l,:),t,wl_ave(2,:));

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

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

figure(2)

plot(t,el_ave);

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

figure⑶

plot(t,w2(l,:,100),t,w2(2,:,100),t,w2(3,:,100),t,w2_ave(l,:),t/w2_ave(2,:),t,w2_ave(3,:));

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

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

figure(4)

plot(t,e2_ave);

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

仿真波形图:

M1=2,步长0.001的权向量收敛曲线

1

C

5.9

.78

O.

O.

O.

O..6

O..45

O.

O.3

O.

2

V

10002000300040005000600070008000900010000

迭代次数n

M2=3,步长0.0005的权向量收敛曲线

10002000300040005000600070008000900010000

迭代次数n

M1=2,步长0.001的学习曲线

10002000300040005000600070008000900010000

迭代次数n

另附课本中5.2.4burg算法功率谱估计仿真代码及波形

仿真程序：％%%%%%%%%%利用Burg算法进行功率谱估计％%%%%%%%%%%%%%%%

clearall;

cic;

N=256;%观测样本数

M=16;%BURG算法的阶数

sigmal_v=0.001;%定义噪声方差

fl=0.1;

f2=0.25;

f3=0.27;%定义3个归一化正弦频率

Al=l;

A2=l;

A3=0.5;%定义三个正弦信号幅度

v=sigmal_v*randn(l,N)+j*sigmal_v*randn(l,N);%构建高斯白噪声；

sl=getSk(Al,fl,N);

s2=getSk(A2,f2,N);

s3=getSk(A3,f3,N);%产生3个复正弦信号

u=sl+s2+s3+v;%构造观测信号

k=zeros(l,M);

a=zeros(M,M);

ef=zeros(N,M);

eb=zeros(N,M);

ef(l:Nzl)=u;

eb(l:N,l)=u;%burg算法的初始化，ef,eb表示格形滤波器的各阶状态，算法中的重要参数

a(l,l)=(-2)*(eb(l:N-l,l)'*ef(2:N,l))/(ef(2:N,l)'*ef(2:N,l)+eb(l:N-l,l)，*eb(l:N-l,l))；

k(l,l)=a(l,l)；

ef(2:N,2)=ef(2:N,l)+k(l,l)*eb(l:N-l,l);

eb(2:N,2)=eb(l:N-l,l)+conj(k(l,l))*ef(2:N,l);%BURG算法步骤2

%%%burg算法步骤3

fori=2:M%计算M阶滤波器参数，实现burg算法，用递归实现

a(i,i)=(-2)*(eb(i:N-l,i)'*ef(i+l:N,i))/(ef(i+l:N,i)'*ef(i+l:N,i)+eb(i:N-l,i),*eb(i:N-l,i));

k(l,i)=a(i,i);

ef(i+l:N,i+l)=ef(i+l:N,i)+k(l,i)*eb(i:N-l,i);

eb(i+l:N,i+l)=eb(i:N-l,i)+conj(k(l,i))*ef(i+l:N,i);

forj=l:i-l

a(ij)=a(i-l,j)+k(l,i)*conj(a(i-l,i-j));

end

end

NF=1024;%%%FFT点数

Sw=sigmal_v./fftshift((abs(fft([l,a(M,:)],NF)).A2));

forw=l:NF

SAR(w)=loglO((abs(Sw(w))A2));

end

d=l/(NF-l);

x=zeros(l,NF);

fori=l:NF

x(i)=-0.5+(i-l)*d;

end

plot(x,SAR);

xlabel('w/2n');

ylabel('归一化的功率谱/dB);

functionSk=getSk(Ak,fk,N)

Sk=Ak*exp(j*2*pi*fk*(O:N-1));

仿真波形：

6.13仿真结果及图形

仿真代码：

%%%%%%%%%%%%%%基于奇异值分解的MVDR方法进行信号频率估计％%%%%%%%%%%%

clearall;

clc;

N=1000;%样本点数

M=32;%FIR阶数

%M=4;

fl=0.1;

f2=0.25;

f3=0.27;%定义3个归一化正弦频率

SNR1=3O;

SNR2=30;

SNR3=27;%定义三个正弦信号信噪比

Al=getAk(SNRl);

A2=getAk(SNR2);

A3=getAk(SNR3);%求得信号的幅度；

sl=getSk(Al,fl,N);

s2=getSk(A2,f2,N);

s3=getSk(A3,f3,N);%产生3个复正弦信号

noise=randn(l,N)+j*randn(l,N);%构建高斯白噪声1x1000；

u=sl+s2+s3+noise;%输入信号

%%构建矩阵

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

forn=l:N-M+l

A(:,n)=u(M+n-l:-l:n);

end

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

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

%%构建向量

P=1024;

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

w=2*pi*f;

aw=zeros(M,P);

fork=l:P

form=l:M

aw(m,k)=exp(-li*w(k)*(m-l));

end

end

%%计算MVDR谱

Pmvdr=zeros(size(w));

fork=l:P

Pmvdr(k)=l/(aw(:,k)'*invphi*aw(:,k));

end

Pmvdr=abs(Pmvdr/max(abs(Pmvdr)));%归一化

Pmvdr=10*logl0(abs(Pmvdr));

x=-0.5+(0:P-l)*(l/P);

figure(l)

plot(x,Pmvdr)

xlabelCw/2Tf);ylabel('归一化频谱/dB?

title('M=4的MVDR谱,)

%%%%%%%基于教材6.11奇异值分解的MVDR方法

fork=l:P

sw=zeros(l,M);

fori=l:M

sw(i)=(aw(:,k),)*V(:j)/S(Li);

end

Psvd(k)=l/sum((abs(sw)).A2);

end

Psvd=abs(Psvd/max(abs(Psvd)));

Psvd=10*loglO(Psvd);

figure(2)

plot(x,Psvd)

xlabelCw/2TT〕ylabel('归一化频谱/dB?

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

functionAK=getAk(SNR)%求得幅度

%%%%%%%%%%%由SNR=101og(AA2/2*oA2)%%%%%%%

AK二((10八(SNR/10))*2)八0.5;

functionSk=getSk(Ak,fk,N)

Sk=Ak*exp(j*2*pi*fk

滤波器抽头个数为4时，仿真波形

M=4的MVDR谱

0

-5

0

5

rn

/p

期0

紧

？

|5

卬

一

“

0

-15

-2

-2

-30

-3

-4a4032

--0.5-O.0.

M=4的基于SVD的MVDR频谱

0

-5

S5

/P-1

融-2

糜-20

W

I-3

S-35

0

5

-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5

W/2TT

由仿真可知，当阶数M=4时，MVDR的频谱估计分辨率较低，不能分辨出0.25和0.27两点

的频率。

滤波器抽头个数为32时，仿真波形

M=32的MVDR谱

0

5

0

5

m-

p

、

卿0

—

-1

-15

—

5-2

-20

-3

-35

-40

-45

-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5

W/2TT

M=32的基于SVD的MVDR频谱

0

5

CPQ

/0

握

糜

^-1

-15

—

®-2

-20

-3

-45^---------E----------E----------c----------c-----------!---------------------1-----------E----------c---------

-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5

W/2TT

当阶数M=32时，具有较高的分辨率，可以分辨出0.25和0.27两点的频率

6.15仿真结果及图形

仿真代码

%%%用RLS算法实现un的线性预测％%%%%%%%%%%%%%%%%%

clearall;

cic;

N=1000;%样本序列长度

P=500;%独立试验次数

sigma_v_2=0.0995;

vn=sqrt(sigma_v_2)*randn(N,l,P);%产生均值为零，方差为0.0731的加性噪声，产生100次

512X1的噪声用于500次独立试验

u0=[00];

a=[10.99];

un=filter(l,a,vn);%构建AR模型，其中vn为输入，un为输出，产生100组独立信

号1000x1

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

n0=l;

M=2;

b=un(nO+l:N,:,:);

L=length(b);

A=zeros(M,L,P);

form=l:P

uni=[zeros(M-Ll).',un(:,:,m).'].‘；

fork=l:L

A(:,k,m)=unl(M-l+k:-l:k);

end

end

%%RLS算法求最优权向量

delta=0.004;

lambda=0.98;

w=zeros(M,L+l,P);

form=l:P

epsilon=zeros(L,l);

Pl=eye(M)/delta;

fork=l:L

Pln=Pl*A(:,k,m);

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

kn=Pln/denok;

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

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

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

end

MSE(m,:)=abs(epsilon).A2;

end

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

form=l:P

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

MSE_mean=MSE_mean+MSE(m,:);

end

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

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

t=l:1000;

figure⑴

plot(t,w(l,:,100)Xw(2,:,100),t,w_mean(l,:),t,w_mean(2,:))

xlabelC迭代次数n');ylabel(权值。

figure(2)

plot(l:L,MSE_mean)

xlabelC迭代次数n'^ylabelCMSE')

仿真波形

500次独立试验权值

7.14仿真结果及图形

仿真代码

clearall;

cic;

N=2000;%样本序列长度

M=100;%独立试验次数

sigma_v_2=0.0332;

vn=sqrt(sigma_v_2)*randn(l,N,M);%产生均值为零，方差为0.0332的加性噪声，产生100次

1X2000的噪声用于100次独立试验

a=[l-1.61.46-0.6160.1525];

un=filter(l,a,vn);%构建AR模型，其中vn为输入，un为输出，产生100组独立信

号1X2000

%display(un);

%卡曼滤波

N2=2000;

Jmin=0.005;

form=l:100

fori=5:N2

U(:/i,m)=[un(l,i-l,m);un(l,i-2,m);un(l,i-3,m);un(l,i-4,m)];

end

end

W_esti=zeros(4,N2,100)；

form=l:100

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

forl=l:N2-l

P_pre=P_esti;

,

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

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

alpha(l)=un(l,l,m)-(U(:J,m)),*W_esti(:,l,m);

W_esti(:,l+l,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(l,N2,100)；d=zeros(l,N2,100)；MSE=zeros(l,N2);

form=l:100

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

forn=5:N2

d(l,n,m)=W_esti(:,n,m)'*un(l,n-l:-l:n-4,m),;

e(l,n,m)=un(l,n,m)-d(l,n,m);

end

MSE=MSE+e(:,:,m),A2;

end

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

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

t=l:N2;

figure(l)

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

Iegend('wl','w2','w3','w4')

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

figure⑵

plot(t,W_esti(l,:,50),t,W_esti(2,:,50),t,W_esti(3,:,50),t,W_esti(4,:,50))

Iegend('wl','w2','w3','w4');xlabel('迭代次数');ylabel('权值')

figure⑶

semilogy(t,MSE)

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

仿真波形

2