基于EMD的非平稳信号去噪方法研究

1、中北大学课 程 设 计 说 明 书学生姓名： 李晓聪 学 号： 1005064125 学 院： 信息与通信工程学院 专 业： 电子信息工程 题 目： 信息处理综合实践: 基于EMD的非平稳信号去噪方法研究 指导教师： 赵英亮 王浩全 职称: 副教授 2013 年 6 月 27 日课 程 设 计 任 务 书1设计目的：1、掌握USB总线或PCI总线的基本结构，了解基于USB总线或PCI总线A/D卡的通用结构。写出关于基于USB总线或PCI总线A/D卡的报告。2、通过A/D卡，利用高级语言编写信号的采集、存储和显示程序。 3、要求同学在学习EMD基本理论的基础上，对一维信号进行各种降噪方法的研究，

2、实现信号去噪。 2设计内容和要求（包括原始数据、技术参数、条件、设计要求等）：1、查阅资料，掌握USB总线或PCI总线的基本结构，了解基于USB总线或PCI总线A/D卡的通用结构。2、通过A/D卡，利用高级语言编写信号的采集、存储和显示程序。3、查阅资料，掌握EMD基本原理。4、利用编程语言实现EMD滤波方法设计。5、对原始信号进行不同噪声叠加，并利用EMD去噪方法对其实现处理，将处理后结果与原信号进行比较。对其结果进行理论及结果分析，比较其优缺点，及对各种噪声的适用范围。6、要求每位同学各自应用不同的方法来实现。3设计工作任务及工作量的要求包括课程设计计算说明书(论文)、图纸、实物样品等：课

3、程设计说明书一份仿真结果课 程 设 计 任 务 书4主要参考文献：1 徐明远,刘增力.MATLAB仿真在信号处理中的应用.西安:西安电子科技大学出版社,2007.11.2 王玉顺.MATLAB实践教程. 西安:西安电子科技大学出版社,2012.8.3 唐向宏，岳恒立，郑雪峰.MATLAB及在电子信息类课程中的应用. 北京:电子工业出版社,2009.6.4 李迎，王帮峰，孙亚飞，施益峰.经验模式分解与Hilbert谱的分析及应用. 南京:南京航空航天大学智能材料与结构航空科技重点实验室,21 00 16.5 李建伟，曾杰，许宝杰.基于EMD的非平稳信号分析方法研究及其在机电系统信号分析中的应用.

4、5设计成果形式及要求：课程设计说明书仿真结果6工作计划及进度：2013年6月7日 6月10日：查资料；6月11日 6月13日：在指导教师指导下设计方案；6月14日 6月25日：完成设计内容；6月26日 6月27日：撰写课程设计说明书； 6月28日：答辩系主任审查意见： 签字： 年 月 日目录1.设计目的12.总体设计方案分析12.1 PCI总线基本结构12.2 基于PCI总线的实时数据采集系统32.3基于PCI总线AD卡的通用结构42.3.1 A/D 转化过程52.3.2 存储格式52.4 基于EMD的非平稳信号分析方法62.4.1 IMF分量需满足的条件62.4.2 EMD的主要分解过程62

5、.4.3 对各个IMF分量的后续分析73.主要算法及程序83.1 主要算法步骤83.2 Main程序84.算法结果及比较分析125.设计评述166.参考文献171.设计目的(1).掌握USB总线或PCI总线的基本结构，了解基于USB总线或PCI总线A/D卡的通用结构。写出关于基于USB总线或PCI总线A/D卡的报告。(2).通过A/D卡，利用高级语言编写信号的采集、存储和显示程序。 (3).要求在学习EMD基本理论的基础上，对一维信号进行各种降噪方法的研究，实现信号去噪。 2.总体设计方案分析我们要对一个原始信号进行不同噪声的叠加，并利用基于EMD去噪方法对其实现处理。利用EMD对一维信号进行

6、去噪包含了三种EMD去噪方法：1、直接EMD方法去噪；2、EMD结合小波阈值去噪；3、EMD结合互相关函数的方法。直接EMD方法去噪是直接抽取部分与原始信号频率相似的低频IMF分量进行信号重构，抛弃高频IMF分量与无效分量。这种方法效果不好，而且对的能处理的信号有限制。EMD结合小波阈值去噪是如果噪声幅度过大，就先进行EMD分解，再进行小波阈值去噪，因为方案原始信号较简单，频率都为低频分量，所以可以直接对IMF高频分量进行小波阈值去噪，再进行信号重构。它和直接EMD方法去噪的区别在于用什么方法区别高频IMF分量和高频噪声。EMD结合互相关函数的方法也是一种较粗糙的去噪方法。2.1 PCI总线基

