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

《2022年数字信号处理实研究.docx》由会员分享，可在线阅读，更多相关《2022年数字信号处理实研究.docx（39页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选学习资料 - - - - - - - - - 试验一 基于 LMS算法的自适应滤波器设计一、自适应算法的概括包括最小均方算法LMS、最小高阶均方算法LMP、最小平方算法 OLS、递推最小算法 RLS；自适应算法主要依据滤波器输入的统计特性进行处理,存在开环算法和闭环算法；开环算法的掌握输出仅取决于滤波器的输入和其他输入数据；闭环的掌握输出就是滤波器输出及其他输入信号的函数；闭环掌握利用输出反馈,不但能在滤波器输入信号变化时保持正确输出,元件参数的变化和误差以及运算误差；二、自适应滤波器的结构而且在某种程度上补偿滤波自适应滤波器由参数可调的数字滤波器和自适应算法两部分组成；如图 1.1所示为

2、自适应滤波器的一般结构； 输入信号 x n 通过参数可调数字滤波器后产生输出信号（或响应）y n ,将其与参考信号（或期望响应）d n 进行比较,形成误差信号 e n ,e n 通过某种自适应算法对滤波器参数进行调整,最终使 e n 得均方值最小； 所以,自适应滤波器实际上是一种能够自动调整本身参数的特别的 Wiener 模型；dn滤波器结构+-e n xny n自适应滤波器算法图 1.1 自适应滤波器的一般结构上面的自适应滤波器设计不需要知道关于输入信号和噪声的统计特性,能够在工作过程中估量出所需的统计特性,并以此为依据自动调整自己的参数,以达到正确的滤波成效； 一旦输入信号统计特性发生变化

3、,动调整参数,从而使滤波器性能达到正确成效；三、滤波器采纳的结构其又能跟踪这种变化, 自名师归纳总结 - - - - - - -第 1 页,共 20 页精选学习资料 - - - - - - - - - 采纳 FIR 横向滤波器（由于 IIR 滤波器存在稳固性问题）作为自适应滤波器结构,如图 1.2 所示；w 0 nxnZ1Z1.Z1xnN1 w 1nw 2nwN1 nyn 图 1.2 FIR 横向滤波器结构图中所示自适应滤波器的输入矢量：X n x n , x n 1 ,., x n N 1 T,权重系数矢量：W n w 0 n , w 1 n ,., w N 1 n T,即自适应滤波器的冲击

4、响应；y n 为自适应滤波器的输出矢量；n为时间序列；N为滤波器的阶数；其均方误差为式（ 1-1）所示：VWn Ee 2 n 11 滤波器输出为式（ 1-2）所示：T Ty n W n X n X n W n 1 2 x n 的自相关函数为式（ 1-3）所示：rxx i E x n x n i 1 3 x n 和期望输出信号 d n 的相互关函数为式（ 1-4）所示：rxy E d n x n i 1 4 误差信号 e n 定义为期望输出 d n 与滤波器的实际输出之间的误差,即式（1-3）所示：名师归纳总结 e n d ny ndn WnTX n15第 2 页,共 20 页- - - - -

5、 - -精选学习资料 - - - - - - - - - 在Wn 是常数矢量的情形下,由式（1-1 ）可得出2 T 2 n E e n E d n W n X n 2 T TE d n 2 p W n W n R xx W n 1 6 上述为均方误差性能函数,其中 R是输入矢量 x n 的自相关矩阵, P 是输入矢量与期望响应的相互关矢量 P E X n d n 为使均方误差最小化,对n关于 W n 求导使其值为 0,得到式（ 1-7）所示：V W 2 P 2 R xx W n 1 7 w n 在最小均方误差下,滤波器的权向量 W opt,肯定满意正就方程：R xx W opt P 1 8 当

