数字信号实验报告.docx
《数字信号实验报告.docx》由会员分享，可在线阅读，更多相关《数字信号实验报告.docx（28页珍藏版）》请在冰豆网上搜索。
数字信号实验报告
测试信号处理实验报告
班级：
姓名：
专业：
学号：
实验目录
1.有限IIR滤波器的设计1
1.1IIR数字滤波器设计原理1
1.2步骤及内容1
1.3实验结果如下运行结果：
2
1.4IIR数字滤波器源程序2
2．维纳滤波器3
2.1维纳滤波器设计原理3
2.2.步骤及内容3
2.3实验结果4
2.4维纳滤波器源程序6
3卡尔曼滤波9
3.1卡尔曼滤波设计原理9
3.2步骤及内容10
3.3实验结果：
10
3.4卡尔曼滤波程序11
4.自适应信号滤波13
4.1自适应信号滤波器设计原理13
4.2设计步骤与方法14
4.3实验结果15
4.4自适应滤波器源程序17
1.有限IIR滤波器的设计
1.1IIR数字滤波器设计原理
利用双线性变换设计IIR滤波器，首先要设计出满足指标要求的模拟滤波器的传递函数
，然后由
通过双线性变换可得所要设计的IIR滤波器的系统函数
。
如果给定的指标为数字滤波器的指标，则首先要转换成模拟滤波器的技术指标，这里主要是边界频率
的转换，对
指标不作变化。
边界频率的转换关系为
。
接着，按照模拟低通滤波器的技术指标根据相应设计公式求出滤波器的阶数
和
截止频率
；根据阶数
查巴特沃斯归一化低通滤波器参数表，得到归一化传输函数
；最后，将
代入
去归一，得到实际的模拟滤波器传输函数
。
之后，通过双线性变换法转换公式
，得到所要设计的IIR滤波器的系统函数
。
1.2步骤及内容
1）用双线性变换法设计一个巴特沃斯IIR低通数字滤波器。
设计指标参数为：
在通带内频率低于
时，最大衰减小于
；在阻带内
频率区间上，最小衰减大于
。
2）以
为采样间隔，绘制出数字滤波器在频率区间
上的幅频响应特性曲线。
1.3实验结果如下运行结果：
1.4IIR数字滤波器源程序
clear
rp=1;rs=15;
wp=.2*pi;ws=.3*pi;
wap=tan（wp/2）;was=tan（ws/2）;
[n,wn]=buttord（wap,was,rp,rs,'s'）;
[z,p,k]=buttap（n）;
[bp,ap]=zp2tf（z,p,k）;
[bs,as]=lp2lp（bp,ap,wap）;
[bz,az]=bilinear（bs,as,.5）;
[h,f]=freqz（bz,az,256,1）;
plot（f,abs（h）,'r*'）
title（'双线性z变换法获得数字低通滤波器，归一化频率轴'）;
xlabel（'\omega/2\pi'）;
ylabel（'低通滤波器的幅频响应'）;grid;
figure;
[h,f]=freqz（bz,az,256,100）;
ff=2*pi*f/100;
absh=abs（h）;
plot（ff（1:
128）,absh（1:
128）,'r*'）;
title（'双线性z变换法获得数字低通滤波器,频率轴取[0,\pi/2]'）;
xlabel（'\omega'）;
ylabel（'低通滤波器的幅频响应'）;gridon;
2．维纳滤波器
2.1维纳滤波器设计原理
维纳滤波是一种从噪声中提取信号的最佳线性方法，假定一个随机信号下
具有以下形式：
（1）
其中
表示有用信号，
表示噪声。
当其输入一个单位脉冲响应为
的线性系统时，其输出为：
（2）
2.2.步骤及内容
实验中，
其中
，
是零均值方差为
的均匀分布白噪声，
是与
互不相关的均匀分布白噪声，零均值方差为1。
维纳最佳滤波器为
（3）
冲激响应：
则
其中
为
的最佳估值
利用
（4）
其中L为滤波数据长度。
2.3实验结果
1.维纳预测
输入样本点个数4000
输入FIR滤波器阶数10
滤波前的均方误差：
0.997589515425905
FIR滤波器的均方误差：
0.250057137233902
维纳滤波器的均方误差：
0.115678278671965
通过程序运行结果的均方误差比较，可以看出信号x（n）经过维纳过滤后，噪音已基本消除，与
趋于一致。
L固定，则N越大，滤波效果越好；N固定，则L越大，滤波效果越好。
2.4维纳滤波器源程序
clear
%%ProgramFilterinMatlabAuthor:
huanbomaTime:
6/10/2009
L=4000;%Inputthenumberofsignalsamples
N=10;%InputFIRFilterorder
a=0.95;
K=40;
sigma_a2=1-a^2;
a_=[1,-a];
%%Generaterandomsignalx（n）=s（n）+v（n）andtherequirementisrou_xx<0.03,rou_xs<0.01
while
（1）%s（n）isSourcesignal,v（n）isNoise,x（n）isthereceivedsignal
w=sqrt（sigma_a2）*（randn（L,1））;%wisthewhitenoise,itsvarianceissigma_a2
v=randn（L,1）;%visthewhitenoise,itsvarianceis1
s
（1）=w
（1）;%x（n）isgotthroughtheformula:
s（n）=a*s（n-1）+wn（n）
fori=1:
L-1
s（i+1）=a*s（i）+w（i+1）;
x（i）=s（i）+v（i）;
end
x（L）=s（L）+v（L）;
%Correlationfunctionr_xxandr_xs
fori=1:
K+1
rxx（i）=sum（x（i:
L）.*x（1:
L-i+1））/（L-i+1）;
end
r_xx_g=[rxx（K+1:
-1:
2）,rxx（1:
K+1）];
fori=1:
K+1
rxs（i）=sum（x（i:
L）.*s（1:
L-i+1））/（L-i+1）;
end
r_xs_g=[rxs（K+1:
-1:
2）,rxs（1:
K+1）];
%Whetherornottomeetrou_xx<0.03,rou_xs<0.01
r_xx_t=a.^abs（[-K:
K]）;
r_xx_t（K+1）=r_xx_t（K+1）+1;
r_xs_t=a.^abs（[-K:
K]）;
rou_xx=（sum（（r_xx_g-r_xx_t）.^2））/sum（r_xx_t.^2）;
rou_xs=（sum（r_xs_g-r_xs_t）.^2）/sum（r_xs_t.^2）;
ifrou_xx<0.03&rou_xs<0.01
break;
end
end
%%FigureTheTheoreticalvalueandReceivedvalueofx（n）correlationfunction
figure
（1）;
stem（r_xx_t,'r'）;
holdon
stem（r_xx_g,'o','b'）;
legend（'theoreticalvalue','receivedvalue'）;
title（'TheoreticalvalueandReceivedvalueofx（n）correlationfunction'）;
%Fromthefollowingfigure,thetheoreticalvalueandreceivedvalueofx（n）'scorrelationfunctionarealmostthesame.
%TheoriginalsignalcanachievedthroughWienerfilter.
%%FigureThelast100s（n）andx（n）
figure
（2）;
stem（s（L-100:
L）,'r'）;
holdon
stem（x（L-100:
L）,'o','b'）;
legend（'s（n）','x（n）'）;
title（'Thelast100s（n）andx（n）'）;
%Fromthefollowingfigure,therearesomedifferencesbetweenthesourcesignals（n）andreceivedsignalx（n）.
%Buttheyareroughlythesame.Becausex（n）hasnoiseinsametime.
%%FigureTheTheoreticalvalueandEstimatedvalueofh（n）
n=0:
N-1;
h_t=0.238*（0.724）.^n;%h（n）'stheoreticalvalue
fori=1:
N
xx（i）=sum（x（i:
L）.*x（1:
L-i+1））/（L-i+1）;
end
fori=1:
N
Rxx（i,1:
N）=[xx（i:
-1:
1）,xx（2:
N+1-i）];
end
fori=1:
N
xs（i）=sum（x（i:
L）.*s（1:
L-i+1））/（L-i+1）;
end
invRxx=inv（Rxx）;
h_g=invRxx*xs';%h（n）'sestimatedvalue
figure（3）;
stem（h_t,'r'）;
holdon
stem（h_g,'o','b'）;
legend（'Theoreticalvalue','Estimatedvalue'）;
title（'TheoreticalvalueandEstimatedvalueofh（n）'）;
%Fromthefollowingfigure,thevalueofh（n）canbeestimatedfromWienerFilter.
%Anditissimilartothetheoreticalvalueofh（n）.
%%FigureThelast100s（n）andWienerFilters_w（n）
y_l
（1）=s
（1）;%s_w（n）fromWienerfliter
fori=1:
L-1
y_l（i+1）=0.724*s（i）+0.238*x（i+1）;
end
figure（4）;
stem（s（L-100:
L）,'r'）;
holdon
stem（y_l（L-100:
L）,'o','b'）;
legend（'s（n）','WienerFilterS_l（n）'）;
title（'Thelast100s（n）（red）andWienerFilters_w（n）'）;
%Fromthefollowingfigure,thesourcesignals（n）canbeestimatedfromWienerFilter,itiscalleds_w（n）.
%Anditissimilartotheoriginalsourcesignals（n）.
%%FigureThelast100s（n）andFIRfliters_f（n）
fori=1:
9%s_f（n）fromFIRfliter
y_F（i）=sum（h_g（1:
i）'.*x（i:
-1:
1））;
end
fori=10:
L
y_F（i）=sum（h_g'.*x（i:
-1:
i-9））;
end
figure（5）;
stem（s（L-100:
L）,'r'）;
holdon
stem（y_F（L-100:
L）,'o','b'）;
legend（'s（n）','FIRfliterS_R（n）'）;
title（'Thelast100s（n）andFIRfliters_f（n）'）;
%Fromthefollowingfigure,thesourcesignals（n）canbeestimatedfromFIRFilter,itiscalleds_f（n）.
%Anditisalsosimilartotheoriginalsourcesignals（n）.
%%MeanSquareError
e_x=sum（（x-s（1:
L））.^2）/L;%MSEbeforeFilter
e_w=sum（（y_l-s（1:
L））.^2）/L;%MSEfromWienerFilter
e_f=sum（（y_F-s（1:
L））.^2）/L;%MSEfromFIRFilter
3卡尔曼滤波
3.1卡尔曼滤波设计原理
卡尔曼滤波是一种递推估计方法，它也是基于最小均方误差准则，但它不需要全部过去的观测数据，只是根据前一个的估计值和最近一个观测数据来估计信号的当前值。
在稳态情况下，卡尔曼滤波结果和维纳滤波结果相同。
假设已知动态系统的状态方程和量测方程：
首先我们要利用系统的过程模型，来预测下一状态的系统。
假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态：
（1）
式
（1）中，X（k|k-1）是利用上一状态预测的结果，X（k-1|k-1）是上一状态最优的结果，U（k）为现在状态的控制量，如果没有控制量，它可以为0。
到现在为止，结果已经更新，可是，对应于X（k|k-1）还没更新。
我们用P表示：
（2）
式
（2）中，P（k|k-1）是X（k|k-1）对应P（k-1|k-1）是X（k-1|k-1）对应的A’表示A的转置矩阵，Q是系统过程的式子1，2就是卡尔曼滤波器5个公式当中的前两个，也就是对系统的预测。
现在我们有了现在状态的预测结果，然后我们再收集现在状态的测量值。
结合预测值和测量值，我们可以得到现在状态（k）的最优化估算值X（k|k）：
（3）
其中Kg为卡尔曼增益（KalmanGain）：
（4）
目前为止，我们已经得到了k状态下最优的估算值X（k|k）。
但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束，我们还要更新k状态下X（k|k）
（5）
其中I为1的矩阵，对于单模型单测量，I=1。
当系统进入k+1状态时，P（k|k）就是式子
（2）的P（k-1|k-1）。
这样，算法就可以自回归的运算下去。
卡尔曼滤波器的原理基本描述了，式子（1-5）5个基本公式。
根据这5个公式，可以很容易的实现计算机的程序。
3.2步骤及内容
1.编制卡尔曼滤波器的通用程序
2.运行卡尔曼滤波器程序，其中U（k）=0，选择L=100，观察并记录实验结果，分析滤波效果
3.3实验结果：
3.4卡尔曼滤波程序
clear
L=200;
Rw=1;
Rv=1;
xe
（1）=0;
Pe
（1）=1e-12;
P
（1）=100;
fori=1:
L
A（i）=0.95;
B（i）=0;
W（i）=1;
C（i）=1;
V（i）=1;
u（i）=0;
end
r1=randn（1,L）;
r2=randn（1,L）;
w=sqrt（Rw）*r1;
v=sqrt（Rv）*r2;
r=randn
（1）;
x
（1）=xe
（1）+sqrt（Pe
（1））;
fork=2:
L
x（k）=0.95*x（k-1）+w（k-1）;
end
y=x+v;
fork=1:
L
yI（k）=C（k）*x（k）;
end
I=eye（1,1）;
t=1:
L;
figure
（1）
plot（t,yI,t,y,':
'）
xlabel（'时间t'）
title（'量测值yI...y'）;
fork=2:
L
Q（k-1）=W（k-1）*W（k-1）'*Rw;
R（k）=V（k）*V（k）'*Rv;
P1（k）=A（k）*P（k-1）*A（k）'+Q（k-1）;
H（k）=P1（k）*C（k）'*inv（C（k）*P1（k）*C（k）'+R（k））;
P（k）=（I-H（k）*C（k））*P1（k）;
xe（k）=A（k）*xe（k-1）+H（k）*（y（k）-C（k）*A（k）*xe（k-1））+B（k）*u（k-1）;
ye（k）=C（k）*xe（k）;
end
figure
（2）
plot（t,P）;
xlabel（'时间t'）
title（'P'）
figure（3）
plot（t,H）;
xlabel（'时间t'）
title（'H'）
figure（4）
plot（t,x,t,xe,':
'）;
xlabel（'时间t'）
title（'纯信号x...滤波结果xe'）;
figure（5）
plot（t,yI,t,ye,':
'）
xlabel（'时间t'）
title（'量测值yI...滤波结果ye'）;
4.自适应信号滤波
4.1自适应信号滤波器设计原理
脉冲响应不变法是实现模拟滤波器数字化的一种直观而常用的方法。
他特别适合用于那些对滤波器的时域特性有一定要求的场合。
具体的说他可以保证所设计的ＩＩＲ滤波器的脉冲响应和相应模拟滤波器的冲激响应应在采样点上完全一致。
脉冲响应不变法也由此得名。
一个模拟滤波器的传递函数可以用有理式表示为
（4-1）
通过拉式变换，我们就可以得到他的冲激响应
（4-2）
脉冲响应不变法就是要保持脉冲响应不变，即
（4-3）
对上式中的冲激响应序列h（n）作z变换，就可以得到数字滤波器的传递函数
（4-4）
上式（4-3）和式（4-4），我们已在一些讨论采样序列Z变换的章节中反复接触到。
H（Z）就是h（t）理想采样的Z变换。
因此，脉冲响应不变的法也成为标准Z变换法。
一般来说，h（s）的分母多项式阶次总是大于分子多项式的阶次。
假定h（s）没有多重极点（有重极点时要对以下讨论作适当的修正），则式（4-1）可以展成部分分式的形式：
（4-5）
其中
，
均为复数，
是
的极点。
其拉式变换为
（4-6）
由此得到数字滤波器的单位脉冲响应序列
（4-7）
数字滤波器的传递函数则为
（4-8）
上式经过合并简化，成为一般形式的有理分式传递函数
（4-9）
4.2设计步骤与方法
自适应FIR维纳滤波器的LMS算法公式：
因此给定初始值h0,每得到一个样本yn,可以得到一组新的滤波器权系数hn。
如果信号yn为一个M阶的AR模型，则利用LMS算法可以对AR模型参数进行自适应估计。
综上所述，我们可以将脉冲响应不变法设计IIR滤波器的过程归纳为如下几步骤：
a）根据要求，确定数字滤波器和相应模拟滤波器的个临界频率
，
和衰减
，其中
，T为采样周期。
b）安给定的
和
设计模拟滤波器传递函数
c）将
展成部分分式形式，见式（4-5），求
d）将
和
代入式（4-8），并化简得到式（4-9），从而完成数字滤波器传递函数H（z）的设计。
e）分析所设计滤波器的频域特性，检查它是否满足要求。
上述涉及过程的计算比较复杂，对于一个高阶的滤波器的设计（5阶以上），一般只能借助计算机来完成。
本实验要求实验者编制脉冲响应不变数字滤波器设计程序，在计算机上完成上述涉及过程。
4.3实验结果
1．自适应滤波
输入样本点个数200，输入w的方差10，输入hI=1.6，输入步长u=0.06
输入初始值h0=0
可见，u越大，h（n）和s（n）的速度越快。
可见a1（n）、a2（n）随着n的增加趋近于a1、a2.
4.4自适应滤波器源程序
clear
L=input（'输入样本点个数'）;
a=input（'输入w的方差'）;
hI=input（'输入hI'）;
u=input（'输入步长u'）;
h0=input（'输入初始值h0'）;
n=1:
L;
x=randn（1,L）;
w=sqrt（a）*x;
s=hI*x;
y=s+w;
figure
（1）
plot（n,s,n,y,':
'）;
title（'曲线s－曲线y....'）;
ylabel（'幅度'）;
xlabel（'迭代次数n'）;
hh
（1）=0;
h
（1）=0;%hI+（1-2*u）*（h0-hI）;
ss
（1）=0;
sss
（1）=0;
fori=2:
L
hh（i）=hh（i-1）+2*u*（s（i-1）-hh（i-1）*x（i-1））*x（i-1）;
ss（i）=hh（i）*x（i）;
h（i）=hI+（（1-2*u）^i）*（h0-hI）;
sss（i）=h（i）*x（i）;
end
figure
（2）
plot（n,h,n,hh,':
'）;
title（'曲线h－曲线hh....'）;
ylabel（'幅度'）;
xlabel（'迭代次数n'）;
figure（3）
plot（n,s,n,ss,':
'）;
title（'曲线s－曲线ss....'）;
ylabel（'幅度'）;
xlabel（'迭代次数n'）;
%%lmsProgramstart
clearall
%channelsystemorder
sysorder=5;
%Numberofsystempoints
N=2000;
inp=randn（N,1）;
n=randn（N,1）;
[b,a]=butter（2,0.25）;
Gz=tf（b,a,-1）;
%ThisfunctionissubmittedtomakeinverseZ-transform（Matlabcentralfileexchange）
%Thefirstsysorderweightvalue
%h=ldiv（b,a,sysorder）';
%ifyouuseldivthiswillgiveh:
filterweightstobe
h=[0.0976;
0.2873;
0.3360;
0.2210;
0.0964];
y=lsim（Gz,inp）;
%%Addsomenoise
n=n*std（y）/（10*std（n））;
d=y+n;
totallength=size（d,1）;
%Take60pointsfor
升级会员 