7、本结构PCI总线作为处理器系统的局部总线，主要目的是为了连接外部设备，而不是作为处理器的系统总线连接Cache和主存储器。但是PCI总线、系统总线和处理器体系结构之间依然存在着紧密的联系。PCI总线作为系统总线的延伸，其设计考虑了许多与处理器相关的内容，如处理器的Cache共享一致性和数据完整性，以及如何与处理器进行数据交换等一系列内容。其中Cache共享一致性和数据完整性是现代处理器局部总线的设计的重点和难点。独立地研究PCI总线并不可取，因为PCI总线仅是处理器系统的一个组成部分。PCI总线作为处理器系统的局部总线，是处理器系统的一个组成部件，讲述PCI总线的组成结构不能离开处理器系统这个

8、大环境。在一个处理器系统中，与PCI总线相关的模块如图11所示。图（1）基于PCI总线的处理器系统如图（1）所示在一个处理器系统中，与PCI总线相关的模块包括，HOST主桥、PCI总线、PCI桥和PCI设备。PCI总线由HOST主桥和PCI桥推出，HOST主桥与主存储器控制器在同一级总线上，PCI设备可以方便地通过HOST主桥访问主存储器，即进行DMA操作。在处理器系统中，含有PCI总线和PCI总线树这两个概念。这两个概念并不相同，在一颗PCI总线树中可能具有多条PCI总线，而具有血缘关系的PCI总线组成一颗PCI总线树。如在图11所示的处理器系统中，PCI总线x树具有两条PCI总线，分别为P

9、CI总线x0和PCI总线x1。而PCI总线y树中仅有一条PCI总线。在PCI总线中有三类设备，PCI主设备、PCI从设备和桥设备。其中PCI从设备只能被动地接收来自HOST主桥，或者其他PCI设备的读写请求；而PCI主设备可以通过总线仲裁获得PCI总线的使用权，主动地向其他PCI设备或者主存储器发起存储器读写请求。而桥设备的主要作用是管理下游的PCI总线，并转发上下游总线之间的总线事务。一个PCI设备可以即是主设备也是从设备，但是在同一个时刻，这个PCI设备或者为主设备或者为从设备。PCI总线规范将PCI主从设备统称为PCI Agent设备。在处理器系统中常见的PCI网卡、显卡、声卡等设备都属

10、于PCI Agent设备。在PCI总线中，HOST主桥是一个特殊的PCI设备，该设备可以获取PCI总线的控制权访问PCI设备，也可以被PCI设备访问。但是HOST主桥并不是PCI设备。PCI规范也没有规定如何设计HOST主桥。在PCI总线中，还有一类特殊的设备，即桥设备。桥设备包括PCI桥、PCI-to-(E)ISA桥和PCI-to-Cardbus桥。本篇重点介绍PCI桥，而不关心其他桥设备的实现原理。PCI桥的存在使PCI总线极具扩展性，处理器系统可以使用PCI桥进一步扩展PCI总线。PCI桥的出现使得采用PCI总线进行大规模系统互连成为可能。但是在目前已经实现的大规模处理器系统中，并没有使

11、用PCI总线进行处理器系统与处理器系统之间的大规模互连。因为PCI总线是一个以HOST主桥为根的树型结构，使用主从架构，因而不易实现多处理器系统间的对等互连。2.2 基于PCI总线的实时数据采集系统PCI-20612是一种基于32位PCI总线的数据采集卡。它采用了大规模门阵列设计技术、表面贴装工艺等先进手段制作而成，广泛应用于实验室、生产现场或野外对冲击、爆炸等原因引起的压力、电磁波、温度等各种动态过程进行实时数据采集。图（2）PCI20612数据采集卡工作原理图信号连接器有通道信号，PCI总线接口，外触发，多卡扩展同步接口等4种，以下分别介绍。通道信号接口位于采集卡的尾部，信号用同轴电缆接入