6、 R为满秩时,正就方程有唯独解,称为维纳解,即 W opt R 1xx P；就得到均方误差的最小值表达式为维纳误差,式（1-9）所示：min E d 2 n P TW opt 1 9 四、LMS 自适应滤波算法LMS 算法使用的是最陡下降法, 它是取单个误差样本的平方e 2 n的梯度作为均方误梯度的估量, 就由式（1-1）和式（1-5）可推出梯度公式为式 （1-10 ）所示；ne2n 2e n Xn 2Xndn2XnXTnWn 110Wn由此得到一个新的权矢量递推公式,即 所示；LMS 算法的递推公式为式（ 1-11）Wn1Wn 1 2u n Wn uXn e n 111 五、LMS 算法仿真

7、名师归纳总结 - - - - - - -第 3 页,共 20 页精选学习资料 - - - - - - - - - 1 、LMS算法步骤：（1）、设置变量和参量：x n 为输入向量,或称为训练样本,W n 为权值向量,d n 为期望输出；y n 为实际输出,为学习速率, n为迭代次数；（2）、初始化,赋给 W n 各一个较小的随机非零值,（3）、对于一组输入样本 x n 和对应的期望输出 d n ,运算e n d n x T n W n W n 1 W n x n e n （4）、判定 e n 是否满意条件,如满意算法终止,如否就 n 1,转入第 3步连续执行；2、基本 LMS算法程序,见附录

8、1；3、程序仿真结果分析；20100200300信 号 s 时 域 波 形70080090010000400500600-2信 号 s加 噪 声 后 的 时 域 波 形20-2501002003004005006007008009001000自 适 应 滤 波 后 输 出 的 时 域 曲 线001002003004005006007008009001000-510100200收 敛 曲 线8009000.53004005006007000图1.3 u=0.0002.K=128时 的 滤 波 情 况步长因子 u=0.0002 时的滤波情形结果分析：从图中可以看到,在加入信噪比为5dB噪声后输出的

9、正弦波形有明显的抖动失真, 在通过步长因子 u=0.0002,时域抽头阶数 K=128的自适应滤波后,一开头由于迭代次数较少,权重次数W仍没有得到完全修正,误差较大,在N在0-120之间输出的波形仍有较大的失真,经过迭代次数的增加,均方误差逐 渐趋于 0,在 N=350之后可以看到输出的波形有了明显的改善,从而可以得出,名师归纳总结 自适应滤波器运用LMS算法能对加噪的信号进行修正, 得到的输出信号有较高的第 4 页,共 20 页- - - - - - -精选学习资料 - - - - - - - - - 精确度；20100200信 号 s时 域 波 形800900100003004005006

10、00700-250100200信 号 s加 噪 声 后 的 时 域 波 形80090010000300400500600700-550100200自 适 应 滤 波 后 输 出 的 时 域 曲 线80090010000300400500600700-510100200收 敛 曲 线7008009000.53004005006000u=0.005,K=128 时 的 滤 波 情 况图1.4 步长因子 u=0.005 时的滤波情形结果分析,在步长因子 u 增加后,可以看到收敛速度加快, 自适应时间缩短,收敛曲线经过较少的时间就趋于0,但是从滤波输出曲线可以看到在N=0-100之间输出波形并没有得到

11、很好的矫正, 这说明步长因子增加后滤波性能有肯定程度上的降低；20100200信 号 s 时 域 波 形80090010000300400500600700-250100200信 号 s加 噪 声 后 的 时 域 波 形80090010000300400500600700-550100200自 适 应 滤 波 后 输 出 的 时 域 曲 线80090010000300400500600700-510100200收 敛 曲 线80090010000.53004005006007000u=0.0002,K=32 时 的 滤 波 情 况图1.5 滤波器阶数 K=32 时的滤波情形结果分析：当滤波器阶

12、数K削减后,比较图 1.3 可以看到其收敛速度变慢；名师归纳总结 - - - - - - -第 5 页,共 20 页精选学习资料 - - - - - - - - - 从图中可以看出,误差曲线收敛缓慢,知道在N接近1000时,仍有较小的误差；但当阶数 K 增大道肯定程度, 收敛速度变化不明显, 且可能引起系数迭代过程不 收敛；六、小结LMS算法最核心思想是用平方误差代替均方误差,自适应步长的挑选 u越大,自适应时间越短,信噪比当信噪比上升时,LMS算法性能将急剧恶化；反映明白调器的抗噪声性能, LMS算法的优点为运算简洁,每次迭代只需要2N+1次运算,缺点是收敛缓慢, 另外一种自适应算法RLS算

