数字信号实验报告Word下载.docx
《数字信号实验报告Word下载.docx》由会员分享,可在线阅读,更多相关《数字信号实验报告Word下载.docx(28页珍藏版)》请在冰豆网上搜索。
代入
去归一,得到实际的模拟滤波器传输函数
之后,通过双线性变换法转换公式
,得到所要设计的IIR滤波器的系统函数
1.2步骤及内容
1)用双线性变换法设计一个巴特沃斯IIR低通数字滤波器。
设计指标参数为:
在通带内频率低于
时,最大衰减小于
在阻带内
频率区间上,最小衰减大于
2)以
为采样间隔,绘制出数字滤波器在频率区间
上的幅频响应特性曲线。
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),'
双线性z变换法获得数字低通滤波器,频率轴取[0,\pi/2]'
\omega'
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维纳滤波器源程序
%%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
K+1
rxx(i)=sum(x(i:
L).*x(1:
L-i+1))/(L-i+1);
r_xx_g=[rxx(K+1:
-1:
2),rxx(1:
K+1)];
rxs(i)=sum(x(i:
L).*s(1:
r_xs_g=[rxs(K+1:
2),rxs(1:
%Whetherornottomeetrou_xx<
r_xx_t=a.^abs([-K:
K]);
r_xx_t(K+1)=r_xx_t(K+1)+1;
r_xs_t=a.^abs([-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<
break;
end
%%FigureTheTheoreticalvalueandReceivedvalueofx(n)correlationfunction
figure
(1);
stem(r_xx_t,'
r'
holdon
stem(r_xx_g,'
o'
'
b'
legend('
theoreticalvalue'
receivedvalue'
TheoreticalvalueandReceivedvalueofx(n)correlationfunction'
%Fromthefollowingfigure,thetheoreticalvalueandreceivedvalueofx(n)'
scorrelationfunctionarealmostthesame.
%TheoriginalsignalcanachievedthroughWienerfilter.
%%FigureThelast100s(n)andx(n)
figure
(2);
stem(s(L-100:
L),'
stem(x(L-100:
s(n)'
x(n)'
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
N
xx(i)=sum(x(i:
Rxx(i,1:
N)=[xx(i:
1),xx(2:
N+1-i)];
xs(i)=sum(x(i:
invRxx=inv(Rxx);
h_g=invRxx*xs'
;
sestimatedvalue
figure(3);
stem(h_t,'
stem(h_g,'
Theoreticalvalue'
Estimatedvalue'
TheoreticalvalueandEstimatedvalueofh(n)'
%Fromthefollowingfigure,thevalueofh(n)canbeestimatedfromWienerFilter.
%Anditissimilartothetheoreticalvalueofh(n).
%%FigureThelast100s(n)andWienerFilters_w(n)
y_l
(1)=s
(1);
%s_w(n)fromWienerfliter
fori=1:
y_l(i+1)=0.724*s(i)+0.238*x(i+1);
figure(4);
stem(y_l(L-100:
WienerFilterS_l(n)'
Thelast100s(n)(red)andWienerFilters_w(n)'
%Fromthefollowingfigure,thesourcesignals(n)canbeestimatedfromWienerFilter,itiscalleds_w(n).
%Anditissimilartotheoriginalsourcesignals(n).
%%FigureThelast100s(n)andFIRfliters_f(n)
9%s_f(n)fromFIRfliter
y_F(i)=sum(h_g(1:
i)'
.*x(i:
1));
fori=10:
L
y_F(i)=sum(h_g'
i-9));
figure(5);
stem(y_F(L-100:
FIRfliterS_R(n)'
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:
%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):
目前为止,我们已经得到了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.4卡尔曼滤波程序
L=200;
Rw=1;
Rv=1;
xe
(1)=0;
Pe
(1)=1e-12;
P
(1)=100;
A(i)=0.95;
B(i)=0;
W(i)=1;
C(i)=1;
V(i)=1;
u(i)=0;
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:
x(k)=0.95*x(k-1)+w(k-1);
y=x+v;
fork=1:
yI(k)=C(k)*x(k);
I=eye(1,1);
t=1:
L;
figure
(1)
plot(t,yI,t,y,'
:
'
时间t'
量测值yI...y'
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);
figure
(2)
plot(t,P);
P'
figure(3)
plot(t,H);
H'
figure(4)
plot(t,x,t,xe,'
纯信号x...滤波结果xe'
figure(5)
plot(t,yI,t,ye,'
量测值yI...滤波结果ye'
4.自适应信号滤波
4.1自适应信号滤波器设计原理
脉冲响应不变法是实现模拟滤波器数字化的一种直观而常用的方法。
他特别适合用于那些对滤波器的时域特性有一定要求的场合。
具体的说他可以保证所设计的IIR滤波器的脉冲响应和相应模拟滤波器的冲激响应应在采样点上完全一致。
脉冲响应不变法也由此得名。
一个模拟滤波器的传递函数可以用有理式表示为
(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自适应滤波器源程序
L=input('
输入样本点个数'
a=input('
输入w的方差'
hI=input('
输入hI'
u=input('
输入步长u'
h0=input('
输入初始值h0'
n=1:
x=randn(1,L);
w=sqrt(a)*x;
s=hI*x;
y=s+w;
plot(n,s,n,y,'
曲线s-曲线y....'
幅度'
迭代次数n'
hh
(1)=0;
h
(1)=0;
%hI+(1-2*u)*(h0-hI);
ss
(1)=0;
sss
(1)=0;
fori=2:
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);
plot(n,h,n,hh,'
曲线h-曲线hh....'
plot(n,s,n,ss,'
曲线s-曲线ss....'
%%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