12、。一张采集卡最多只有4个通道，也可能少于4个通道，视用户的需求而定。PCI总线接口是采集卡与上位PC机的数据和控制信号通信接口。其采用32位PCI总线，全兼容PCI Specification Version 2.2。外触发接口用于接外触发信号用。多卡扩展同步接口用于多卡扩展。其传输的信号包括主卡的同步时钟源、触发以及状态信号。一般用16针带缆连接。2.3基于PCI总线AD卡的通用结构PCI-20612是一种基于32位PCI总线的数据采集卡。它采用了大规模门阵列设计技术、表面贴装工艺等先进手段制作而成，广泛应用于实验室、生产现场或野外对冲击、爆炸等原因引起的压力、电磁波、温度等各种动态过程进行

13、实时数据采集。这里介绍PCI-20612数据采集卡的工作原理，其中包括了寄存器地址分配表，信号连接器，模拟输入部分，A/D转换过程，数字I/O接口及多卡同步扩展接口总线等等。2.3.1 A/D 转化过程当进行A/D转换时，根据所测的信号的特性，程控A/D转换设置，如：通道信号输入带宽、偏置和增益调节（自校准）等。图（3） A/D转换过程图如图（3）所示，当启动采集后，A/D转换后的数据经锁存，然后保存于每通道独立的卡上存储器中。卡上存储器相当于环形缓冲，如果A/D转换的数据样点数超过了卡上存储器的最大容量，新数据会覆盖旧数据。这个过程是周而复始的，只有当触发条件满足后，门阵列开始计数，计数达到

14、指定值（该值由采集长度决定）后，采集结束，卡上存储器保存了满足用户需要的采集数据。上位机通过PCI接口门阵列经由门阵列控制核心取得卡上存储器样点数据。上过程每通道是独立进行的。每通道卡上存储器最多可存储8M样点，且此参数可由DIP开关设置。2.3.2 存储格式图（4） 每通道卡上存储器的数据存储格式每通道卡上存储器的数据存储格式如图4-6所示。样点起点并不是从物理存储单元的0地址开始，而是由触发起点和延时长度决定。2.4 基于EMD的非平稳信号分析方法经验模态分解(Empirical Mode Decomposition, 简称EMD)方法是由美国NASA的黄锷博士提出的一种信号分析方法。它依

15、据数据自身的时间尺度特征来进行信号分解, 无须预先设定任何基函数。2.4.1 IMF分量需满足的条件1）其极值点和过零点数目相等或最多相差一个。2）连接其局部极大值和局部极小值所形成的两条包络线的均值在任一点处为零或近似为零。2.4.2 EMD的主要分解过程分解过程基于以下假设：（1）信号最少有一个极大值和一个极小值；（2）时域特性由极值间隔决定；（3）如果数据序列完全缺乏极值但是仅包含拐点，那么它也可通过求导一次或多次来揭示极值点，而最终结果可以由这些成分求积分来获得。具体方法是由一个“筛选”过程完成的：（1）首先找出所有的极大值点并将其用三次样条函数拟合成原数据序列的上包络线：以及所有的极

16、小值点并将其用三次样条函数拟合成原数据序列的下包络线；（2）计算上下包络线的均值，记为；将原数据序列减去该均值即可得到一个去掉低频的新数据序列： （3）因为一般仍不是一个IMF分量序列，为此需要对它重复进行上述处理过程。重复进行上述处理过程七次，直到符合IMF的定义要求，所得到的均值趋于零为止，这样就得到了第1个IMF分量，它代表信号中最高频率的分量： （4）将从中分离出来，即得到一个去掉高频分量的差值信号，即有 将作为原始数据，重复步骤（1）、（2）（3），得到第二个IMF分量，重复n次，得到n个IMF分量。这样就有： 1）找出待分析信号的全部极大值和极小值点，利用三次样条函数分别把它们拟和

17、为该信号的上下包络线，计算出两包络线的均值，进而求出待分析信号和均值的差值h；2）若h不满足本征模函数的要求，则对其重复上述过程k次，使得新的h满足本征模函数的条件；若h满足本征模函数的要求，则令其为原信号的第一个本征模函数，并求出原信号与该本征模函数的差值r；3）将r作为待分解信号，重复以上过程，直至所剩余的r不可分解或研究意义不大为止。经过EMD分解后，信号变为一系列IMF分量和一个趋势项之和，各IMF分量按频率段由高到低依次分解得到。每一IMF分量所包含的频率成分随信号本身变化而变化，因此EMD方法是一种自适应的信号分析方法，非常适用于非线性和非平稳过程，具有很高的信噪比。2.4.3 对

