共模噪声是电子电路中的一个基本概念，首先来介绍几个名词。
1、共模与差模电流
根据基尔霍夫定律，流入一个节点的电流总量等于流出这个节点的电流总量。很多电路图
上只画出源发射到负载的信号线路，看不出流回信号源的回路；实际上，电流是经过地线流回
到源的，地线作为电流返回到源的低阻抗路径。 
差模电流：在一对差分信号线上，大小相同，方向相反的一对信号，一般是电路中的工作
电流，对于信号线就是信号线与信号地线之间流动的电流。
 共模电流：在一对差分信号线上，大小相同，方向相同的一对信号(或噪音)。在电路中，
一般对地噪音一般都是以共模电流的方式传输的，所以又称为共模噪声。 
2、共模电流的产生
共模电流既可以由设备外部的因素在电缆上产生，也可以由设备自身的因素在电缆上产生；在高频开关电源设计中，广泛利用共模电感产生共模电流，来供给初级变压器。
外部因素：外部导致共模电流产生的主要因素是地线噪声电压和电磁场，次要原因则包括
电磁波、雷电等空间干扰在电缆上感应出的共模电压感应出的共模电流。 
地线噪声电压：电流总是沿着阻抗最低的线路流动，理想情况下地线是电流返回源的阻抗
最低路径，然而实际的情况是，地线作为导线依然存在着阻抗，电流并不总是一定按照地线返
回。特别是对于某些频率高的电流，因为杂散电容的存在，电流通过这一耦合通道流回源，而
与地线上的电压形成了一个共模电流的驱动源。 
电磁场：两根导线上感应出的共模电压是相同的，但由于两个导线对地的阻抗不同，使得
产生的共模电流不同，导致差模电压产生，从而对原有的差模信号造成干扰。 
设备产生较强电磁辐射的一个重要原因是共模电流，而共模电流来自于设备内部产生的共
模电压，来源主要是以下几种： 
工作的信号电压：因为电路周围导体存在杂散电容，电压通过杂散电容产生电流，大小取
决于电压的幅度和回路阻抗的大小； 
线路板的地线噪声：线路板上的地线作为信号回流线，总是会存在电压噪声，当有线缆连
接到电路板上时就产生了共模电流； 
设备内部的电磁耦合：当设备内部的电路产生较强的电磁场，处于这个电磁场中的导体都
会感应出电压，而这个电压就会产生电流，从而对线路造成影响； 
开关器件上的脉冲电压：主要是开关电源、逆变器等部件中的IGBT等器件，这些期间上
具有较高的脉冲电压，其中含有大量高频成分，通过空间杂散电容在线路中形成共模电流。 
3、共模电流的抑制 
为了减小电路的电磁干扰，需要对共模电流进行抑制，常见的抑制方法可分为以下几类：
减小地线阻抗：地线上更小的阻抗意味着更低的噪声电压，高频下线路阻抗与长度的
关系更密切，因此减小地线阻抗需要缩短地线长度。另外还可以采用多层线路板，其中一层专
门作为地线，也是为了便于接地； 
更好的接地：对于较低频线路，使用单点接地的方式，减小电路的面积，从而减少电
路所包围的磁通量，减少电磁感应；对于高速线路则采用多点接地，所有电路就近连接到地，
保持地线最短以减小地线上的阻抗； 
将两个电路隔离：有时因为安全的目地，设备必须接地，单点接地不可能。这时可以
使用隔离变压器或光耦隔离器来连接两台设备，切断两个设备之间的电气连接，差模信号通过
磁场或光传送，同时切断地线产生的共模干扰； 
强干扰电路远离 I/O 端口：高速数字脉冲电路、时钟电路、振荡器电路等在工作时会
产生较强的干扰，这些电路要尽量远离I/O接口电路，防止干扰耦合到I/O电缆上。 
使用旁路电容：有时为了抑制线路中的高频共模成分造成的辐射，会在线路中增加旁
路电容作为高频噪声流入地的通道。但由于加入电容，使得共模回路的阻抗减小，会增大线路
中感应出的共模电流，进而增大负载上的共模电压； 
使用共模扼流器：当无法对地线等部分进行修改，或者旁路电容使线路的共模抗干扰
能力受影响时，可以在线路中串联共模扼流器件。这样做的目地是增大共模回路的阻抗，使得
共模电流被扼流器所消耗，从而抑制线路中的共模噪声，这在各种电子设备的电源插头线上很常见。 
 共模扼流器的原理：