13、法收敛速度和均衡成效都比LMS算法好,但其运算复杂度高, 一般我们依据实际需要挑选相应的自适应算法；七、参考文献1 张立萍 .LMS 自适应滤波器的 Matlab设计与仿真J. 赤峰学院学报（自然科学版） ,2022,265 2 贺宽 , 黄涛 . 基于Matlab 的自适应滤波器设计J. 武汉理工高校学报 信息与治理工程版 ,2022,301:70-73. 附录 1 基本 LMS 算法滤波器程序% %最小均方算法 LMS 的程序%程序部分参考自然科学报LMS 自适应滤波器的Matlab 设计与仿真%备注：%编写：%参数说明：Wn: 滤波器的权向量序列 % clc clear all K=12

14、8; %时域抽头 LMS 算法滤波器阶数可变： K=32 u=0.0002; %滤波器步长因子可变： u=0.005 pp中,以便后面对其平均N=1000; %时域 t的范畴 ,抽样点数a=200; %统计仿真次数 a all=zerosa,N-K; %将每次独立循环的误差结果存于矩阵for b=1:a %求每次的误差循环t=1:N; s=sin0.05*pi*t;% 标准正弦信号figure1 subplot4,1,1 plott,s; 名师归纳总结 - - - - - - -第 6 页,共 20 页精选学习资料 - - - - - - - - - title 信号 s时域波形 ; setg

15、ca,YLim,-2,2; % 设置 y轴的范畴为 -2 2 sn=awgns,3; %加入均值为零的高斯白噪声,信噪比为3dB subplot4,1,2 plotsn; title 信号 s加噪声后的时域波形 ; setgca,YLim,-5,5; %参数初始化 w=zeros1,K; %设置抽头加权初值 y=zeros1,N; %设置输出初值 y1:K=sn1:K;% 将输入信号 sn的前 k个值作为输出 y的前 k个值 e=zeros1,N; %误差信号 %用LMS 算法迭代滤波 for i=K+1:N SN=sni-K+1:i; % 输入初始值 yi=w*SN; %输出值 ei=si-

16、yi; %si 减去滤波器输出信号后得到的误差信号 w=w+u*ei*SN; %LMS 算法的递推公式 end allb,:=eK+1:N.2; % end for c=1:N-K aic=sumall:,c/a; % 求误差的统计平均 end subplot4,1,3 ploty; title 自适应滤波后输出的时域曲线 ; subplot4,1,4 t=1:N-K; plott,ai; title 收敛曲线 ; 名师归纳总结 xlabelu=0.0002 ,K=128 时的收敛情形 ; % x轴注解第 7 页,共 20 页hold off %将每次循环的图形显示结果储存下来- - - -

17、- - -精选学习资料 - - - - - - - - - 试验二功率谱估量与MATLAB仿真一、信号功率谱介绍对于常见的具有各态历经的平稳随机信号,不行能用清晰的数学关系式来描述,但可以利用给定的 N 个样本数据估量一个平稳随机信号的功率谱密度叫做功率谱估量 PSD；对确定性信号傅里叶变换是在频域分析讨论的理论基础,但对于随机信号, 其傅里叶变换并不存在, 因此讨论其功率谱对分析信号有重要作用；依据 Weiner-Khintchine 定理,信号的功率谱和其自相关函数听从一对傅里叶变换公式,如式2-1、式 2-2、式 2-3 所示；21P xxejwrxxm ejwnmrxxm 1P xxe

