2022年数字信号处理实研究 .pdf

实验一基于 LMS算法的自适应滤波器设计一、自适应算法的概括包括最小均方算法LMS 、最小高阶均方算法LMP 、最小平方算法 OLS 、递推最小算法 RLS 。自适应算法主要根据滤波器输入的统计特性进行处理，存在开环算法和闭环算法；开环算法的控制输出仅取决于滤波器的输入和其他输入数据；闭环的控制输出则是滤波器输出及其他输入信号的函数。闭环控制利用输出反馈，不但能在滤波器输入信号变化时保持最佳输出，而且在某种程度上补偿滤波元件参数的变化和误差以及运算误差。二、自适应滤波器的结构自适应滤波器由参数可调的数字滤波器和自适应算法两部分组成。如图1.1所示为自适应滤波器的一般结构。 输入信号)(nx通过参数可调数字滤波器后产生输出信号（或响应）)(ny，将其与参考信号（或期望响应）)(nd进行比较，形成误差信号)(ne，)(ne通过某种自适应算法对滤波器参数进行调整，最终使)(ne得均方值最小。 所以，自适应滤波器实际上是一种能够自动调整本身参数的特殊的 Wiener 模型。)(nd自适应滤波器算法滤波器结构)(nx)(ne)(ny+-图 1.1 自适应滤波器的一般结构上面的自适应滤波器设计不需要知道关于输入信号和噪声的统计特性，能够在工作过程中估计出所需的统计特性，并以此为依据自动调整自己的参数，以达到最佳的滤波效果。 一旦输入信号统计特性发生变化，其又能跟踪这种变化， 自动调整参数，从而使滤波器性能达到最佳效果。三、滤波器采用的结构精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 20 页采用 FIR 横向滤波器（由于 IIR 滤波器存在稳定性问题）作为自适应滤波器结构，如图 1.2 所示。)(nx)(0nw)(1nw)(2nw)(1nwN1Z1Z.1Z)1(Nnx)(ny图 1.2 FIR 横向滤波器结构图中所示自适应滤波器的输入矢量：TNnxnxnxnX)1(),.,1(),()(，权重系数矢量：TNnwnwnwnW)(),.,(),()(110，即自适应滤波器的冲击响应。)(ny为自适应滤波器的输出矢量；n为时间序列；N为滤波器的阶数。其均方误差为式（ 1-1）所示：)11()()()(2neEnWV滤波器输出为式（ 1-2）所示：)21()()()()()(nWnXnXnWnyTT)(nx的自相关函数为式（ 1-3）所示：)31()()()(inxnxEirxx)(nx和期望输出信号)(nd的互相关函数为式（ 1-4）所示：)41()()(inxndErxy误差信号)(ne定义为期望输出)(nd与滤波器的实际输出之间的误差，即式（1-3）所示：)51 ()()()()()()(nXnWndnyndneT精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 20 页在)(nW是常数矢量的情况下，由式（1-1）可得出)61()()()(2)()()()()()(222nWRnWnWpndEnXnWndEneEnxxTTT上述为均方误差性能函数，其中xxR是输入矢量)(nx的自相关矩阵，P是输入矢量与期望响应的互相关矢量)()(ndnXEP为使均方误差最小化，对)(n关于)(nW求导使其值为 0，得到式（ 1-7）所示：)71()(22)()(nWRPWVnwxx在最小均方误差下，滤波器的权向量optW，一定满足正则方程：)81(PWRoptxx当xxR为满秩时，正则方程有唯一解，称为维纳解，即PRWxxopt1。则得到均方误差的最小值表达式为维纳误差，式（1-9）所示：)91()(2minoptTWPndE四、LMS 自适应滤波算法LMS 算法使用的是最陡下降法， 它是取单个误差样本的平方)(2ne的梯度作为均方误梯度的估计， 则由式（1-1）和式（1-5）可推出梯度公式为式 （1-10）所示。)101()()()(2)()(2)()(2)()()(2nWnXnXndnXnXnenWnenT由此得到一个新的权矢量递推公式，即LMS 算法的递推公式为式（ 1-11）所示。)111()()()()()() 1(21nenuXnWnunWnW五、LMS 算法仿真精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 20 页 1 、LMS算法步骤：（1） 、设置变量和参量：)(nx为输入向量，或称为训练样本，)(nW为权值向量，)(nd为期望输出；)(ny为实际输出，为学习速率，n为迭代次数。（2） 、初始化，赋给)(nW各一个较小的随机非零值，（3） 、对于一组输入样本)(nx和对应的期望输出)(nd，计算)()()()(nWnxndneT)()()()1(nenxnWnW（4） 、判断)(ne是否满足条件，若满足算法结束，若否则1n，转入第 3步继续执行。2、基本 LMS算法程序，见附录 1。3、程序仿真结果分析。01002003004005006007008009001000-202信 号 s 时 域 波 形01002003004005006007008009001000-202信 号 s加 噪 声 后 的 时 域 波 形01002003004005006007008009001000-505自 适 应 滤 波 后 输 出 的 时 域 曲 线010020030040050060070080090000.51收 敛 曲 线u=0.0002.K=128时 的 滤 波 情 况图1.3 步长因子u=0.0002 时的滤波情况结果分析：从图中可以看到，在加入信噪比为5dB噪声后输出的正弦波形有明显的抖动失真， 在通过步长因子 u=0.0002，时域抽头阶数 K=128的自适应滤波后，一开始由于迭代次数较少，权重次数W还没有得到完全修正，误差较大，在N在0-120之间输出的波形还有较大的失真，经过迭代次数的增加，均方误差逐渐趋于 0，在 N=350 之后可以看到输出的波形有了明显的改善，从而可以得出，自适应滤波器运用LMS算法能对加噪的信号进行修正， 得到的输出信号有较高的精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 20 页精确度。01002003004005006007008009001000-202信 号 s时 域 波 形01002003004005006007008009001000-505信 号 s加 噪 声 后 的 时 域 波 形01002003004005006007008009001000-505自 适 应 滤 波 后 输 出 的 时 域 曲 线010020030040050060070080090000.51收 敛 曲 线u=0.005,K=128 时 的 滤 波 情 况图1.4 步长因子u=0.005 时的滤波情况结果分析，在步长因子 u 增加后，可以看到收敛速度加快， 自适应时间缩短，收敛曲线经过较少的时间就趋于0，但是从滤波输出曲线可以看到在N=0-100之间输出波形并没有得到很好的矫正， 这说明步长因子增加后滤波性能有一定程度上的降低。01002003004005006007008009001000-202信 号 s 时 域 波 形01002003004005006007008009001000-505信 号 s加 噪 声 后 的 时 域 波 形01002003004005006007008009001000-505自 适 应 滤 波 后 输 出 的 时 域 曲 线0100200300400500600700800900100000.51收 敛 曲 线u=0.0002,K=32 时 的 滤 波 情 况图1.5 滤波器阶数K=32 时的滤波情况结果分析：当滤波器阶数K减少后，比较图 1.3 可以看到其收敛速度变慢。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 20 页从图中可以看出，误差曲线收敛缓慢，知道在N接近1000时，还有较小的误差。但当阶数 K增大道一定程度， 收敛速度变化不明显， 且可能引起系数迭代过程不收敛。六、小结LMS算法最核心思想是用平方误差代替均方误差，自适应步长的选择u越大，自适应时间越短，信噪比当信噪比升高时，LMS 算法性能将急剧恶化；反映了解调器的抗噪声性能， LMS算法的优点为计算简单，每次迭代只需要2N+1 次计算，缺点是收敛缓慢， 另外一种自适应算法RLS算法收敛速度和均衡效果都比LMS算法好，但其计算复杂度高， 一般我们根据实际需要选择相应的自适应算法。七、参考文献1 张立萍.LMS 自适应滤波器的 Matlab设计与仿真J. 赤峰学院学报（自然科学版） ,2010,26(5) 2 贺宽 , 黄涛 . 基于Matlab的自适应滤波器设计J. 武汉理工大学学报( 信息与管理工程版) ,2008,30(1):70-73. 附录1 基本 LMS 算法滤波器程序% %最小均方算法 LMS 的程序%程序部分参考自然科学报LMS 自适应滤波器的Matlab设计与仿真%备注：%编写：%参数说明：W(n): 滤波器的权向量序列% clc clear all K=128; %时域抽头 LMS 算法滤波器阶数可变： K=32 u=0.0002; %滤波器步长因子可变： u=0.005 N=1000; %时域 t的范围 ,抽样点数a=200; %统计仿真次数a all=zeros(a,N-K); %将每次独立循环的误差结果存于矩阵pp中，以便后面对其平均for b=1:a %求每次的误差循环t=1:N; s=sin(0.05*pi)*t);%标准正弦信号figure(1) subplot(4,1,1) plot(t,s); 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 20 页title( 信号 s时域波形 ); set(gca,YLim,-2,2); % 设置 y轴的范围为 -2 2 sn=awgn(s,3); %加入均值为零的高斯白噪声，信噪比为3dB subplot(4,1,2) plot(sn); title( 信号 s加噪声后的时域波形); set(gca,YLim,-5,5); %参数初始化w=zeros(1,K); %设置抽头加权初值y=zeros(1,N); %设置输出初值y(1:K)=sn(1:K);% 将输入信号 sn的前 k个值作为输出 y的前 k个值e=zeros(1,N); %误差信号%用LMS 算法迭代滤波for i=(K+1):N SN=sn(i-K+1):(i); % 输入初始值y(i)=w*SN; %输出值e(i)=s(i)-y(i); %s(i)减去滤波器输出信号后得到的误差信号w=w+u*e(i)*SN; %LMS 算法的递推公式end all(b,:)=(e(K+1:N).2; % end for c=1:N-K ai(c)=sum(all(:,c)/a; % 求误差的统计平均end subplot(4,1,3) plot(y); title( 自适应滤波后输出的时域曲线); subplot(4,1,4) t=1:N-K; plot(t,ai); title( 收敛曲线 ); xlabel(u=0.0002 ，K=128 时的收敛情况 ); % x轴注解hold off %将每次循环的图形显示结果保存下来精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 20 页实验二功率谱估计与MATLAB仿真一、信号功率谱介绍对于常见的具有各态历经的平稳随机信号，不可能用清楚的数学关系式来描述，但可以利用给定的N 个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计 (PSD)。对确定性信号傅里叶变换是在频域分析研究的理论基础，但对于随机信号， 其傅里叶变换并不存在， 因此研究其功率谱对分析信号有重要作用。按照 Weiner-Khintchine 定理，信号的功率谱和其自相关函数服从一对傅里叶变换公式，如式2-1、式 2-2、式 2-3 所示。)32()()()()22(d)(21)() 12()()(mnxnxEmrweePmremrePxxjwnjwxxxxmjwnxxjwxx二、功率谱估计方法精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 20 页功率谱的估计方法有很多种， 一般分成两大类， 一类是经典谱估计， 也称线性谱估计；另一类是现代谱估计，也称非线性谱估计。其中，经典功率谱估计方法分为：相关函数法 (BT 法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法。现代谱估计方法分为：AR 模型谱估计、 Welch 法谱估计、 MUSIC 法谱估计。1、相关函数法 (BT 法) 该方法先由序列 x(n)估计出自相关函数R(n)，然后对 R(n)按式(2-1)进行傅立叶变换，便得到 x(n)的功率谱估计。一般我们采用有偏自相关函数估计，如式2-4 所示。)42()()(1)(1|0mNnxxmnxnxNmr2、周期图法周期图法是把随机序列x(n)的 N 个观测数据视为一能量有限的序列，直接计算 x(n)的离散傅立叶变换，得x(k)，然后再取其幅值的平方，并除以N，作为序列 x(n)真实功率谱的估计，式 (2-5)所示。)52()(1)(210NnjwnjwxxenxNeP计算框图如图 2.1 所示：观测数据 x(n)FFT取模的平方N1)(jwxxeP图 2.1 周期图法计算功率谱框图3、平均周期图法对于周期图的功率谱估计，当数据长度N 太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好， 因此需要改进。 两种改进的估计法是平均周期图法和平滑平均周期图法。平均周期图法，也称Bartlett 法，方法是将 N 点的有限长序列 x(n)分段求周期图再平均。4、Welch 法Welch 法对 Bartlett 法进行了两方面的修正，一是选择适当的窗函数w(n)，精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 20 页并在周期图计算前直接加进去， 加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。5、AR 模型谱估计AR 模型功率谱估计又称为自回归模型，它是一个全极点的模型，要利用AR 模型进行功率谱估计须通过列文森Levinson_dubin 递推算法由 YuleWalker方程求得 AR 的参数：p212.,。伯格 (Burg)递推法则由信号观测数据直接计算 AR 模型的参数。6、MUSIC 法谱估计MUSIC 作为一种高分辨率的子空间方法，首先其主要应用于离散谱的估计，比如混叠在一起的单频信号； 其频谱峰值反映了这些主要信号成分所在的频率位置，但是其并不能反映各信号成分之间的幅度比值（相对强度），也反映不出信噪比水平， 所以 MUSIC 算法所得到的 “ 谱” 被称为伪谱。 可直接在 MATLAB中调用函数 pmusic 来实现 MUSIC 算法的功率谱估计。三、功率谱估计 MATLAB仿真实现 （相关代码见附录2）1、相关函数法 (BT 法) 实现步骤： (1)、求信号序列 x(n)的自相关函数；cxn=xcorr(xn,unbiased); (2)、对自相关函数进行傅里叶变换；CXk=fft(cxn,nfft); (3)、求傅里叶变换CXk 幅值的平方，并转换成dB 表示，即对复数信号操作如下：Px_dB = CXk*conj(CXk); %求包络的平方Px_dB = 10*log10(Px_dB/nfft+eps); 仿真信号：采用课本 170 页的上机作业 1 中的信号，即正弦波加正态零均值白噪声，信噪比为10dB，信号频率为2kHz，取样频率为100kHz，信号数据长度 N=256。仿真结果：其信号时域图如图2.2，功率谱估计如图2.3 所示；精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 20 页00.511.522.53x 10-3-1-0.500.51 t x时 域 信号00.511.522.53x 10-3-2-1012 t xn加 噪 信号-5-4-3-2-1012345x 104-80-70-60-50-40-30-20-1001020 FrequencyHz MagnitudedB图 2.2 信号时域图图 2.3 BT 法功率谱估计图2、周期图法仿真结果：其信号时域图如图2.2，功率谱估计如图2.5 所示。-5-4-3-2-1012345x 104-40-30-20-1001020 FrequencyHz MagnitudedB-5-4-3-2-1012345x 104-20-15-10-505101520 FrequencyHz MagnitudedB图 2.4 周期图法功率谱估计图 2.5 修正的周期图法功率谱估计3、平均周期图法将 N 点的有限长序列x(n)分段求周期图再平均，仿真结果：其信号时域图如图 2.2，功率谱估计如图 2.6 所示。4、Welch 法选择了适当的窗函数w(n)，在分段时使各段之间有重叠。仿真结果：其信号时域图如图2.2，功率谱估计如图2.7 所示。 （分段长度为L=64，加汉明窗）精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 20 页-5-4-3-2-1012345x 104-30-25-20-15-10-505 FrequencyHz MagnitudedB05101520253035404550-65-60-55-50-45-40-35-30-25Frequency (kHz)Power/frequency(dB/Hz)Burg Power Spectral Density Estimate图 2.7 Welch 法功率谱估计图 2.8 AR 模型谱估计5、AR 模型谱估计在MATLAB 中，主要由基于 AR模型的功率谱估计 pyulear函数、 pburg函数、pcov函数、pmov函数和pmem函数，此处采用基于 burg算法求出的 AR 功率谱曲线，如图2.8 所示。四、小结通过上面对单音正弦加噪信号的不同功率谱估计方法仿真分析，我们可以看出，相关函数 (BT)法和周期图法得到的谱线图离散性大，曲线粗糙，方差较大。平均周期图法的收敛性较好， 曲线比较平滑， 估计的结果方差较小， 但是其功率谱主瓣较 BT 法更宽，分辨率减少了。采用Welch法可以看到其功率谱曲线平滑，由于用了窗函数 (汉明窗 ) 极大地压低了旁瓣， 但是从图中可以看出， 其主瓣宽度加宽了，功率谱波峰变宽了，降低了分辨率。采用AR 模型的 burg算法求得的功率谱曲线图主瓣窄，频谱方差小，分辨率高。注意：在分析数据后面填充0，是分析点数增多，可以对信号的频谱观测变得更细，但是根本上不能提高频率分辨率，因为在数据后面加0，其傅里叶变换结果没有受到影响，因此也不能提高分辨率。五、参考文献1 王凤瑛 ,张丽丽 .功率谱估计及其MATLAB仿真 J.微计算机信息,2006,22(31). 2 王福杰，潘宏侠.MATLAB 中几种功率谱估计函数的比较分析与选择J.电子产品可靠性与环境试验 ,2009,27(6). 3 丁玉美等 . 数字信号处理 . 西安：西安电子科技大学出版社，2002. 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 20 页附录2 功率谱估计仿真程序% %功能描述：功率谱估计(Power Spectrum Density Estimation) ，% clc clear all close all %= 信号产生 = f = 2000; %信号频率为 2kHZ fs = 10e4; %采样频率为 100kHz N = 256; %数据长度t = (0:N-1)/fs; x = cos(2*pi*f*t); xn = awgn(x,10); %加10dB噪声%= 各种谱估计方法= Type = 5; % 取不同值采用不同方法谱估计switch Type case 1 % 相关函数 (BT) 法nfft = 512; % FFT点数cxn = xcorr(xn,unbiased); %计算序列的自相关函数CXk = fft(cxn,nfft); f1 = 0:nfft-1/nfft; Px_dB = 10*log10(CXk.*conj(CXk)/nfft+eps); Px_dB = fftshift(Px_dB); figure(1) subplot(211) plot(t,x,r),xlabel( t ),ylabel( x ), title( 时域信号 ),grid on subplot(212) plot(t,xn),xlabel( t ),ylabel( xn ), title( 加噪信号 ),grid on figure(2) plot(fs*f1-fs/2,Px_dB,b), xlabel( FrequencyHz ),ylabel(Magnitude dB),grid on case 2 %周期图法nfft = 512; f1 = 0:nfft-1/nfft; Px = fft(xn,nfft); Px_dB = 10*log10(abs(Px).2/nfft+eps); Px_dB = fftshift(Px_dB); plot(fs*f1-fs/2,Px_dB,b), xlabel( FrequencyHz ),ylabel(Magnitude dB),grid on case 3 % 修正的周期图法L = 32; % 数据分段长度(书中要求 ) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 20 页M = N/L; % 分段数n = 1 : L; L_nfft = 512; %每一段的 FFT点数sum_P = zeros(1,L_nfft); % 对每一段数据进行FFT，并求和for k = 1 : M x = xn(n:(k-1)*L); sum_P = sum_P + (abs(fft(x,L_nfft).2/L_nfft; %对每段求 FFT end f1 = 0:L_nfft-1/L_nfft; Px_dB = 10*log10(sum_P+eps); Px_dB = fftshift(Px_dB); plot(fs*f1-fs/2,Px_dB,b),xlabel( FrequencyHz ), ylabel(Magnitude dB),grid on case 4 % Welch法%对数据进行补零M = ceil(log2(length(xn); N = 2M; % 求xn的长度对应的 2的最低幂次 m %若xn的长度不是 2的幂，补零到 2的整数幂if length(xn)=0.5L, 即最多重叠 50% D = L - K; % 重叠部分长度 D Mmax = floor(N-L)/K+1); %分段数计算M = Mmax ; % 分段数windowb = hamming(L); % 汉 明 窗 系 数blackman(L)=blackman窗 ，boxcar(L)= 矩形窗% beta =20; % windowb = kaiser(L,beta); % 加窗提高分辨率U = sum(windowb.2)/L; n = 1:L; L_nfft = 512; sumP = zeros(1,L_nfft); for i = 1:M X = xn(n+(i-1)*K); % 对每段数据进行提取X = X.*windowb; % 数据加窗sumP = sumP+(abs(fft(X,L_nfft).2/L_nfft ; %对每段求 FFT end sumP = sumP/(M*U); Px_dB = 10*log10(sumP+eps); Px_dB = fftshift(Px_dB); 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 20 页f1 = (0:L_nfft-1)/L_nfft; plot(fs*f1-fs/2,Px_dB,b), xlabel( FrequencyHz ),ylabel(Magnitude dB),grid on case 5 % AR模型谱估计基于 burg算法求出的 AR 功率谱曲线order = 50; range = half; nfft = 512; figure pburg(xn,order,nfft,fs,range); end精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 20 页实验三卡尔曼滤波器的设计一、卡尔曼滤波器介绍卡尔曼滤波是用状态空间法描述系统的，由状态方程和测量方程所组成。卡尔曼滤波用前一个状态的估计值和最近一个观测数据来估计状态变量的当前值，并以状态变量的估计值的形式给出。卡尔曼滤波的状态方程和测量方程分别为式(3-1)、式(3-2)所示。)23() 13(1kkkkkkkkvxCywxAx其中， k 表示时间，输入信号kw是一白噪声，输出信号的观测噪声kv也是一个白噪声， A表示状态变量之间的增益矩阵，可随时间变化，kA表示第 k 次迭代的取值， C 表示状态变量与输出信号之间的增益矩阵，可随时间变化， 其信号模型如图 3.1 所示( k 用1k代替)。1Z+1kw1kAkC+1kxkxkvky图 3.1 卡尔曼滤波器的信号模型二、卡尔曼滤波的递推算法卡尔曼递推公式如下式 (3-3a)、(3-3b)、(3-3c)及(3-3d)所示。)33()()33()33()()33()(11111dPCHIPcQAPAPbRCPCCPHaxACyHxAxkkkkkTkkkkkTkkkTkkkkkkkkkkk假设初始条件11,kkkkkkkPxyRQCA已知，其中var,0000 xPxEx，那么递推流程见图3.2 所示。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 16 页，共 20 页1kx1kP式(3-3c)kP式(3-3b)式(3-3a)kH式(3-3d)kxkP图 3.2 卡尔曼滤波递推过程三、卡尔曼滤波仿真1、 仿真条件：采用课本 61 页上机作业 1 中的条件，|)(erx，sTs02.0。2、假设一连续平稳的随机信号)(tx，自相关函数|)(erx，信号)(tx为加性噪声所干扰， 噪声是白噪声， 以此并结合观测值， 设计卡尔曼滤波递推程序，估计信号)(tx的波形。信号)(tx只在加性白噪声下通过， 因此，输出kkkvxy，所以1kC，我们可以求出各个参数如下，由已知条件可得：)53()1)(1(1)43()()(1112|zezeezezmrzSmmmmmxxxx)63()()()(12zBzBzSwxx所以：)73()1()1()(111111eAnwnxenxzezBkz由取样时间可得：02. 0eeAsTk。又有式 (3-8)式：ssssssTxxTxxTxxTxxTTxxwkerererernxenxnxenxEnwnwErQ22210 1 10)83()()1()() 1()()(0因此04.0211eeQsTk，又由测量方程可得：)93(1)0(2vvvkSR我们采用的信号初值01x， 则1var11xP。 因此将解得的参数代入式 (3-3)可以得到卡尔曼递推公式，如式(3-10)所示。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 17 页，共 20 页)103()()103(1)103()1()103()(04. 0104. 01102. 0102. 0dP?HIPcePePb? P?PHaxeyHxexkkkkkkkkkkkkk由上式可得递推流程为图3.3 所示，则代入初值迭代后就可得到结果。1kP式(3-10c)kP式(3-10b)式(3-10a)kH式(3-10d)kxkP图 3.2 仿真递推过程3、由书中的测量值的离散值)(kz为：-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，则通过 MATLAB 实现，程序见附录3。仿真结果，如图 3.3 所示。00.10.20.30.40.50.60.70.80.9-20-15-10-5051015 t/s z(t),x(t)x(t)的 估 计 波 形 与 观 测 值 z(t) 波 形 比 较观 测 数 据 z(t)信 号 估 计 值 x(t)图 3.3 卡尔曼滤波估计信号与观测信号比较结果分析：从图3.3 可以看出，卡尔曼滤波输出的估计信号)(tx与实际观测到的离散值)(tz还是存在一定的误差， 卡尔曼滤波是从初始状态就采用递推方法精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 20 页进行滤波，那么在初值迭代后的一段时间内可能会出现较大的误差，随着迭代进行，各参数逐渐趋于稳定，后面的估计值与观察值的误差就减少了。四、小结通过对卡尔曼滤波的仿真，我们可以知道，卡尔曼滤波要求已知前一个时刻的状态估计值)1(kx和当前的观测值ky，由状态方程和观测方程递推得到结果。维纳滤波的解以)(zH的形式给出，卡尔曼滤波是以状态变量的估计值给出解的形式， 它们都采用均方误差最小准则，但是卡尔曼有一个过渡过程，其解与维纳滤波不完全相同，但达到稳态后结果是相同的。五、参考文献1 丁玉美等 . 数字信号处理. 西安：西安电子科技大学出版社，2002. 2 百度网页 http:/ 附录 3 卡尔曼滤波仿真程序% %设计卡尔曼滤波估计信号的波形% clc clear all close all %=参数数值 = Ak = exp(-0.02); %各系数由前面确定；Ck = 1; Qk = 1-exp(-0.04); Rk = 1; p(1) = 1; 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 为测量出来的离散值 = 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 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); %增益方程；精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 20 页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 ),ylabel( z(t),x(t) ); title(x(t) 的估计波形与观测值 z(t)波形比较 ); legend( 观测数据 z(t) , 信号估计值 x(t) ,2); grid on; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 20 页，共 20 页