1、h1自适应滤波器导引自适应滤波器导引Introduction of Adaptive FilterIntroduction of Adaptive Filterh2内容提要最优滤波维纳滤波器卡尔曼滤波器自适应滤波自适应滤波原理最速下降法最小均方算法自适应滤波器的应用h3学习要求了解 （1）噪声和干扰 （2）学习自适应滤波的意义掌握 （1）维纳滤波的信号模型 （2）最小均方误差准则和正交原理 （3）维纳霍夫方程及其最优解的矩阵形式 （4）卡尔曼滤波的信号模型 （5）新息过程与卡尔曼递推公式 （6）自适应滤波原理h4第一节 概述h5噪声、干扰与信号随机信号或随机过程(random process)
2、是普遍存在的。噪声按功率谱密度划分可以分为白噪声（white noise）和色噪声（color noise），我们把均值为零的白噪声叫纯随机信号（pure random signal）。任何其它随机信号都可看成是纯随机信号与确定性信号并存的混合随机信号，或简称为随机信号。h6白噪声白噪声的定义除了要求均值为零、功率谱密度为常数外，并没有对其应当服从哪种概率分布作出任何假设，常见有高斯白噪声、泊松白噪声、柯西白噪声等。根据中心极限定理，高斯白噪声是许多现实世界过程的一个很好的近似，并且能够生成数学上可以跟踪的模型。这些模型用得如此频繁，以至于加性高斯白噪声成了一个标准的缩写词：AWGN（Addi
3、tive white Gaussian noise）。此外，高斯白噪声有着非常有用的统计学特性，因为高斯变量的独立性与不相关性等价。h7白噪声白噪声的自相关函数为狄拉克函数：由于随机过程的功率谱密度是其自相关函数的傅里叶变换，而函数的傅里叶变换为常数，因此白噪声的功率谱密度是平坦的。1,0()()()()()0,0nnnnRE n t n tR h8高斯白噪声(White Gaussian Noise)如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声，如热噪声；MATLAB中的wgn函数可以生成高斯分布的随机白噪声序列，如 生成一个长度为50的实高斯白
4、噪声列向量，功率为2dBW，y1=wgn(50,1,2)，h9加性高斯白噪声Additive white Gaussian noise xn=sin(2*pi*50*t)+randn(1,N);randn是一种高斯分布的随机数发生函数y=awgn(x,snr)t=0:0.1:10;x=sawtooth(t);%Create sawtooth signal.y=awgn(x,10);%Add white Gaussian noise.plot(t,x,t,y)%Plot both signals.legend(Original signal,Signal with AWGN);1010log20
5、log10ssnnPVSNRPVh10要区别干扰（interference）和噪声(noise)两种事实和两个概念。非目标信号（nonobjective signal）都可叫干扰。干扰可以是确定信号，如国内的50Hz工频干扰。干扰也可以是噪声，纯随机信号（白噪声）加上一个直流成分（确定性信号），就成了最简单的混合随机信号。干扰与噪声的区别干扰与噪声的区别h11通常滤波器设计理论仅仅在频域设计一个特定频率响应的IIR和FIR滤波器，即在处理输入信号的过程中滤波器的参数是固定的；当环境发生变化时，仍然会有噪声通过滤波器。通常这类滤波器并非是适应我们需要的最优滤波器。学习自适应滤波的意义h12自适应
6、系统根据当前自身的状态和环境调整自身的参数以达到预先设定的目标。自适应滤波器的系数是根据输入信号，通过自适应算法自动调整的。h13Noise Cancelationh14 最优滤波估计误差定义为期望响应与滤波器输出之差。对滤波器的要求是使估计误差在某种统计意义下“尽可能小”。最优化滤波也称为维纳滤波。输入输入x(0),x(1),x(2),线性离散线性离散时间滤波器时间滤波器H0,H1,H2,输出输出y(n)+期望响应期望响应s(n)或或d(n)估计误差估计误差e(n)-+h15优化统计准则使某个代价函数或性能指标最小化，其中估计误差的均方值的计算简单，实际中使用最广泛使估计误差均方值最小化的准
7、则称为最小均方误差（minimum mean square error,MMSE）准则滤波器的脉冲响应类型多用FIR型IIR滤波器在计算上更为简单一些FIR滤波器稳定性好。大多数应用中，更倾向于使用FIR滤波器22E()E()()mine ns ns n最优滤波的滤波准则最优滤波的滤波准则h16线性最优滤波器两点约束线性最优滤波器两点约束 要使估计在统计优化准则要使估计在统计优化准则最小均方误差下得到最优化，最小均方误差下得到最优化，对滤波器有两点约束：对滤波器有两点约束：滤波器是线性的滤波器是线性的（一方面是为了使信号通过滤波器后不致发生（一方面是为了使信号通过滤波器后不致发生“畸变畸变”，
8、另一方面是为了方便对滤波器的数学分析）；，另一方面是为了方便对滤波器的数学分析）；滤波器是离散时间的滤波器是离散时间的，这将使滤波器可以利用数字硬件或软件，这将使滤波器可以利用数字硬件或软件实现。实现。由于由于FIRFIR滤波器是滤波器是IIRIIR滤波器的特例，这里以滤波器的特例，这里以IIRIIR滤波器作为讨滤波器作为讨论对象。论对象。h17最优化理论基础正交原理定义代价函数为下列均方误差使代价函数最小的条件是实际上就是指，使代价函数最小化的充分必要条件是估计误差e(n)与输入信号x(0),x(1),x(n)正交，这就是著名的正交原理。2()()J nE e n()()0,0,1,2,E
9、x nk e nkh18维纳滤波与卡尔曼滤波的异同(1)维纳滤波和卡尔曼滤波都是解决线性滤波和预测问题的方法，并且都是以均方误差最小为准则的，在平稳条件下两者的稳态结果是一致的。(2)但是它们解决问题的方法有很大区别。(3)维纳滤波只适用于平稳随机过程，卡尔曼滤波就没有这个限制。h19关键词：维纳滤波h20关键词：卡尔曼滤波h21关键词：自适应h22第一节第一节 维纳滤波器维纳滤波器Wiener filterWiener filterh23简介设计维纳滤波器的过程就是寻求在最小均方误差(minimum mean square error,MMSE)下滤波器的单位脉冲响应或传递函数的表达式，其实
10、质就是解维纳霍夫（Wiener-Hopf）方程。这里只讨论因果可实现滤波器的设计。h24设有一个线性系统，它的单位脉冲响应是设有一个线性系统，它的单位脉冲响应是h(n)h(n)，当输入一个观测，当输入一个观测到的随机信号到的随机信号x(n)x(n)，简称，简称观测值观测值，且该信号包含噪声，且该信号包含噪声w(n)w(n)和有用和有用信号信号s(n)s(n)，即，即则输出则输出y(n)y(n)为为()()()x ns nw n (1)()()()()()my nx nh nh m x n m(2)信号与系统信号与系统h25我们希望滤波输出得到的信号与有用信号尽量接近，因此称y(n)为s(n)的
11、估计值，用 来表示，从而就有了维纳滤波器的系统框图，如图1。这个系统的单位脉冲响应h(n)也称为对于s(n)的一种估计器。)()()(nwnsnx)()(nsny)(nh图图1 1 维纳滤波器的输入输出关系维纳滤波器的输入输出关系()s n维纳滤波器维纳滤波器h26)()()(nsnsne(3)从图从图1 1的系统框图中估计到的信号和我们期望得到的有用信号的系统框图中估计到的信号和我们期望得到的有用信号可能不完全相同，这里用可能不完全相同，这里用e(n)来表示真值和估计值之间的误差来表示真值和估计值之间的误差 e(n)e(n)显然是随机变量，维纳滤波和卡尔曼滤波的误差准则就是显然是随机变量，维
12、纳滤波和卡尔曼滤波的误差准则就是最小均方误差准则最小均方误差准则：维纳滤波器维纳滤波器22()()()E e nEs ns n(4)h27一、因果维纳滤波器一、因果维纳滤波器设设h(n)h(n)是物理可实现的，也即是因果序列：是物理可实现的，也即是因果序列：因此，从式因此，从式(1)(1)、(2)(2)、(3)(3)、(4)(4)推导：推导：()0,0h nn0()()()()my ns nh m x nm(5)220()()()()mE e nEs nh m x nm(6)h28因此要使得均方误差最小，根据因此要使得均方误差最小，根据正交原理正交原理可得：可得：即即用相关函数来表达上式，则得
13、到用相关函数来表达上式，则得到维纳霍夫方程的离散形式：维纳霍夫方程的离散形式：从维纳霍夫方程中解出的从维纳霍夫方程中解出的 h就是最小均方误差下的最佳就是最小均方误差下的最佳(optimal)的的h(n)，即，即hopt(n)，此时，此时均方误差为最小：均方误差为最小：0()()()()00,1,2optmEs nhm x nmx njj(7)0()()()()()0optmE s n x njhm E x nm x njj(8)0()()()0 xsoptxxmRjhm Rjmj(9)维纳维纳-霍夫方程霍夫方程22min00()()()()(0)()()optssoptxsmmE e nEs
14、 nhm x nmRhm Rm(10)h29二、有限冲击响应法求解维纳霍夫方程二、有限冲击响应法求解维纳霍夫方程要得到最小均方意义下的线性最优IIR滤波器hopt(n)，现在转化为如何去求解维纳霍夫方程的问题：设hopt(n)是一个因果序列且可以用有限长（N点长）的序列去逼近它，则维纳霍夫方程发生变化 0)()()(0jmjRmhjRmxxoptxs1,2,1,0)()()(10NjmjRmhjRNmxxoptxs(15)h30于是得到于是得到N N个线性方程个线性方程：写成矩阵形式有：写成矩阵形式有：)0()1()2()1()1()0()1(1)2()1()0()1()1()0()1(1)1
15、()1()1()1()0()0()0(0 xxxxxxxsxxxxxxxsxxxxxxxsRNhNRhNRhNRNjNRNhRhRhRjNRNhRhRhRj)1()1()0()1()1()0()0()2()1()2()0()1()1()1()0(NRRRNhhhRNRNRNRRRNRRRxsxsxsxxxxxxxxxxxxxxxxxx维纳霍夫方程的矩阵形式维纳霍夫方程的矩阵形式自相关矩阵自相关矩阵 互相关向量互相关向量h31这样就得到这样就得到维纳霍夫维纳霍夫方程简化矩阵形式：方程简化矩阵形式：(17)式中式中，Hh(0),h(1)h(N-1)，是待求的单位冲击响应；是待求的单位冲击响应；xx
16、xsR HR维纳霍夫方程的矩阵形式解维纳霍夫方程的矩阵形式解只要只要Rxx是非奇异的，就可以求到是非奇异的，就可以求到H：(18)1xxxsHRR求得求得hopt(n)后，这时的均方误差为最小：后，这时的均方误差为最小：12min0()(0)()()NssoptxsmE e nRhm Rm (19)h32关于维纳滤波器的两个重要结论：（1）计算维纳滤波器最优权系数需要预知以下统计量：输入向量的自相关矩阵Rxx 输入向量与期望响应的互相关向量Rxs（2）维纳滤波器实际上是无约束优化(unconstrained optimization)最优滤波器问题的解h33 用有限长的用有限长的h(n)h(n
17、)来实现维纳滤波时，当已知观测值来实现维纳滤波时，当已知观测值的自相关函数和观测值与信号的互相关函数时就可以按的自相关函数和观测值与信号的互相关函数时就可以按照式（照式（1818）在时域里求解；）在时域里求解；但是当但是当N N比较大时，比较大时，计算量很大计算量很大，并且涉及到，并且涉及到求自求自相关矩阵的逆矩阵相关矩阵的逆矩阵问题时，时域求解已不能很好的解决问题时，时域求解已不能很好的解决问题。问题。新的问题：新的问题：h34若信号 s(n)与噪声w(n)互不相关，即则有0)()(mRmRwssw()()()()()()()ssE x n s nmE s n s nmw n s nmRm(
18、)xsRm()()()()()()sswwEs nw ns nmw nmRmRm()xxRm解决办法之一：利用相关函数的性质解决办法之一：利用相关函数的性质h35则式（则式（1515）和式（）和式（1919）化简为：）化简为：10()()()()0,1,2,1NssoptsswwmRjhmRjmRjmjN(20)10min2)()()0()(NmssoptssmRmhRneE(21)维纳霍夫方程的简化形式维纳霍夫方程的简化形式h36mssmR6.0)(例1已知x(n)=s(n)+w(n)且s(n)与w(n)统计独立，其中s(n)的自相关序列为 w(n)是方差为1的白噪声，试用2阶（N=2）维纳
19、滤波器来估计 s(n)，并求最小均方误差。h37mssmR6.0)()()(mmRww代入式（代入式（2020）得得 )1(2)0(6.06.01)1(6.0)0(210hhjhhj解得：解得：h(0)0.451，h(1)0.165。将上述结果代入式（将上述结果代入式（2121），求得最小均方误差：），求得最小均方误差：45.0)1(6.0)0(1)()()0()(10min2hhmRmhRneEmssss 若要进一步减小误差可以适当增加维纳滤波的阶数，但相若要进一步减小误差可以适当增加维纳滤波的阶数，但相应的计算量也会增加。应的计算量也会增加。解解依题意，已知信号的自相关和噪声的自相关函数分
20、别为：依题意，已知信号的自相关和噪声的自相关函数分别为：10()()()()NssoptsswwmRjhm RjmRjmh38三、预白化法求解维纳霍夫方程三、预白化法求解维纳霍夫方程从上面分析知求解维纳霍夫方程比较复杂。本节用波德（Bode）和香农（Shannon）提出的预白化法求解维纳霍夫方程，得到系统传递函数 H(z)。Claude ShannonHendrik Wade Bodeh39随机信号都可以看成是由一白色噪声w1(n)激励一个物理可实现的系统或模型A(z)的响应，如图2所示。由于x(n)=s(n)+w(n)，图3在图2的基础上给出x(n)的信号模型。把这两个模型合并，最后得到维纳
21、滤波器的信号模型，如图4所示，其中传递函数用B(z)表示。)(ns)(zA1()w n2()w n)(nx)(nx)(zB)(1nw图图3 x(n)3 x(n)的信号模型的信号模型 图图4 4 维纳滤波器的输入信号模型维纳滤波器的输入信号模型信号模型信号模型)(ns)(zA1()w n图图2 s(n)2 s(n)的信号模型的信号模型h401 11 11 111()()()()()()()()()()()()()()()krwwkrl r kwwklwwlkE s n s n mEa k w n ka r w n m ra ka r Rm kra ka kl Rm lRm la k a kl )
22、(mRss()()()()()kf la k a kla lal令令 图图2 2中输出信号中输出信号s(n)s(n)的自相关函数为的自相关函数为1 11 11 1()()()()()()()()sswwwwwwlR mRm l f lRmf mRm a m a m(22)则则 121()()()sswRzA z A z白噪声的自相关函数为 它的z变换就等于 21w)()(2111mmRwwwZ变换 (23)h41同样，对图同样，对图4 4所示的输入输出关系进行所示的输入输出关系进行z z变换得变换得利用卷积性质还可以找到互相关函数之间的关系：利用卷积性质还可以找到互相关函数之间的关系：两边两边
23、z z变换得到：变换得到：121()()()xxwRzB z B z(24)111()()()()()()()()()kw skw sE x n s nmEb k w nks nmb k RkmRmbm()xsRm11()()()xsw sRzRz B z(25)h42 如果已知观测信号的自相关函数，如果已知观测信号的自相关函数，（1 1）首先求它的）首先求它的z z变换，变换，（2 2）然后找到该函数的成对零点、极点，）然后找到该函数的成对零点、极点，（3 3）取其中在单位圆内的那一半零点、极点构成）取其中在单位圆内的那一半零点、极点构成B(z)B(z)，另外在，另外在单位圆外的零、极点构成
24、单位圆外的零、极点构成B(zB(z-1-1)。这样就保证了这样就保证了B(z)B(z)是因果的，并且是最小相位系统。是因果的，并且是最小相位系统。零极点分布与最小相位系统零极点分布与最小相位系统P.S.P.S.最小相位系统：零极点都在最小相位系统：零极点都在z z平面单位圆内的因果系统称为平面单位圆内的因果系统称为最最小相位系统小相位系统。在具有相同幅频特性的同阶系统中，最小相位系统具。在具有相同幅频特性的同阶系统中，最小相位系统具有有最小的延时。最小的延时。h43 从图从图4 4可得可得 由于系统函数由于系统函数B B(z)(z)的零点和极点都在单位圆内，即是一个物理的零点和极点都在单位圆内
25、，即是一个物理可实现的最小相位系统，则可实现的最小相位系统，则 也是一个物理可实现的最小相移网也是一个物理可实现的最小相移网络。络。我们就可以利用式（我们就可以利用式（2626）对）对x(n)进行白化，即把进行白化，即把x(n)当作输入，当作输入，w w1 1(n)(n)当作输出，当作输出，是系统传递函数。是系统传递函数。)()(1)(1zXzBzW(26)(1zB滤波器输入信号滤波器输入信号x(n)的白化的白化)(nx)(zB)(1nw图图4 4 维纳滤波器的输入信号模型维纳滤波器的输入信号模型)(1zBh44)(1zB)()()(zBzGzH(27)()()(nwnsnx)()(nsny)
26、(nh(a)(nx)()(nsny1()B z)(zG)(1nw(b)图图5 5（a a）维纳滤波问题）维纳滤波问题（b b）白化法求解模型）白化法求解模型白化法求解模型白化法求解模型 将图将图1 1重新给出，待求的问题就是最小均方误差下的最佳重新给出，待求的问题就是最小均方误差下的最佳H Hoptopt(z)(z)，如图如图5(a)5(a)所示，为了便于求这个所示，为了便于求这个HHoptopt(z)(z)，将图，将图5(a)5(a)的滤波器分解成两的滤波器分解成两个级联的滤波器：个级联的滤波器：G(z)G(z)和和 h451）对观测信号x(n)的自相关函数 Rxx(m)求z变换得到Rxx(
27、z)；2）利用等式 ，找到最小相位系统B(z)；3）利用均方误差最小原则求解因果的G(z)，由于G(z)的激励源是白噪声，求解变得容易多了；4），即得到维纳霍夫方程的系统函数解。)()()(121zBzBzRwxx)()()(zBzGzH白化法的步骤白化法的步骤建立了上述的模型后，白化法求解维纳霍夫方程步骤如下：建立了上述的模型后，白化法求解维纳霍夫方程步骤如下：)(nx)()(nsny1()B z)(zG)(1nwh46下面我们分析步骤3求解因果的G(z)的过程。(28)01)()()()(mmnwmgnsny按白化法求解模型有按白化法求解模型有:均方误差为：均方误差为：11 1221021
28、11000000()()()()()2()()()()()()()(0)2()()()()()mmmrssw sw wmmrE e nEs ng m w nmE s ns ng m w nmg m w nm g r w nrRg m Rmg mg r Rmr 由于由于 ，代入上式，并且进行配方得：，代入上式，并且进行配方得：)()(2111mmRwww11111222200()1()(0)()()w sssww smmwwRmE e nRg mRm h47均方误差最小也就是上式的中间一项最小，所以均方误差最小也就是上式的中间一项最小，所以0,)()(211mmRmgwswopt(30)注意，这
29、里的注意，这里的 g(m)g(m)是因果的。对该式求是因果的。对该式求z z变换，得到变换，得到112()()w soptwRzGz (31)(1zRsw)(1mRsw表示对表示对求单边求单边z z变换。变换。维纳霍夫方程的系统函数解维纳霍夫方程的系统函数解由式（由式（2525）上式可以表示为上式可以表示为11()()()xsw sRzRz B z(25)111122()()/()()w sxsoptwwRzRzB zGzh48所以维纳所以维纳-霍夫方程的系统函数解表示为：霍夫方程的系统函数解表示为：111122()()()/()()()()()w soptxsoptwwRzGzRzB zHz
30、B zB zB z因果的维纳滤波器的最小均方误差为：因果的维纳滤波器的最小均方误差为：111122222min011()(0)()(0)()()ssw sssw smmwwE e nRRmRRm u m(33)利用帕塞伐尔利用帕塞伐尔(Parseval)(Parseval)定理，上式可用定理，上式可用z z域来表示：域来表示：21min1()()()()2ssoptxscdzE e nRzHz Rzjz(34)*围线积分可以取单位圆。围线积分可以取单位圆。(32)h49【例2】已知图1中 ，且 s(n)、w(n)统计独立，其中 s(n)的自相关序列为 是方差为1的单位白噪声，试设计一个物理可实
31、现的维纳滤波器来估计 ，并求最小均方误差。)()()(nwnsnx)(nwmssmR8.0)()(ns，。，h50步骤步骤1 1：)()()(mRmRmRwwssxx求求z z变换变换1110.36()1(1 0.8)(1 0.8)(1 0.5)(1 0.5)1.6(1 0.8)(1 0.8)0.81.25xxRzzzzzzzzmssmR8.0)()()(mmRww0)(mRsw)()(mRmRssxs解解依题意，已知依题意，已知h51步骤步骤2 2：由于由于)()()(121zBzBzRwxx容易找到最小相位系统和白噪声方差容易找到最小相位系统和白噪声方差zzzzB8.0,8.015.01)
32、(1125.1,8.015.01)(1zzzzB6.121wh52步骤步骤3 3：)(zHopt111211()/()1 0.80.36()1.6(1 0.5)(1 0.8)(1 0.5)xswRzB zzB zzzz对括号里面求反变换，注意括号内的收敛域为对括号里面求反变换，注意括号内的收敛域为28.0 z110.360.6(0.8)()0.6(2)(1)(1 0.8)(1 0.5)nnZu nunzz 取因果部分，也就是第一项，所以取因果部分，也就是第一项，所以118.0116.0)5.01)(8.01(36.0zzz0,)5.0(375.0)(nnhn0.375()1-0.5optHzz
33、h53步骤步骤4 4：最小均方误差为最小均方误差为 11()()()210.45(0.6250.5)2(0.8)(1.25)(0.5)ssoptxsccdzRzHz Rzjzzdzjzzz取单位圆为积分围线，有两个单位圆内的极点：取单位圆为积分围线，有两个单位圆内的极点：0.80.8和和0.50.5，求它们的留数和，所以求它们的留数和，所以0.45(0.625*0.80.5)(0.8 1.25)(0.80.5)0.45(0.625*0.50.5)(0.50.8)(0.5 1.25)0.3752min()E e n2min()E e nh54第第二节节 维纳滤波器的应用维纳滤波器的应用Appli
34、cation of Wiener Filter h55诱发电位中的P300分量是受试者辨认靶刺激时，在头顶中缝区域内诱发出的300ms左右潜伏期的正向电位，它反映脑功能状态，尤其与认知心理过程关系密切。诱发电位诱发电位(Evoked potentials，EP)h56 自自19651965年年SuttonSutton首次报告首次报告P300P300电位以来，引起了国内电位以来，引起了国内外学者的极大关注，外学者的极大关注，P300P300作为反映认知功能的内源性成分已作为反映认知功能的内源性成分已广泛用于神经科学、神经医学和认知心理学等方面的研究。广泛用于神经科学、神经医学和认知心理学等方面的
35、研究。h57诱发电位（EP）提取 作为一种EP，P300经常被淹没在很强的背景噪声中，如眼电、肌电、自发脑电等，所以它的提取有一定难度。已有很多方法用于提高EP的信噪比,它们可进一步用于从背景噪声中提取P300。基本的方法有两种：平均和滤波。前者包括相干平均、加权平均等，它们的主要特点是信噪比的提高与刺激次数的平方根成正比，但过多的刺激次数易导致神经疲劳，损失EP的动态变化及细节信息。h58维纳滤波而滤波技术是抑制噪声，提高信噪比，降低刺激次数的有效方法。将维纳滤波应用于EP的提取：一定数量刺激平均后的EP作为已知信号s，不包括在求平均的刺激中的一次EP作为输入信号x，设计维纳滤波器输出 逼近
36、已知信号s：ys1xxxsHRRh59与ICA结合考虑到维纳滤波器本身不能够有效处理EP与噪声在频域上的相叠成分，独立分量分析(Independent component analysis,ICA)是近年来发展起来的一种多道信号处理方法，它将多道观察信号按照统计独立的原则通过优化算法分解为若干独立成分，并可去除伪迹噪声，实现信号的增强和分析。研究表明，由头皮记录到的EP基本满足线性混合ICA模型，主要成分间相互独立。h60EP仿真信号提取结果比较h61维纳滤波和FastICA结合提取P300整个提取过程的算法实现如下:(1)由一次刺激的P300及不包括此次刺激的少量其他P300平均值构造维纳滤
37、波器，对所有参与提取的少量P300重复上述过程，并将滤波后的结果进行叠加平均，得到新的P300。(2)利用FastICA进行伪迹噪声与原始P300的分离。参考文献：李晓欧,张笑微,冯焕清.基于维纳滤波和快速独立分量分析的P300提取方法J.数据采集与处理,2004,19(3):317-322h62P300提取结果比较h63频谱比较Fz，Cz和和Pz导联导联150次相干平均与次相干平均与17次刺激提取结果的频谱比较次刺激提取结果的频谱比较结论：将维纳滤波与结论：将维纳滤波与FastICAFastICA结合，取得了较好的结合，取得了较好的EPEP提取效果，同时提取效果，同时由较少的刺激次数就可获得
38、易于辨识和测量的由较少的刺激次数就可获得易于辨识和测量的P300P300，避免了过多刺避免了过多刺激导致的神经疲劳及波形失真激导致的神经疲劳及波形失真，是诱发电位少次提取的一种有效方法。，是诱发电位少次提取的一种有效方法。h64第三节第三节 卡尔曼滤波卡尔曼滤波Kalman FilteringKalman Filteringh65维纳滤波器的局限计算维纳滤波器最优权系数需要预知以下统计量：（1）输入向量的自相关矩阵Rxx；（2）输入向量与期望响应的互相关向量Rxs；若期望相应s(n)未知，则无法进行线性最优滤波。实时信号处理时，每到一个新的样本，维纳滤波器需使用以往所有数据重新计算全部自相关和
39、互相关项，计算量很大。h66背景介绍：背景介绍：Kalman，匈牙利数学家。卡尔曼滤波器源于他的博士论卡尔曼滤波器源于他的博士论 文和文和19601960年发表的论文年发表的论文A New A New Approach to Linear Filtering Approach to Linear Filtering and Prediction Problems and Prediction Problems（线（线 性滤波与预测问题的新方法）。性滤波与预测问题的新方法）。h67卡尔曼滤波器是维纳滤波理论的发展，具有如下特点：卡尔曼滤波器是维纳滤波理论的发展，具有如下特点：用状态空间概念描述用
40、状态空间概念描述：卡尔曼提出的递推最优估计理论，采卡尔曼提出的递推最优估计理论，采用状态空间信号模型；用状态空间信号模型；递推计算维纳解递推计算维纳解：在算法采用递推形式，卡尔曼滤波能处理在算法采用递推形式，卡尔曼滤波能处理多维和非平稳的随机过程。多维和非平稳的随机过程。卡尔曼滤波理论的提出，克服了卡尔曼滤波理论的提出，克服了维纳维纳滤波理论的局限性滤波理论的局限性，使其在工程上得到了广泛的应用，尤其在控制、制导、导使其在工程上得到了广泛的应用，尤其在控制、制导、导航、通讯等现代工程方面。航、通讯等现代工程方面。卡尔曼滤波卡尔曼滤波h68前面分析了期望响应s(n)存在的情况下的线性最优滤波器，
41、得到了Wiener滤波器。一个自然的问题是：若期望响应未知，又如何进行线性最优滤波呢？本节将从维纳滤波的信号模型引入到卡尔曼滤波的状态空间模型，来回答这个问题。要给出卡尔曼滤波的信号模型，先来讨论过程方程和观测方程。回顾回顾h69)(ns)(zA1()w n2()w n)(nx图图11 11 维纳滤波的信号模型维纳滤波的信号模型图11是维纳滤波的模型，期望响应s(n)可以认为是由白噪声w1(n)激励一个线性系统A(z)的响应，即s(n)可以由以下差分方程表示则A(z)可表示为1()()()w z A zs z1111()()()()pqijijs na s niw nb w nj11111()
42、1ppqqa za zA zb zb zARMAARMA过程过程bj=0ARAR过程过程 Wold分解定理h70一、一、状态方程和观测方程方程和观测方程在卡尔曼滤波中，为了实现递推，假设响应和激励的时域关系可以用一阶AR过程逼近：在卡尔曼滤波中，期望响应s(n)被称为是状态变量(state variable)，在k时刻的状态用S(k)表示，在k-1时刻的状态用S(k-1)表示。)1()1()(1nwnasns(52)h71 激励信号激励信号w1(n)可用向量表示为可用向量表示为WW1 1(k)(k)，称为，称为过程噪声向量过程噪声向量，激励，激励和响应之间的关系用和响应之间的关系用状态转移矩阵
43、状态转移矩阵A(k)A(k)来表示，描述动态系统从时间来表示，描述动态系统从时间n-1n-1的状态到时间的状态到时间n n的状态之间的转移，它应该是已知的。的状态之间的转移，它应该是已知的。有了这些假设后，我们给出状态方程：有了这些假设后，我们给出状态方程：W1S(k)A(k)S(k1)(k1)(53)上式的含义就是在上式的含义就是在k k时刻的状态时刻的状态S(K)S(K)，可以由它的前一个时刻的状，可以由它的前一个时刻的状态态S(k-1)S(k-1)来求得，即认为来求得，即认为k-1k-1时刻以前的各状态都已记忆在时刻以前的各状态都已记忆在S(k-1)S(k-1)状态状态中了中了 状态方程
44、（状态方程（state equation）h722()()()x ns nw n 还是从维纳滤波的观测信号模型入手（图还是从维纳滤波的观测信号模型入手（图1111的右图），观测数据的右图），观测数据和信号的关系为：和信号的关系为：假设假设WW2 2(k)(k)是均值为零的白噪声的是均值为零的白噪声的观测噪声向量观测噪声向量，则观测向量，则观测向量X(k)X(k)与状态向量与状态向量S(k)S(k)之间可表示为之间可表示为上式所表示的卡尔曼滤波一维信号模型和维纳滤波的信号模型是一致上式所表示的卡尔曼滤波一维信号模型和维纳滤波的信号模型是一致的。的。2WX(k)S(k)(k)观测方程（观测方程（m
45、easurement equation）h73把式把式(54)(54)推广就得到更普遍的多维观测方程推广就得到更普遍的多维观测方程2WX(k)C(k)S(k)(k)(55)上式中的上式中的C(k)C(k)称为称为量测矩阵量测矩阵，由于系统状态，由于系统状态S(k)S(k)并不一定已知，并不一定已知，该矩阵的意义就是：描述系统状态经过其作用，变成可观测的，要该矩阵的意义就是：描述系统状态经过其作用，变成可观测的，要求它是已知的。求它是已知的。假设假设X(k)X(k)是是mm 1 1的向量，的向量，S(k)S(k)是是n n 1 1的向量，的向量，则则C(k)C(k)就是就是mm n n的矩阵，的
46、矩阵，WW2 2(k)(k)是是mm 1 1的向量。的向量。h74二、信号模型二、信号模型有了状态方程 和观测方程 后，我们就能给出卡尔曼滤波的信号模型，如图12所示：W1S(k)A(k)S(k1)(k1)2WX(k)C(k)S(k)(k)S(k)C(k)1)A(k 1z2W(k)W1(k)X(k)1)S(k 图图12 卡尔曼滤波的信号模型卡尔曼滤波的信号模型)(ns)(zA1()w n)(nw)(nx维纳滤波的信号模型维纳滤波的信号模型h752WX(k)S(k)(k)10.36()0.81.25(1 0.8)(1 0.8)ssRzzzz，)()(mmRwwA(k)C(k)【例3】设卡尔曼滤波
47、中观测方程为，已知信号的自相关函数的z变换为噪声的自相关函数为：信号和噪声统计独立。求卡尔曼滤波信号模型中的和。h76)()(8.0)1(1kwksks8.0A(k)2WX(k)S(k)(k)=1C(k)变换到时域得：因此 又因为所以112110.36()()()(1 0.8)(1 0.8)sswz zRzA z A zzz解根据等式：可以求得：)()(8.01)(111zWzSzzzAh77卡尔曼滤波建立好了卡尔曼滤波的信号模型以及状态方程、观测方程后，要解决的问题就是要寻找在最小均方误差下信号S(k)的估计值。其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻的估计值和现时刻的观测值
48、来更新对状态变量的估计，求出现在时刻的估计值，从而寻求一套递推估计的算法。它适合于实时处理和计算机运算。h78一、卡尔曼滤波的一步递推法模型一、卡尔曼滤波的一步递推法模型把状态方程和观测方程重新给出：W1S(k)A(k)S(k1)(k1)(56)2WX(k)C(k)S(k)(k)(57)1)(kS(k)S 上式中状态转移矩阵上式中状态转移矩阵A(k)A(k)和量测矩阵和量测矩阵C(k)C(k)是已知的，是已知的，X(k)X(k)是是观测到的数据，也是已知的，假设信号的上一个估计值观测到的数据，也是已知的，假设信号的上一个估计值 已知，现在的问题就是如何来求当前时刻的估计值已知，现在的问题就是如
49、何来求当前时刻的估计值 。h79(k)S(k)X 假设暂不考虑假设暂不考虑w w1 1(k)(k)与与w(k)w(k)，此时得到的，此时得到的S(k)S(k)的估计和的估计和 X(k)X(k)的估计的估计 分别用分别用 和和 表示，得：表示，得：1)(kSA(k)(k)S (58)1)(kSC(k)A(k)(k)SC(k)(k)X(59)(k)X 必然，观测值必然，观测值X(k)X(k)和估计值和估计值 之间有误差，它们之之间有误差，它们之间的差称为间的差称为新息（新息（innovationinnovation）：X(k)X(k)X(k)(60)新息（新息（innovation）一步预测一步预
50、测h80 显然，新息的产生是由于我们前面忽略了显然，新息的产生是由于我们前面忽略了WW1 1(k)(k)与与WW2 2(k)(k)所引起的，也就是说新息里面包含了所引起的，也就是说新息里面包含了w w1 1(k)(k)与与w(k)w(k)的信息成分。的信息成分。因而我们用新息因而我们用新息 乘以一个修正矩阵乘以一个修正矩阵H(k)H(k)，用它来代，用它来代替式（替式（5656）的）的 WW1 1(k)(k)来对来对S(k)S(k)进行估计：进行估计：由于由于H(k)H(k)是作用于新息的是作用于新息的“增益增益”，因此被称为，因此被称为卡尔曼增益卡尔曼增益。(k)X(61)新息与状态方程新息
