数字信号处理实验

实验一 自适应滤波器

一、实验目的

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 法是等价的，相应的功率谱估计如下所示：

2

1

1?

()()01

(13)N jw jwn xx n P e x n e n N 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计算正弦信号的幅度

wn1=randn(1,N1);xn1=a*sin(2*pi*f*n1./fs)+wn1; Pxx1=10*log10(abs(fft(xn1).^2)/N1);%周期法求功率谱f1=((0:length(Pxx1)-1))/length(Pxx1);

wn2=randn(1,N2);xn2=a*sin(2*pi*f*n2./fs)+wn2; Pxx2=10*log10(abs(fft(xn2).^2)/N2);

f2=((0:length(Pxx2)-1))/length(Pxx2);

wn3=randn(1,N3);xn3=a*sin(2*pi*f*n3./fs)+wn3; Pxx3=10*log10(abs(fft(xn3).^2)/N3);

f3=((0:length(Pxx3)-1))/length(Pxx3);

wn4=randn(1,N4);xn4=a*sin(2*pi*f*n4./fs)+wn4; Pxx4=10*log10(abs(fft(xn4).^2)/N4);

f4=((0:length(Pxx4)-1))/length(Pxx4);

subplot(2,2,1);

plot(f1,Pxx1);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=128');

subplot(2,2,2);

plot(f2,Pxx2);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=256');

subplot(2,2,3);

plot(f3,Pxx3);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=512');

subplot(2,2,4);

plot(f4,Pxx4);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=1024');

B、实验仿真结果：

（2）采用汉明窗，分段长度L=32，用修正的周期图求平均法进行谱估计

A：实验程序:

clear all;

N=512;%数据长度

Ns=32;%分段长度

f1=2;%正弦波频率，单位为kHZ

fs=100;%抽样频率,单位为kHZ

n=0:N-1;

a=sqrt(20);%由信噪比为10dB计算正弦信号的幅度

wn=randn(1,N);

xn=a*sin(2*pi*f1*n./fs)+wn;

w=hamming(32)';%汉明窗

Pxx1=abs(fft(w.*xn(1:32),Ns).^2)/norm(w)^2;

Pxx2=abs(fft(w.*xn(33:64),Ns).^2)/norm(w)^2;

Pxx3=abs(fft(w.*xn(65:96),Ns).^2)/norm(w)^2;

Pxx4=abs(fft(w.*xn(97:128),Ns).^2)/norm(w)^2;

Pxx5=abs(fft(w.*xn(129:160),Ns).^2)/norm(w)^2;

Pxx6=abs(fft(w.*xn(161:192),Ns).^2)/norm(w)^2;

Pxx7=abs(fft(w.*xn(193:224),Ns).^2)/norm(w)^2;

Pxx8=abs(fft(w.*xn(225:256),Ns).^2)/norm(w)^2;

Pxx9=abs(fft(w.*xn(257:288),Ns).^2)/norm(w)^2;

Pxx10=abs(fft(w.*xn(289:320),Ns).^2)/norm(w)^2;

Pxx11=abs(fft(w.*xn(321:352),Ns).^2)/norm(w)^2;

Pxx12=abs(fft(w.*xn(353:384),Ns).^2)/norm(w)^2;

Pxx13=abs(fft(w.*xn(385:416),Ns).^2)/norm(w)^2;

Pxx14=abs(fft(w.*xn(417:448),Ns).^2)/norm(w)^2;

Pxx15=abs(fft(w.*xn(449:480),Ns).^2)/norm(w)^2;

Pxx16=abs(fft(w.*xn(481:512),Ns).^2)/norm(w)^2;

Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8+Pxx9 +Pxx10+Pxx11+Pxx12+Pxx13+Pxx14+Pxx15+Pxx16)/16);

f=((0:length(Pxx)-1))/length(Pxx);

plot(f,Pxx);xlabel('频率');ylabel('功率(dB)');

title('加窗平均周期图法功率谱Pxx，N=512');

grid on;

B：实验仿真结果：

五．参考文献：

[1]丁玉美，阔永红，高新波.数字信号处理-时域离散随机信号处理[M].西安：

西安电子科技大学出版社，2002.

[2]万建伟，王玲.信号处理仿真技术[M].长沙：国防科技大学出版社，2008.

实验二 卡尔曼滤波器的设计

一．实验目的

1.熟悉并掌握卡尔曼滤波、自适应滤波和谱估计的原理。

2.可以仿真符合要求的卡尔曼滤波器、自适应滤波器和各种谱估计方法。

3.掌握卡尔曼滤波器的递推公式和仿真方法。

4.熟悉matlab 的用法。 二．实验原理