18、jwejwndw22 2rxxm Exnxnm 23二、功率谱估量方法名师归纳总结 - - - - - - -第 8 页,共 20 页精选学习资料 - - - - - - - - - 功率谱的估量方法有许多种, 一般分成两大类, 一类是经典谱估量, 也称线性谱估量；另一类是现代谱估量,也称非线性谱估量；其中,经典功率谱估量方法分为：相关函数法 BT 法、周期图法以及两种改进的周期图估量法即平均周期图法和平滑平均周期图法；现代谱估量方法分为：谱估量、 MUSIC 法谱估量；1、相关函数法 BT 法 AR 模型谱估量、 Welch 法该方法先由序列 xn估量出自相关函数 Rn,然后对 Rn按式 2

19、-1进行傅立 叶变换,便得到 xn的功率谱估量；一般我们采纳有偏自相关函数估量,如式 2-4 所示；rxxm 1Nn|m|1nx nm 240xN2、周期图法周期图法是把随机序列xn的 N 个观测数据视为一能量有限的序列,直接运算 xn的离散傅立叶变换,得xk,然后再取其幅值的平方,并除以N,作为序列 xn真实功率谱的估量,式 2-5所示；Pxxejw1N1xn ejwn225 Nn0运算框图如图 2.1 所示：观测数据 xnFFT取模的平方1NPxx ejw图 2.1 周期图法运算功率谱框图3、平均周期图法对于周期图的功率谱估量,当数据长度N 太大时,谱曲线起伏加剧,如N太小,谱的辨论率又不

20、好, 因此需要改进； 两种改进的估量法是平均周期图法和平滑平均周期图法；平均周期图法,也称 列 xn分段求周期图再平均；4、Welch 法Bartlett 法,方法是将 N 点的有限长序名师归纳总结 Welch 法对 Bartlett 法进行了两方面的修正,一是挑选适当的窗函数wn,第 9 页,共 20 页- - - - - - -精选学习资料 - - - - - - - - - 并在周期图运算前直接加进去, 加窗的优点是无论什么样的窗函数均可使谱估量非负；二是在分段时,可使各段之间有重叠,这样会使方差减小；5、AR 模型谱估量AR 模型功率谱估量又称为自回来模型,它是一个全极点的模型,要利用

21、AR 模型进行功率谱估量须通过列文森Levinson_dubin 递推算法由 YuleWalker方程求得 AR 的参数：2,1,2.p；伯格 Burg递推法就由信号观测数据直接运算 AR 模型的参数；6、MUSIC 法谱估量MUSIC 作为一种高辨论率的子空间方法,第一其主要应用于离散谱的估计,比如混叠在一起的单频信号； 其频谱峰值反映了这些主要信号成分所在的频率位置,但是其并不能反映各信号成分之间的幅度比值（相对强度）,也反映不出信噪比水平, 所以 MUSIC 算法所得到的 “谱”被称为伪谱； 可直接在 MATLAB中调用函数 pmusic 来实现 MUSIC 算法的功率谱估量；三、 功率

22、谱估量 MATLAB仿真实现 （相关代码见附录2）1、相关函数法 BT 法 实现步骤： 1、求信号序列 xn的自相关函数；cxn=xcorrxn,unbiased; 2、对自相关函数进行傅里叶变换；CXk=fftcxn,nfft; 3、求傅里叶变换CXk 幅值的平方,并转换成dB 表示,即对复数信号操作如下：Px_dB = CXk*conjCXk; %求包络的平方Px_dB = 10*log10Px_dB/nfft+eps; 仿真信号：采纳课本 170 页的上机作业 1 中的信号, 即正弦波加正态零均值白噪声,信噪比为10dB,信号频率为2kHz,取样频率为100kHz,信号数据长度 N=25

23、6；名师归纳总结 仿真结果：其信号时域图如图2.2,功率谱估量如图2.3 所示；第 10 页,共 20 页- - - - - - -精选学习资料 - - - - - - - - - 1时 域 信号20xn x0.500.511.522.5x 103 dB agnitude M10-4-3-2-101234x 10500-0.5-10-1-20 t -3-30200.5加 噪 信号2.5x 103-401-500-6011.52-1-70-2-80-5 t -3 FrequencyHz 4图 2.2 信号时域图图 2.3 BT 法功率谱估量图2、周期图法仿真结果：其信号时域图如图2.2,功率谱估