以共模磁环为例。在一铁氧体磁环上绕上同向的一对线圈，当交变电流通过时，因为电磁
感应而在线圈中产生磁通量。对于差模信号，产生的磁通量大小相同方向相反，两者相互抵消，
因而磁环产生的差模阻抗非常小；而共模信号产生的磁通量大小方向均相同，两者相互叠加从
而使磁环产生了较大的共模阻抗。这一特性使得共模扼流器对于差模信号的影响较小而对共模
噪声具有很好的滤波性能。
针对一些数据传输接口，例如 USB2.0及以上的接口，其传输速率很高，此时信号和噪音
的频率差别不大。在这种情况下若使用磁珠滤波，虽然滤除了噪音，但是也削弱了有用信号，
而使用共模扼流器则能够对差模信号与共模噪声进行分别，滤除噪音保留信号。

小波分析的重要应用之一就是用于信号消噪，首先简要地说明一下小波变换实现信号消噪的基本原理。

一、基本原理

含噪的一维信号模型可以表示如下



式中s(k)为含噪信号，f(k)为有用信号，e(k)为噪声信号。这里假设e(k)是一个高斯白噪声，通常表现为高频信号，而实际工程中f(k)通常为低频信号或者是一些比较平稳的信号。

因此我们可按如下方法进行消噪处理:

首先对信号进行小波分解，由于噪声信号多包含在具有较高频率的细节中，从而可以利用门限、阈值等形式对分解所得的小波系数进行处理，然后对信号进行小波重构即可达到对信号消噪的目的。

对信号消噪实质上是抑制信号中的无用部分，增强信号中有用部分的过程。

一般地，一维信号消噪的过程可分为如下3个步骤。

步骤1:一维信号的小波分解。选择一个小波并确定分解的层次，然后进行分解计算。

步骤2:小波分解高频系数的阈值量化。对各个分解尺度下的高频系数选择一个阈值进行软阈值量化处理。

步骤3:一维小波重构。根据小波分解的最底层低频系数和各层高频系数进行一维小波重构。

这3个步骤中，最关键的是如何选择阈值以及进行阈值量化。在某种程度上，它关系到信号消噪的质量。

二、噪声在小波分解下的特性

总体上，对于一维离散信号来说，其高频部分所影响的是小波分解的第一层细节，其低频部分所影响的是小波分解的最深层和低频层。如果对一个仅由白噪声所组成的信号进行分析，则可得出这样的结论:高频系数的幅值随着分解层次的增加而迅速地衰减，且其方差也有同样的变化趋势。

用C(j,k)表示噪声经过小波分解后的系数，其中j表示尺度，k表示时间。

下面将噪声看做普通信号，分析它的相关性、频谱和频率分布这3个主要特征。

(1)如果所分析的信号s是一个平稳的、零均值的白噪声，那么它的小波分解系数是相互独立的。

(2)如果信号s是个高斯型噪声，那么其小波分解系数是互不相关的，且服从高斯分布。

(3)如果信号s是一个平稳、有色、零均值的高斯型噪声序列，那么它的小波分解系数也是高斯序列，并且对每一个分解尺度j，其相应的系数也是一个平稳、有色的序列。如何选择对分解系数具有相关性的小波是一个很困难的问题，在目前也没有得到很好的解决。进一步需要指出的是，即使存在这样的一个小波，但它对噪声的解相关性还取决于噪声的有色性。为了利用小波计算噪声的解相关性，必须事先知道噪声本身的特性。

(4)如果信号s是一个固定的、零均值ARMA模型，那么对每一个小波分解尺度j，C(j,k)(k∈z)也是固定的、零均值ARMA模型，且其特性取决于尺度j。

(5)如果信号s是一个一般的噪声

①若它的相关函数已知，则可计算系数序列C(j,k)和C(j,k’);

②若它的相关函数谱已知，则可计算C(j,k)(k∈z)的谱及尺度j和j’的交叉谱。

三、应用一维小波分析进行信号消噪处理

小波分析工具箱中用于信号消噪的一维小波函数是wden.m的wdencmp.m。

小波分析进行消噪处理一般有下述3种方法。

(1)默认阈值消噪处理。该方法利用函数ddencmp生成信号的默认阈值，然后利用函数wdencmp进行消噪处理。

(2)给定阈值消噪处理。在实际的消噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh。

(3)强制消噪处理。该方法是将小波分解结构中的高频系数全部置为0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且消噪后的信号比较平滑，但是容易丢失信号中的有用成分。

