基本滤波算法-维纳滤波+卡尔曼滤波+自适应滤波

矩阵在信息处理中的应用课程研讨报告---基本滤波算法的原理及应用目录3.第一部分：维纳滤波原理及应用——蒋芮第二部分：常规Kalman滤波原理及应用——王丹第三部分：Kalman滤波处理语音信号——程谦1.2.目录第四部分：扩展卡尔曼滤波（EKF)与无迹卡尔曼滤波（UKF）原理及应用——宋其岩第五部分：自适应滤波算法原理及应用——李宏伟4.5.目录3.第一部分：维纳滤波原理及应用——蒋芮第二部分：常规Kalman滤波原理及应用——王丹第三部分：Kalman滤波处理语音信号——程谦1.2.目录3.维纳滤波的背景及一般结构维纳-霍夫方程求解维纳滤波的实际应用1.2.维纳滤波理论是由数学家N∙维纳(NorbertWiener,1894~1964)于第二次世大战期间提出的。这一科研成果是这一时期重大科学发现之一，他提出了线性滤波的理论和线性预测的理论，对通信工程理论和应用的发展起了重要的作用。维纳滤波就是为纪念他的重要贡献而命名的。1.维纳滤波的背景1.维纳滤波的一般结构2.维纳霍夫方程求解mmnxmhnsEneE22)()()()(0)]())()()([(20jnxmnxmhnsEmopt)()()(0mjRmhjRxxmoptxs0j)()()()()(jnxmnxEmhjnxnsEopt0j2.维纳霍夫方程求解3.维纳滤波的应用维纳滤波在运动观测中的应用目录3.第一部分：维纳滤波原理及应用——蒋芮第二部分：常规Kalman滤波原理及应用——王丹第三部分：Kalman滤波处理语音信号——程谦1.2.维纳滤波存在的问题：适用于平稳随机信号，时域求解所需数据存储量和计算量较大。解决方法：利用状态方程和递推方法寻找最小均方误差下状态变量的估计值，即：𝑥𝑘=𝑥𝑘−𝑥𝑘min⁡(𝐸𝑥𝑘𝑥𝑘𝑇)𝑥𝑘目录3.引言——温度问题Kalman滤波的递推原理温度问题的Matlab实现1.2.目录3.引言——温度问题Kalman滤波的递推原理温度问题的Matlab实现1.2.考虑一个房间的温度，假设该房间温度内的各处温度是相等的，该房间内没有任何发热或散热的设备，根据经验，我们可以认定该房间的温度满足以下等式：𝑥𝑘=𝐴𝑘𝑥𝑘−1+𝑤𝑘−1，其中𝐴𝑘=1，即下一刻的温度等于现在的温度，然而存在着一定误差𝑤𝑘。另外，在房间内放置一个温度计，但由于制造工艺等方面的因素，该测量值也并不是完全可靠的，即测量值满足以下式子：z𝑘=𝐶𝑘𝑥𝑘+𝑣𝑘，其中𝐶𝑘=1，误差为𝑣𝑘。引言——温度问题因此，针对任意时刻，我们可以得到这个房间的两个温度值：1、根据生活经验得到的估计值（系统预测值）；2、根据温度计得出来的测量值（系统测量值）。下面我们就可以用这两个值结合他们各自的噪声来估算出房间的实际温度。但究竟实际温度是多少？是相信自己的判断还是相信温度计的测量？相信那个比较多一些呢？kalman滤波器就给我们了一种判断。引言——温度问题目录3.引言——温度问题Kalman滤波的递推原理温度问题的Matlab实现1.2.Kalman滤波的系统过程模型和测量模型：𝑥𝑘=𝐴𝑘𝑥𝑘−1+𝑤𝑘−1⁡//系统过程模型z𝑘=𝐶𝑘𝑥𝑘+𝑣𝑘//系统测量模型式中，𝑥𝑘为待估量；z𝑘为通过传感器得到的测量数据；𝐴𝑘为过程矩阵，描述系统由𝑥𝑘−1到𝑥𝑘的状态之间的转移；𝐶𝑘为测量矩阵，描述状态经其作用，变成可量测或可观测的；𝑤𝑘为过程噪声；𝑣𝑘为测量噪声。Kalman滤波的递推原理Kalman滤波的递推原理假设：状态变量的过程矩阵A不随时间发生变化，𝑤𝑘,𝑣𝑘都是零均值白噪声，方差分别是𝑄𝑘和𝑅𝑘⁡，并且初始状态𝑥0、𝑤𝑘、𝑣𝑘彼此都不相关，即：𝑤𝑘~(0，𝑄𝑘)𝑣𝑘~0，𝑅𝑘𝐸𝑤𝑘𝑤𝑇𝑗=𝑄𝑘𝛿𝑘−𝑗𝐸𝑣𝑘𝑣𝑇𝑗=𝑅𝑘𝛿𝑘−𝑗𝐸𝑤𝑘𝑣𝑇𝑗=0Kalman滤波的递推原理基本思想：先不考虑噪声𝑤𝑘、𝑣𝑘的影响，得到𝑥𝑘和z𝑘的估计值𝑥′𝑘和𝑧′𝑘，再用输出信号的估计误差𝑧(𝑘)加权后校正状态变量的估计值𝑥𝑘，使得状态变量的估计误差𝑥𝑘的均方值最小。𝑥′𝑘测量方程𝑧′𝑘𝑧(𝑘)=z𝑘−𝑧′𝑘⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡校正𝑥𝑘𝑥𝑘=𝑥𝑘−𝑥𝑘min𝐸𝑥𝑘𝑥𝑘𝑇𝑥𝑘具体过程：当不考虑噪声𝑤𝑘、𝑣𝑘的影响时，系统过程模型和系统测量模型变为：𝑥′𝑘=𝐴𝑘𝑥𝑘−1𝑧′𝑘=𝐶𝑘𝑥′𝑘=𝐶𝑘𝐴𝑘𝑥𝑘−1则𝑧′𝑘的估计误差为：𝑧(𝑘)=z𝑘−𝑧′𝑘；为了进一步提高估计状态的准确性，我们用该误差𝑧(𝑘)来校正状态变量𝑥𝑘=𝐴𝑘𝑥𝑘−1+𝐻𝑘z𝑘−𝑧′𝑘=𝐴𝑘𝑥𝑘−1+𝐻𝑘z𝑘−𝐶𝑘𝐴𝑘𝑥𝑘−1式中，𝐻𝑘为增益矩阵，本质上就是加权矩阵。Kalman滤波的递推原理具体过程：校正后的状态变量的估计误差和其均方值为：𝑥𝑘=𝑥𝑘−𝑥𝑘𝑃𝑘=𝐸𝑥𝑘𝑥𝑇𝑘=𝐸(𝑥𝑘−𝑥𝑘)(𝑥𝑘−𝑥𝑘)𝑇未经校正的状态变量估计误差的均方值为：𝑃′𝑘=𝐸(𝑥𝑘−𝑥′𝑘)(𝑥𝑘−𝑥′𝑘)𝑇kalman滤波要求状态变量的估计误差的均方值𝑃𝑘最小，因此kalman滤波的关键是得到𝑃𝑘与𝐻𝑘的关系，即通过选择合适的𝐻𝑘使得𝑃𝑘取得最小值。Kalman滤波的递推原理kalman递推步骤为：a、任意选取估计方差的初值𝑃0；foriteration=1:numberb、计算未校正时的估计方差𝑃′𝑘；c、计算增益矩阵𝐻𝑘；d、计算校正过后的状态变量误差的均方差𝑃𝑘−1；end校正状态变量𝑥𝑘。Kalman滤波的递推原理kalman递推公式为：𝑃′𝑘=𝐴𝑘𝑃𝑘−1𝐴𝑇𝑘+𝑄𝑘−1𝐻𝑘=𝑃′𝑘𝐶𝑇(𝑘)(𝐶𝑘𝑃′𝑘𝐶𝑇𝑘+𝑅𝑘)−1𝑃𝑘−1=(𝐼−𝐻𝑘𝐶𝑘)𝑃′𝑘𝑥𝑘=𝐴(𝑘)𝑥𝑘−1+𝐻𝑘(z𝑘−𝐶𝑘𝐴𝑘𝑥𝑘−1)Kalman滤波的递推原理目录3.引言——温度问题Kalman滤波的递推原理温度问题的Matlab实现1.2.现在我们回归到刚开始的温度问题，其中𝐴𝑘=𝐶𝑘=1，假设：A、前一时刻温度为23度，且𝑄𝑘=9,。B、温度计得到的此时温度值为25，且𝑅𝑘=4温度问题的Matlab实现温度问题的Matlab实现温度问题的Matlab实现可以看出增益𝐻𝑘=0.75，根据公式𝑥𝑘=𝐴(𝑘)𝑥𝑘−1+𝐻𝑘(z𝑘−𝐶𝑘𝐴𝑘𝑥𝑘−1)即现在的温度为𝑥𝑘=1∗23+0.75∗25−1∗1∗23=23.5°𝐶目录3.第一部分：维纳滤波原理及应用——蒋芮第二部分：常规Kalman滤波原理及应用——王丹第三部分：Kalman滤波处理语音信号——程谦1.2.目录3.背景介绍工程实现卡尔曼和维纳滤波比较1.2.目录3.背景介绍工程实现卡尔曼和维纳滤波比较1.2.背景介绍目的：从带噪语音信号中尽可能提取干净语音信号，提高信噪比，改善语音质量周围环境噪声干扰传输媒介电气设备背景介绍--语音特性语音信号是非平稳信号，但在短时间内，其频谱是稳定的。即在短时间内可以用平稳随机过程方法来分析语音信号。语音信号是随机的，可以利用许多统计分析特征进行分析。但由于语音信号非平稳、非遍历，因此长时间时域统计特性对语音增强算法的意义不大。在高斯模型的假设中，认为傅立叶展开系数是独立的高斯随机变量，均值为零，而方差是时变的。在有限帧长时这种高斯模型是一种近似的描述。背景介绍--噪声特性噪声通常可以定义为通信、测量以及其他信号处理过程中的无用信号成分。(1)周期性噪声其特点是频谱上有许多离散的线谱，主要来源于发动机等周期运转的机械设备。(2)脉冲噪声脉冲噪声表现为时域波形中突然出现的窄脉冲，主要来源于爆炸、撞击、放电及突发性干扰。(3)宽带噪声宽带噪声的来源很多，热噪声、气流噪声及各种随机噪声源，量化噪声都可视为宽带噪声。宽带噪声与语音信号在时域和频域上完全重叠，只有在无声期间，噪声分量才单独存在。对于平稳的高斯噪声，通常可以认为是高斯白噪声。而不具有白色频谱的噪声，可以进行白化处理或者采取特殊的处理方法。如本文后面介绍的建模的方法。语音信号全极点模型(Autoregressivemodel，AR模型)状态空间形式：卡尔曼模型带噪语音模型语音模型的假设1.语音信号x(n)在短时窗内(20～30ms)假定是短时平稳的，这使得在短时窗内，其A参数保持不变。2.W(n)、v(n)始终是均值为0、方差为𝛿𝑤2和𝛿𝑣2的不相关白噪声，并且在短时窗内其方差不变。3.随机系统初始状念x(0)与w(n)、v(n)不相关。卡尔曼模型卡尔曼模型卡尔曼模型n-1时刻对n时刻状态的预测值n时刻结合观测值对真实状态的估计预测的误差协方差矩阵滤波估计的误差协方差矩阵卡尔曼增益框架采样分帧卡尔曼滤波组帧y(n)𝒚𝒏(𝒏,𝒍)𝒙(𝒏)𝒙(𝒏,𝒍)参数估计语音分帧滑动窗帧内估计加窗合帧目录3.背景介绍工程实现卡尔曼和维纳滤波比较1.2.m1=wavread('bingyu.wav');y=m1(1:10240*20)';Fs=22050;x=y+0.6*randn(1,length(y));加载语音数据Fs为采样速率，可以听到的声音频率为20HZ~20kHZ。根据奈奎斯特采样定理，采样速率为原信号频谱两倍即可无失真恢复但人耳分辨率有限，一般取22050，44100为CD音质。figure(1);Fs=22050;plot((1:length(y(500:900)))/Fs,y(500:900),'r');holdon;plot((1:length(x(500:900)))/Fs,x(500:900),'g');原音频信号和加噪声后的信号00.0020.0040.0060.0080.010.0120.0140.0160.0180.02-3-2-10123原始信号加噪声信号Fs=22050;%信号采样频率bits=16;%信号采样位数N=256;%帧长m=N/2;%每帧移动的距离lenth=length(x);%输入信号的长度count=floor(lenth/m)-1;%处理整个信号需要移动的帧数p=11;%AR模型的阶数音频处理相关参数w=hamming(N);%汉明窗，汉明窗点数为256F=zeros(11,11);%转移矩阵F(1,2)=1;F(2,3)=1;F(3,4)=1;F(4,5)=1;F(5,6)=1;F(6,7)=1;F(7,8)=1;F(8,9)=1;F(9,10)=1;F(10,11)=1;定义窗函数和转移矩阵a=zeros(1,p);H=zeros(1,p);S0=zeros(p,1);P0=zeros(p);S=zeros(p);H(11)=1;s=zeros(N,1);G=H';P=zeros(p);矩阵初始化置零y_temp=0;y_temp=cov(x(1:7680));x_frame=zeros(256,1);x_frame1=zeros(256,1);T=zeros(lenth,1);求测试噪声方差x_frame=x((r-1)*m+1:(r+1)*m);分帧forr=1:countx_frame=zeros(256,1);x_frame1=zeros(256,1);T=zeros(lenth,1);ifr==1[a,VS]=lpc(x_fram