24、量如图2.5 所示；2020101510dB agnitude M0dB agnitude M5-100-20-5-10-30-4-3-2-101234x 105-15-4-3-2-101234x 105-40-5-20-5 FrequencyHz 4 FrequencyHz 4图 2.4 周期图法功率谱估量 3、平均周期图法图 2.5 修正的周期图法功率谱估量将 N 点的有限长序列 xn分段求周期图再平均,仿真结果：其信号时域图如图 2.2,功率谱估量如图 2.6 所示；4、Welch 法挑选了适当的窗函数wn,在分段时使各段之间有重叠；仿真结果：其信名师归纳总结 号时域图如图2.2,功率谱

25、估量如图2.7 所示；（分段长度为L=64,加汉明窗）第 11 页,共 20 页- - - - - - -精选学习资料 - - - - - - - - - 5-25Burg Power Spectral Density EstimatedB agnitude M0z /H dB ncy er/freque ow P-30-35-5-10-40-15-45-20-50-55-25-4-3-2-101234x 105-6005101520253035404550-30-5-65 FrequencyHz 4Frequency kHz图 2.8 AR 模型谱估量图 2.7 Welch 法功率谱估量5、

26、AR 模型谱估量在MATLAB 中,主要由基于 AR 模型的功率谱估量 pyulear函数、 pburg函数、pcov函数、pmov函数和 pmem函数,此处采纳基于 burg算法求出的 AR功率谱曲线,如图 2.8 所示；四、小结通过上面对单音正弦加噪信号的不同功率谱估量方法仿真分析,我们可以看出,相关函数 BT法和周期图法得到的谱线图离散性大,曲线粗糙,方差较大；平均周期图法的收敛性较好, 曲线比较平滑, 估量的结果方差较小, 但是其功率谱主瓣较 BT法更宽,辨论率削减了；采纳Welch法可以看到其功率谱曲线平滑,由于用了窗函数 汉明窗 极大地压低了旁瓣, 但是从图中可以看出, 其主瓣宽度

27、加宽了,功率谱波峰变宽了,降低了辨论率；采纳 谱曲线图主瓣窄,频谱方差小,辨论率高；AR模型的 burg算法求得的功率留意：在分析数据后面填充 0,是分析点数增多,可以对信号的频谱观测变得更细,但是根本上不能提高频率辨论率,由于在数据后面加 0,其傅里叶变换结果没有受到影响,因此也不能提高辨论率；五、参考文献1 王凤瑛 ,张丽丽 .功率谱估量及其 MATLAB 仿真 J.微运算机信息 ,2006,2231. 2 王福杰,潘宏侠 .MA TLAB 中几种功率谱估量函数的比较分析与挑选 J. 电子产品牢靠性与环境试验 ,2022,276. 名师归纳总结 3 丁玉美等 . 数字信号处理 . 西安：西

28、安电子科技高校出版社,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 = cos2*pi*f*t; xn = awgnx,10; %加10dB噪声 %= 各种谱估量方法 = Typ

29、e = 5; % 取不同值采纳不同方法谱估量 switch Type case 1 % 相关函数 BT 法 nfft = 512; % FFT点数 cxn = xcorrxn,unbiased; %运算序列的自相关函数 CXk = fftcxn,nfft; f1 = 0:nfft-1/nfft; Px_dB = 10*log10CXk.*conjCXk/nfft+eps; Px_dB = fftshiftPx_dB; figure1 subplot211 plott,x,r,xlabel t ,ylabel x , title 时域信号 ,grid on subplot212 plott,xn