[例1]利用小波分析对含噪正弦波进行消噪。

%% 利用小波分析对含噪正弦波进行降噪
clc;clear;
% 1、生成正弦信号
N=1000;
t=1:N;
x=sin(0.03*t);
% 2、加噪声
load noissin;% 这是一个含噪正弦波
ns=noissin;
% 3、显示波形
subplot(3,1,1);
title('含噪正弦波消噪结果');
plot(t,x);
xlabel('样本序号 n');
ylabel('(原始信号)幅值 A');
subplot(3,1,2);
plot(ns);
xlabel('样本序号n');
ylabel('(含噪信号)幅值A');
% 4、小波消噪
xd=wden(ns,'minimaxi','s','one',5,'db3');
% Wden函数：一维信号的小波消噪处理
% [xd，cxd，lxd]=wden（x，tptr，sorh，scal，n，‘wname’）；返回经过小波消噪处理后的信号xd及其小波分解结构。
% 1、x = ns为含噪信号
% 2、tptr = rigrsure：stein无偏估计； heursure：启发式阈值； sqtwolog：固定式阈值； minimaxi：极大极小值阈值
% 3、sorh = s：为软阈值； h： 为硬阈值；
% 4、输入参数scal规定了阈值处理随噪声水平的变化：
%   Scal=one：不随噪声水平变化。
%    Scal=sln：根据第一层小波分解的噪声水平估计进行调整。
%   Scal=mln：根据每一层小波分解的噪声水平估计进行调整。
% 5、N 是小波变换的级别。
% 6、WNAME 是指定为字符串的小波。

subplot(3,1,3);
plot(xd);
xlabel('样本序号n');
ylabel('(消噪信号)幅值A'); 



从上图可以看出:消噪后的信号大体上恢复了原始信号的形状，并明显地除去了噪声所引起的干扰。但是，恢复后的信号和原始信号相比，有明显的改变。这主要是因为在进行消噪处理的过程中所用的分析小波和细节系数阙值不恰当所致。

下面我们再通过一个例程对此进行说明。

[例2]在电网电压值监测过程中，由于监测设备出现了一点故障，致使所采集到的信号受到噪声的污染。现在利用小波分析对污染信号进行消噪处理以恢复原始信号。

%% 分别使用3中阈值来消噪，并分析其适用场合：默认阈值消噪、给定阈值消噪、强制消噪
% 零、信号加载与处理(分解、系数提取、重构)
% 1、装载采集的信号leleccum.mat
load leleccum;
% 将信号中第2000到第3450个采样点赋给s
indx=2000:3450;
s=leleccum(indx);
% 2、画出原始信号
subplot(2,2,1);
plot(s);
title('原始含噪信号');
xlabel('样本序号n');
ylabel('幅值A');

% 3、对原始信号进行分解
% 用db1小波对原始信号进行3层分解并提取系数
[c,l]=wavedec(s,3,'db1');   %这里c的长度为1453，比s（1451）多，l的长度为5

% function [c,l] = wavedec(x,n,IN3,IN4)：多层一维小波分解。
% WAVEDEC 使用一个特定的小波“wname”或一组特定的小波分解滤波器进行多层一维小波分析。
% [C,L] = WAVEDEC(X,N，'wname') 使用'wname'返回信号X在N级的小波分解。
% 输出分解结构包含小波 the wavelet decomposition vector C和 the bookkeeping vector L
% N 必须是一个严格的正整数。

% 4、对[c,l]提取近似系数
a3=appcoef(c,l,'db1',3);
% APPCOEF: 提取一维小波变换近似系数。
% A = APPCOEF(C,L,'wname',N) 使用小波分解结构计算N级的近似系数[C,L]
% Level N must be an integer such that 0 <= N <= length(L)-2. 
% A = APPCOEF(C,L,'wname') extracts（提取） the approximation coefficients（近似系数） at the last level length(L)-2.

% 5、对[c,l]提取细节系数
% 5.1、提取3级细节系数
d3=detcoef(c,l,3);
%  DETCOEF提取一维细节系数。
%  D = DETCOEF(C,L,N)从小波分解结构[C,L]中提取出N级的细节系数(the detail coefficients)。
%  Level N must be an integer such that 1 <= N <= NMAX where NMAX = length(L)-2.
%  D = DETCOEF(C,L)提取最后一级NMAX的细节系数。
%  If N is a vector of integers such that 1 <= N(j) <= NMAX:

% 5.2、提取2级细节系数
d2=detcoef(c,l,2);
% 5.3、提取1级细节系数
d1=detcoef(c,l,1);

% 一、强制消噪
%对信号进行强制性消噪处理并图示结果
dd3=zeros(1,length(d3)); % length(d3) = 182，这句话执行完，生成了182个0
dd2=zeros(1,length(d2)); % length(d2) = 363
dd1=zeros(1,length(dl)); % length(dl) = 726
% 将近似系数a3与其他细节系数（均为0）合并， length(a3) = 182, length(cl) = length(c) = 1453
c1=[a3 dd3 dd2 dd1]; % 由于dd3 dd2 dd1 均为0，那么c1有用信息仅来自使用db1进行的3级分解的近似系数a3
s1=waverec(c1,l,'db1');% 用c1来重构c，命名为s1（重构后的信号s1，与原信号c对应）

% WAVEREC：多层一维小波重建。与WAVEDEC相对应
% WAVEREC使用特定的小波('wname'，参见WFILTERS)或特定的重建滤波器(Lo_R和Hi_R)来执行多层一维小波重构。
% X = WAVEREC(C,L，'wname')根据多级小波分解结构对信号X进行重构

subplot(2,2,2);
plot(s1);
title('强制消噪后的信号');
xlabel('样本序号n');
ylabel('幅值A');


% 二、用默认阈值对信号进行消噪处理并图示结果
% 1、用ddencmp函数获得信号的默认阈值
[thr,sorh,keepapp]=ddencmp('den','wv',s); 
%使用ddencmp函数的去噪功能对原始信号s进行小波分解，并得到阈值thr，阈值方式，是否允许近似保持系数

%  DDENCMP:去噪或压缩的默认值
%  [THR,SORH,KEEPAPP,CRIT] = DDENCMP(IN1,IN2,X)：返回输入向量或矩阵X的默认去噪或压缩值，使用小波或小波包，可以是一维或二维信号。其中：
%  THR是阈值，SORH是函数（自动）选择的阈值方式：软阈值（SORH = s）或硬阈值（SORH = h），KEEPAPP允许保持近似系数（0、1），CRIT(仅用于小波包)是熵名(参见WENTROPY)。
%  IN1是den（去噪）或者cmp（压缩）
%  IN2是wv（小波分解）或者wp（小波包分解）。
 
s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);%利用ddencmp得到的阈值、阈值方式、近似保持系数对小波分解后的信号c进行去噪

% WDENCMP使用小波对信号或图像进行去噪或压缩处理。
% [XC,CXC,LXC,PERF0,PERFL2] = WDENCMP('gbl'，X，'wname'，N,THR,SORH,KEEPAPP)返回输入信号X 的去噪或压缩版本的XC，使用全局正阈值THR进行小波系数阈值化。
%  附加输出参数[CXC,LXC]是XC的小波分解结构，
% PERFL2和PERF0是以百分比表示的L^2恢复和压缩分数。
% PERFL2 = 100*(CXC的向量范数/ C的向量范数)^2其中[C,L]为X的小波分解结构。
% 小波分解在第N层执行，'wname'是一个包含小波名称的字符串。
% SORH ('s'或'h')用于软阈值或硬阈值(有关详细信息，请参阅WTHRESH)。
% 如果KEEPAPP = 1，逼近系数不能设置阈值，否则是可能的。


subplot(2,2,3);
plot(s2);
title('默认阈值消噪后的信号');
xlabel('样本序号n');
ylabel('幅值A');


% 三、用给定的软阈值进行消噪处理 
softd1=wthresh(d1,'s',1.465);% 对d1进行软阈值
% WTHRESH执行软或硬阈值。
% Y = WTHRESH(X,SORH,T)返回输入向量或矩阵X的软(如果SORH = 's')或硬(如果SORH = 'h') 阈值。T是阈值。
% Y = WTHRESH(X,'s',T) returns Y = SIGN(X).(|X|-T)+, soft  thresholding is shrinkage.
% Y = WTHRESH(X,'h',T) returns Y = X.1_(|X|>T), hard  thresholding is cruder.
softd2=wthresh(d2,'s',1.823); % 对d2进行软阈值
softd3=wthresh(d3,'s',2.768);% 对d3进行软阈值

% softdl softd2 softd3为执行阈值处理后的细节系数
c2=[a3 softd3 softd2 softd1];% 将近似系数与（经过阈值处理后的）细节系数合并
s3=waverec(c2,l,'db1'); % 使用c2和级数l来重构s3
subplot(2,2,4);
plot(s3);
title('给定软阈值消噪后的信号');
xlabel('样本序号n');
ylabel('幅值A');


 