18、各个IMF分量的后续分析EMD作为一种信号预分析方法，将一个复杂的时序信号分解为一系列窄带信号。要获得原信号较为完整的时频信息，就要对各窄带信号进行进一步的分析。小波分析理论对非平稳信号的分析提供了有效的途径，小波变换的基本思想与傅里叶变换是一致的，它也是用一族函数来表示信号或函数，这一组函数称为小波函数系，它是一种窗口面积恒定，窗口形状可变的时频局域化分析方法。小波去噪首先对实际信号进行小波分解，选择小波并确定分解层次为N，则噪声部分通常包含在高频中；然后对小波分解的高频系数进行门限阈值量化处理；最后根据小波分解的第N层低频系数和经过量化后的1N层高频系数进行小波重构，达到消除噪声的目的，对

19、信号降噪实质上是抑制信号的噪声，在实际信号中恢复真实信号的过程。3.主要算法及程序3.1 主要算法步骤(1)对信号x(t)进行进行EMD分解得到各模态分量imf1imfk；(2)对高频分量imf1imfh,选择阈值确定准则进行各个分量上的阈值估计；(3)选择一种合适的阈值函数，结合高频分-量imf1imfh上的阈值估计值进行去噪计算，得到去噪后各个分量imf1imfh；(4)以去噪后的分量imf1imfh和原来分解的imfi+himfk来重构信号，该信号即为去噪后的信号x(t); 式中x(t)为去噪后的信号，ci(t)为阈值处理后的imf1imfh分量和原来分解的imfi+himfk分量的和，

20、rn(t)为残余分量。3.2 Main程序clc;close all;clear all;%产生原始信号，噪声信号，加噪信号；N=2048;t=linspace(1,2,N);x1=5*sin(2*pi*10*t); x2=5*sin(2*pi*35*t);x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);figure(1); subplot(3,1,1); plot(x1); title(10hz信号);subplot(3,1,2); plot(x2); title(35hz信号);subplot(3,1,3); plot(x); title(原始信号);xr=rand

21、(1,2048);figure(2);plot(xr); title(噪声信号);y=x+xr;figure(3);plot(y); title(加噪信号);%执行emd算法，对加噪信号进行emd分解，产生emd分量；imf=emd(y); imf1=imf(1,:);figure(4);plot(imf1);title(IMF1); imf2=imf(2,:);figure(5);plot(imf2);title(IMF2); imf3=imf(3,:);figure(6);plot(imf3);title(IMF3); imf4=imf(4,:);figure(7);plot(imf4);

22、title(IMF4); imf5=imf(5,:);figure(8);plot(imf5);title(IMF5); imf6=imf(6,:);figure(9);plot(imf6);title(IMF6); imf7=imf(7,:);figure(10);plot(imf7);title(IMF7); imf8=imf(8,:);figure(11);plot(imf8);title(残余项); %提取高频imf分量，进行小波阈值去噪；C,S=wavedec2(imf1,3,db9); dcoef=C(prod(S(1,:)+1:end); ind=find(abs(dcoef)=

23、8)+prod(S(1,:); C(ind)=sign(C(ind).*sqrt(C(ind).2-2*28); X1=waverec2(C,S,db9); C,S=wavedec2(imf2,3,db9); dcoef=C(prod(S(1,:)+1:end); ind=find(abs(dcoef)=8)+prod(S(1,:); C(ind)=sign(C(ind).* sqrt( C(ind).2-2*28);X2 = waverec2(C, S, db9); %分析结果；figure(12);subplot(2,1,1); plot(x); title(原始信号);subplot(2,1,2); plot(X1+X2+imf3+imf4+imf5+imf6+imf7); title(分析结果);axis(0 2500 -10 10);4.算法结果及比较分析 图（5）原始信号 图（6）噪声信号图（7） 加噪信号图（8） IMF1分量 图（9） IMF2分量 图（10） IMF3分量图（11） IMF4分量图（12） IMF5分量图（13）