30、,xlabel t ,ylabel xn , title 加噪信号 ,grid on figure2 plotfs*f1-fs/2,Px_dB,b, xlabel FrequencyHz ,ylabelMagnitude dB,grid on case 2 %周期图法 nfft = 512; f1 = 0:nfft-1/nfft; Px = fftxn,nfft; Px_dB = 10*log10absPx.2/nfft+eps; Px_dB = fftshiftPx_dB; plotfs*f1-fs/2,Px_dB,b, xlabel FrequencyHz ,ylabelMagnitude

31、 dB,grid on case 3 % 修正的周期图法 L = 32; % 数据分段长度 书中要求 名师归纳总结 - - - - - - -第 13 页,共 20 页精选学习资料 - - - - - - - - - M = N/L; % 分段数 n = 1 : L; L_nfft = 512; %每一段的 FFT点数sum_P = zeros1,L_nfft; % 对每一段数据进行 FFT,并求和 for k = 1 : M x = xnn:k-1*L; sum_P = sum_P + absfftx,L_nfft.2/L_nfft; %对每段求 FFT end f1 = 0:L_nfft-

32、1/L_nfft; Px_dB = 10*log10sum_P+eps; Px_dB = fftshiftPx_dB; plotfs*f1-fs/2,Px_dB,b,xlabel FrequencyHz , ylabelMagnitude dB,grid on case 4 % Welch法 %对数据进行补零 M = ceillog2lengthxn; N = 2M; % 求xn的长度对应的 2的最低幂次 m %如xn的长度不是 2的幂,补零到 2的整数幂 if lengthxn=0.5L, 即最多重叠 50% D = L - K; % 重叠部分长度 D Mmax = floorN-L/K+1

33、; %分段数运算 M = Mmax ; % 分段数windowb = hammingL; % 汉 明 窗 系 数blackmanL=blackman窗 ,boxcarL= 矩形窗 % beta =20; % windowb = kaiserL,beta; % 加窗提高辨论率 U = sumwindowb.2/L; n = 1:L; L_nfft = 512; sumP = zeros1,L_nfft; for i = 1:M X = xnn+i-1*K; % 对每段数据进行提取 X = X.*windowb; % 数据加窗 sumP = sumP+absfftX,L_nfft.2/L_nfft

34、 ; % 对每段求 FFT end sumP = sumP/M*U; Px_dB = 10*log10sumP+eps; Px_dB = fftshiftPx_dB; 名师归纳总结 - - - - - - -第 14 页,共 20 页精选学习资料 - - - - - - - - - f1 = 0:L_nfft-1/L_nfft; plotfs*f1-fs/2,Px_dB,b, xlabel FrequencyHz ,ylabelMagnitude dB,grid on case 5 % AR模型谱估量 order = 50; range = half; nfft = 512; figure 基

35、于 burg算法求出的 AR 功率谱曲线pburgxn,order,nfft,fs,range; end名师归纳总结 - - - - - - -第 15 页,共 20 页精选学习资料 - - - - - - - - - 试验三 卡尔曼滤波器的设计一、卡尔曼滤波器介绍卡尔曼滤波是用状态空间法描述系统的,由状态方程和测量方程所组成；卡尔曼滤波用前一个状态的估量值和最近一个观测数据来估量状态变量的当前值,并以状态变量的估量值的形式给出；卡尔曼滤波的状态方程和测量方程分别为式3-1、式3-2所示；x k 1 A k x k w k 3 1y k C k x k v k 3 2 其中, k 表示时间,输入信号 w 是一白噪声,输出信号的观测噪声 kv 也是一个白噪声, A 表示状态变量之间的增益矩阵,可随时间变化,A 表示第 k 次迭代的取值, C 表示状态变量与输出信号之间的增益矩阵,模型如图 3.1 所示 k 用k1代替 ；可随时间变化, 其信号wk1+xkZ1kxC kvkyk+1A k1图 3.1 卡尔曼滤波器的信号模型二、卡尔曼滤波的递推算法卡尔曼递推公式如下式 3-3a、3-3b、3-3c及3-3d所示；名师归纳总结 假设初始条件xkA kxk1Hkyk1,CkA kxk1303 ax 0,P 0varx 0,第 16 页,共