在本例程中，分别利用前面所提到的了3种消噪方法进行处理。

从上图得到的结果来看：应用强制消噪处理后的信号较为光滑，但是它很有可能丢了信号中的一些有用成分；默认阈值消噪的效果较好，保留了有用信号的高频成分，同时对微小扰动去噪效果明显；给定软阈值消噪需要有一定的经验或者需要多次试验，也可作为备选。

四、消噪阈值选取规则

在Matlab的小波工具箱中，设置软或硬阈值的函数为wthresh.m。

该函数根据参数sorh的值计算分解系数的软阈值或硬阈值。其中，硬阀值对应于最简单的处理方法，而软阙值具有很好的数学特性，并且所得到的理论结果是可用的。

[例3]生成不同阀值下的信号。

% 生成线性信号
y=linspace(-1,1,100);
% 设置T阈值
thr=0.4;
% 计算软、硬阈值
ythard=wthresh(y,'h',thr);	% 对信号y进行硬阈值，阈值设为thr
ytsoft=wthresh(y,'s',thr);  % 对信号y进行软阈值，阈值设为thr
% 显示不同阈值后的信号
subplot(1,3,1);
plot(y);
title('原始信号');
xlabel('样本序号n');
ylabel('幅值A');
subplot(1,3,2);
plot(ythard);
title('硬阈值信号');
xlabel('样本序号n');
ylabel('幅值A');
subplot(1,3,3);
plot(ytsoft);
title('软阈值信号');
xlabel('样本序号n');
ylabel('幅值A');




根据基本的噪声模型，阈值的选取有以下4个规则，其中每一条规则对应于函数thselect中输入参数 tptr的一个选项。

(1)选项tptr='rigrsure'，是一种基于Stein无偏似然估计原理的自适应阙值选择。给定一个阈值t，得到它的似然估计，再将非似然最小化，就可得到所选的阈值。这是一种软件阈值估计器。

(2)选项 tptr='sqtwolog'，是一种固定的阈值形式，它所产生的阈值为 sprt(2*log(length(X)))。

(3)选项 tptr='heursure'，是前两种阈值的综合，所选择的是最优预测变量阈值。如果信噪比很小，而SURE估计有很大的噪声，此时就需要采用这种固定的阈值形式。

(4)选项tptr='minimaxi'，也是一种固定的阅值选择形式，它所产生的是一个最小均方差的极值，而不是无误差。

【例4-6】计算高斯白噪声信号N(0,1)在上述4条规则下的阈值。

%产生白噪声
y=randn(1,1000);
%生成不同的阙值
thr1=thselect(y,'rigrsure');
thr2=thselect(y,'sqtwolog');
thr3=thselect(y,'heursure');
thr4=thselect(y,'minimaxi');


计算结果为：



从计算结果看，rigrsure的结果会随高斯白噪声变化而变化，而其他三种不会变化；

在本例程中，信号y是一个标准的高斯白噪声，所以每一种方法都能粗略地将所有系数剔除。从计算的结果来看，对于Stein的无偏似然估计（SURE）和极大极小（minimaxi）原理的阈值选择规则，仅保存了约3%的系数；而其两种阈值选择规则，将所有的系数都变成了零。

同样地，对噪声进行小波分解时，也会产生高频系数，故一个信号的高频系数向量是有用信号和噪声信号的高频系数的叠加。由于SURE和minimaxi阈值选取规则较为保守(仅将部分系数置为零)，因此在信号的高频信息有很少一部分在噪声范围内时，这两种阈值非常有用，可以将弱小的信号提取出来。其他两种阈值选取规则， 在去除噪声时更为有效，但是也可能将有用信号的高频部分当作噪声信号去除掉。 

在实际的工程应用中，大多数信号可能包含着许多尖峰或突变，而且噪声信号也并不是平稳的白噪声。对这种信号进行消噪处理时，传统的傅里叶变换完全是在频域中对信号进行分析，它不能给出信号在某个时间点上的变化情况，因此分辨不出信号在时间轴上的任何一个突变。但是小波分析能同时在时频域内对信号进行分析，所以它能有效地区分信号中的突变部分和噪声，从而实现对非平稳信号的消噪。

下面通过一个实例，考察小波分析对非平稳信号的消噪。

【例4-7】利用小波分析对一个含噪的矩形波信号进行消噪处理。