卡尔曼滤波是用状态空间法描述系统的，由状态方程和测量方程所组成。卡尔曼滤波用前一个状态的估计值和最近一个观测数据来估计状态变量的当前值，并以状态变量的估计值的形式给出。其状态方程和量测方程如下所示：

1(11)

(12)

k k k k k k k k

x A x w y C x v +=+-=+-

其中，k 表示时间，输入信号k w 是一白噪声，输出信号的观测噪声k v 也是一个白噪声，输入信号到状态变量的支路增益等于1，即B=1；A 表示状态变量之间的增益矩阵，可随时间变化，k A 表示第

k 次迭代的取值，C 表示状态变量与输出信号之间的增益矩阵，可

随时间变化，其信号模型如图1.1所示(k 用1k -代替)。

图1.1 卡尔曼滤波器的信号模型

卡尔曼滤波是采用递推的算法实现的，其基本思想是先不考虑输入信号k w 和观测噪声k v 的影响，得到状态变量和输出信号的估计值，再用输出信号的估计误差加权矫正状态变量的估计值，使状态变量估计误差的均方值最小。其递推公式如下所示：

0.020.02111

0.040.04

1??(y x )

(112)(1)(112)1(112)(I )(112)

k k k k k k k k

k k k k k x e x H e a H P P b P e P e c P H P d --------?=+--?''=+-??

'=+--??

'

=--? 假设初始条件11,,,,,,-∧

-k k k k k k k P x y R Q C A 已知，其中

]var[],[0000x P x E x ==∧

，那么递推流程见图1.2所示。

-k P k

P 式(1-3)k

H k x ∧

k

P

图1.2 卡尔曼滤波递推流程图

三．实验要求

一连续平稳的随机信号x(t),自相关()x r e

τ

τ-=，信号x(t)为加性

噪声所干扰，噪声是白噪声，测量值的离散值y(k)为已知。

Matlab仿真程序如下：

%编卡尔曼滤波递推程序，估计信号x(t)的波形

clear all;

clc;

Ak=exp(-0.02); %各系数由前面确定；

Ck=1; Rk=0.1; p(1)=20; %各初值；

Qk=1-exp(-0.04);

p1(1)=Ak*p(1)*Ak'+Qk; %由p1代表p'；

x(1)=0; %设信号初值为0；

H(1)=p1(1)*Ck'*inv(Ck*p1(1)*Ck'+Rk);

zk=[-3.2,-0.8,-14,-16,-17,-18,-3.3,-2.4,-18,-0.3,-0.4,-0.8,-19,-2.0,-1.2,-11, -14,-0.9,0.8,10,0.2,0.5,-0.5,2.4,-0.5,0.5,-13,0.5,10,-12,0.5,-0.6,-15,-0.7,15 ,0.5,-0.7,-2.0,-19,-17,-11,-14] %zk为测量出来的离散值；

N=length(zk); %要测量的点数；

for k=2:N

p1(k)=Ak*p(k-1)*Ak'+Qk; %未考虑噪声时的均方误差阵；

H(k)=p1(k)*Ck'*inv(Ck*p1(k)*Ck'+Rk); %增益方程；

I=eye(size(H(k))); %产生和H（k）维数相同的单位矩阵；

p(k)=(I-H(k)*Ck)*p1(k); %滤波的均方误差阵；

x(k)=Ak*x(k-1)+H(k)*(zk(k)-Ck*Ak*x(k-1)); %递推公式；

end,

x %显示信号x(k)的数据；

m=1:N;

n=m*0.02;

plot(n,zk,'-r*',n,x,'-bo'); %便于比较zk和x(k)在同一窗口输出；xlabel('t/s','Fontsize',16);

ylabel('z(t),x(t)','fontsize',16);

title('卡尔曼滤波递推——x(t)的估计波形与z(t)波形','fontsize',16) legend('观测数据z(t)','信号估计值x(t)',2);

grid;

四．实验结果

五．实验小结

通过卡尔曼滤波估计信号与观测信号比较知，卡尔曼滤波输出的估计信号)(t x与实际观测到的离散值)(t z还是存在一定的误差，卡尔曼滤波是从初始状态

就采用递推方法进行滤波，那么在初值迭代后的一段时间内可能会出现较大的误差，随着迭代进行，各参数逐渐趋于稳定，后面的估计值与观察值的误差就减少了。

六．参考文献

[1]丁玉美，阔永红，高新波.数字信号处理-时域离散随机信号处理[M].西安：

西安电子科技大学出版社，2002.

[2]万建伟，王玲.信号处理仿真技术[M].长沙：国防科技大学出版社，2008.