%% 利用小波分析对一个含噪的矩形波信号进行消噪处理
% 设置信噪比和随机种子值
snr=4;	
init=2055615866;
% 产生原始信号sref和被高斯白噪卢污染的信号s
[sref,s]=wnoise(1,11,snr,init);

% 用sym8小波对信号s进行3层分解并对细节系数
% 选用sure阈值模式和尺度噪声
xd=wden(s,'heursure','s','one',3,'sym8');%对信号s进行启发式阈值，阈值类型为软阈值，输入参数scal为one，用sym8小波进行3级分解
% Wden函数：一维信号的小波消噪处理
% [xd，cxd，lxd]=wden（x，tptr，sorh，scal，n，‘wname’）；返回经过小波消噪处理后的信号xd及其小波分解结构。
% 1、x = ns为含噪信号
% 2、tptr = rigrsure：stein无偏估计； heursure：启发式阈值； sqtwolog：固定式阈值； minimaxi：极大极小值阈值
% 3、sorh = s：为软阈值； h： 为硬阈值；
% 4、输入参数scal规定了阈值处理随噪声水平的变化：
%   Scal=one：不随噪声水平变化。
%    Scal=sln：根据第一层小波分解的噪声水平估计进行调整。
%   Scal=mln：根据每一层小波分解的噪声水平估计进行调整。
% 5、N 是小波变换的级别。
% 6、WNAME 是指定为字符串的小波。

% 显示信号波形
subplot(3,1,1);
plot(sref);
xlabel('样本序号n');
ylabel('幅值A');
title('原始不含噪信号sref');
subplot(3,1,2);
plot(s);
title('含噪信号s');
subplot(3,1,3);
plot(xd);
title('小波消噪后的信号s');	
xlabel('样本序号n');
ylabel('幅值A');

工程中对时域信号进行处理
0 整个系统的流程框图

1 中值滤波
y=medfilt1(data,300);                 %中值滤波,滤除高频噪声
subplot(312);
plot(t,y);title('中值滤波后');

2 基线校准
%% --------------多项式拟合--基线校准------------------------------

opol = 6;
[p,s,mu] = polyfit(t,y,opol); 
f_y = polyval(p,t,[],mu); 
dt_y = y - f_y;
%将基线在0处的矫正后的数据迁移回原始数据附近(最终的目的相当于只在原始数据附近进行了基线矫正)
ZX = mean(f_y);                       %获得多项式值序列的均值
                                      %基线在0处的,基线矫正后的数据与原始数据对应时间点，在纵坐标上差异(迁移)的均值
dt_data = dt_y+ZX;                    %将基线在0处的矫正后的数据整体迁移回原始数据附近
%在同一个figure中画出原始数据以及矫正后的数据的图像，便于观察矫正效果
subplot(313)
% plot(t,y,'r*');                       %用红色※号线画出原始数据
% hold on;
plot(t,dt_data,'r-');                 %用蓝色线段画出基线在原始数据附近的，矫正后的数据
hold off;
title('基线校准后');

3.1 分区间峰值扫描
%% -------------------------查找波峰点-----------------------------
N=length(dt_data);
M=150;                         % 设置查找的分段数
dlength=floor(N/M);   % 将螺纹轮廓波形等分成 M 段, 每一段的长度
for k=1:M
    DATA{k}=data((k-1)*dlength+1:k*dlength);   % 保存每一段的 x 值
end

for k=1:M
    [bofeng(k),index1]=max(DATA{k});   % 寻找每一段里 y 的最大值
     bofeng(k)=DATA{k}(index1);         % 得到与波峰对应的横坐标 x
end
figure
plot(bofeng);
xlabel('划分区域的个数');
ylabel('波峰峰值');
title('在固定区间内查找峰值点');

4.1 峰值统计
%% ------门限判定，统计在所有划分区域内大于门限的所有波峰数目-----------

bofeng_door11=bofeng>2.5;
bofeng_door12=bofeng>4;
bofeng_door1=sum(sum(bofeng_door11));          %大于2.5的峰值,进行统计，认为是正常的脉搏波
bofeng_door2=sum(sum(bofeng_door12));          %大于4的峰值,进行统计，认为是异常的脉搏波
bofengnum=[bofeng_door11;bofeng_door12];
figure
bar(bofengnum,'stacked');

3.2&4.2 小波分解&小波重构
%% ---------------------小波分解与重构---------------------
[c,l]=wavedec(dt_data,3,'db2');                           %小波分解
approx=appcoef(c,l,'db2');                          
[cd1,cd2,cd3]=detcoef(c,l,[1,2,3]);                 %细节系数
figure
subplot(2,2,1)
plot(approx)
title('Approximation Coefficients-近似系数')
subplot(2,2,2)
plot(cd3)
title('Level 3 Detail Coefficients-3级细节系数')
subplot(2,2,3)
plot(cd2)
title('Level 2 Detail Coefficients-2级细节系数')
subplot(2,2,4)
plot(cd1)
title('Level 1 Detail Coefficients-1级细节系数')
x=waverec(c,l,'db2');                                %小波重建函数
figure
plot(x)
title('小波重建函数');

5 时频分析
本应是用小波进行时域的功率谱分析的，但是博主真的忘了怎么弄了，只能用STFT代替下，毕竟看的是功率。
%-------------STFT时频谱-------------
[S,F,T,P1]=spectrogram(x,window1,noverlap,nfft,Fs);
P1=10*log10(P1);
figure
% subplot(5,1,1)
surf(T,F,P1,'edgecolor','none'); axis tight;
view(0,90);
xlabel('Time (Seconds)'); ylabel('Hz');title('Time-frequency Analyse');

6 归一化功率谱密度
%------------------------绘制归一化功率谱---------------------
nMax = max(plot_Pxx1); nMin = min(plot_Pxx1);
ax_norm = (plot_Pxx1 - nMin)/(nMax - nMin);

figure
plot(ax_norm);
xlabel('频率(Hz)')
ylabel('归一化功率谱密度(dB)')
title('归一化小波功率谱')
grid on 

本次的程序介绍就到这里了，写的不好的地方，欢迎指正，欢迎私信或评论讨论。

	

前言：最近研究汽车碰撞的加速度信号，在信号的采集过程中难免遇到噪音，导致信号偏差，为了更好的反映系统情况，故常需要信号去噪，本文分享一些
常用信号平滑去噪的方法。
关键字：信号；去噪；Matlab
信号在实际测量中，难免会混入各种噪声。通常我们希望去除高频的随机噪声，或者是偏离正常测量太大的离群误差，以获得低频的测量数据。下面介绍几种常用的信号平滑去噪的方法。
1、移动平均法
滑动平均法(moving average)也叫做移动平均法、平均法、移动平均值滤波法等等，是一种时间域思想上的信号光滑方法。算法思路为，将该点附近的采样点做算数平均，作为这个点光滑后的值。
一般窗口为对称窗口，防止出现相位偏差。窗口一般为奇数。
以3点平均(窗口长度为3)公式为例，原数据为x，平滑后的数据为y：
y(n)=1/3∗(x(n−1)+x(n)+x(n+1))
对y(n)和y(n+1)相减，可以得到另一种计算形式：
2、Matlab内自带函数实现移动平均法
matlab有两个函数实现滑动平均法，一个是smoothdata()函数，一个是movmean()函数。
以窗口长度为5为例，smoothdata()函数调用方法为：
y = smoothdata( x , 'movmean' , 5 );
但是这个smoothdata函数实际上是调用了movmean()函数。所以如果直接使用的话，直接用movmean()会更快。
movmean()函数的调用方法为：
y = movmean( x , 5 );
下面以一个加噪声的正弦信号为例：
%移动平均滤波Nber_window = 3;%窗口长度(最好为奇数)t = 0:0.02:3;A = cos(2*pi*2*t)+0.5*rand(size(t));B1 = movmean(A,Nber_window);figure(1)plot(t,A,t,B1)
3、利用卷积函数conv()实现移动平均法
根据之前的公式，我们可以看到
y(n)=1/3∗(x(n−1)+x(n)+x(n+1))
就相当于一个对x和向量[1/3 1/3 1/3]做卷积。可以验证：
%移动平均滤波Nber_window = 3;%窗口长度(最好为奇数)t = 0:0.02:3;A = cos(2*pi*2*t)+0.5*rand(size(t));B1 = movmean(A,Nber_window);F_average = 1/N_window*ones(1,N_window);%构建卷积核B2 = conv(A,F_average,'same');%利用卷积的方法计算plot(t,B1,t,B2)
中间部分两者完全一致，但是两端有所差别。主要是因为，movmean()函数在处理边缘时，采用减小窗口的方式，而conv()相当于在两端补零。所以如何处理边缘也是值得注意的。
4、利用filter滤波函数实现移动平均法
首先介绍一下Z变换。以向前的滑动平均为例(这里中间值不是n而是n+1，所以相位会移动)。
y(n)=1/3∗(x(n)+x(n+1)+x(n+2))
它的Z变换可以简单的理解为，把x(n+k)替换为z^(-k)，即
因此对于filter滤波函数，输入的格式为：
y = filter(b,a,x)
其中b和a的定义为：
其中a1必须为1。所以对应的移动平均滤波可以表示为：
y = filter( [1/3,1/3,1/3] , [1] , x );
它和下面代码的是等价的(在边缘上的处理方式有所不同)
y = movmean( x , [2,0] );
代表有偏移的滑动平均滤波，2是向后2个点的意思，0是向前0个点的意思。
因为 filter滤波器使用有偏移的向后滤波。滤波后，相位会发生改变。所以通常采用零相位滤波器进行滤波，matlab内的函数为filtfilt()。原理从函数名字上就可以体现出来，就是先正常滤波，之后再将信号从后向前再次滤波。正滤一遍反滤一遍，使得相位偏移等于0。
5、移动平均的幅频响应
幅频响应可以通过之前4得到的H(z)函数来得到，在单位圆上采样，也就是把z替换为e^iw。
以中心窗口为例，
H(iw)的绝对值就是该滤波方法的幅频响应。以3点滤波为例，matlab代码为：
%计算频率响应N_window = 3;w = linspace(0,pi,400);H_iw = zeros(N_window,length(w));n=0;for k = -(N_window-1)/2:1:(N_window-1)/2    n = n+1;    H_iw(n,:) =1/3* exp(w.*1i).^(-k);%公式(e^iw)^(-k)endH_iw_Sum = abs(sum(H_iw,1));%相加figure()plot(w/2/pi,H_iw_Sum)
由于H变换在单位圆上的特性相当于傅里叶变换，所以上面代码等效于下面计算方法：
%计算频率响应N_window = 3;Y=zeros(400,1);Y(1:N_window) = 1/3;%设置滤波器Y_F = fft(Y);figure()plot(linspace(0,0.5,200),abs(Y_F(1:200)));
matlab也有自带的函数来看频率特性，freqz()，推荐使用这种。
其中，归一化频率等于信号频率除以采样频率f/Fs，采样频率等于时间采样间隔的倒数1/dt。对比不同窗口长度的幅频响应，可以看到：
1)平均所采用的点数越多，高频信号的滤波效果越好。
2)3点平均对于1/3频率的信号滤波效果最好，5点平均对1/5和2/5频率的信号滤波效果最好。所以根据这个特性，一方面我们要好好利用，一方面也要避免其影响。
举个应用的例子，比如你的采样频率为10Hz，采样3点滑动平均滤波，但是你的信号在3.3hz左右，你就会发现你的信号被过滤掉了，只留下了噪声。
反之，以美国近期的疫情为例，疫情的采样频率为1天一采样，而且显示出很强的7日一周期的特性(病毒也要过周末)。所以，为了消除这个归一化频率为1/7的噪声影响，采样7点的滑动平均滤波。可以看到所有以7天为一变化的信号分量全部被消除掉了。(下面这个图经常被引用，但是很少有人思考为什么用7天平均方法来平滑数据。)
回到原本的幅频特性问题上。当点数较少的时候，比如3个点，高频滤波效果并不是很好。所以当取的点数比较少的时候，需要平滑完一遍之后再平滑一遍，直到满意为止。多次平滑之后，高频的衰减非常明显。这也就是说，即使只有3个点平均，多次平滑之后也可以等效为一个较好的低通滤波器。
所以总结一下，移动平均滤波拥有保低频滤高频的特点，而且对于特点频率的滤波具有良好的效果。但是缺点是所有频率分量的信号都会有不同程度衰减。
6、时域和频域的转换关系
时域上的滤波和频域上的滤波是可以互相转换，且一一对应的。也就是时域上的卷积等于频域上的乘积。
下图为3点移动平均滤波法，时域和频域的转换关系：
虽然前面的 movmean()或者conv()等函数都是用时域实现的信号滤波，但是同样也可以完全在频域上实现。采用ifft(fft(x).*fft(F))实现的滤波效果，和完全时域上的滤波效果是等价的。
这也意味着你也可以在频域上操作，实现想要的滤波。比如想要低频通过高频衰减，就把fft后的信号，高频部分强行等于0即可。比如想要消除某个频率的信号(陷波)，就令fft后那个信号的频率等于0即可。同理，想要把振幅衰减1/2，就在对应频域上乘以0.5.