首页 > 其他分享 >5.信号处理(1) --常用信号平滑去噪的方法

5.信号处理(1) --常用信号平滑去噪的方法

时间:2022-11-12 10:02:41浏览次数:43  
标签:采样 函数 -- 信号处理 滤波 window 信号 movmean 平滑

前言:最近研究汽车碰撞的加速度信号,在信号的采集过程中难免遇到噪音,导致信号偏差,为了更好的反映系统情况,故常需要信号去噪,本文分享一些
常用信号平滑去噪的方法。

关键字:信号;去噪;Matlab


信号在实际测量中,难免会混入各种噪声。通常我们希望去除高频的随机噪声,或者是偏离正常测量太大的离群误差,以获得低频的测量数据。下面介绍几种常用的信号平滑去噪的方法。


1、移动平均法

滑动平均法(moving average)也叫做移动平均法、平均法、移动平均值滤波法等等,是一种时间域思想上的信号光滑方法。算法思路为,将该点附近的采样点做算数平均,作为这个点光滑后的值。

5.信号处理(1) --常用信号平滑去噪的方法_时域

一般窗口为对称窗口,防止出现相位偏差。窗口一般为奇数。

以3点平均(窗口长度为3)公式为例,原数据为x,平滑后的数据为y:

y(n)=1/3∗(x(n−1)+x(n)+x(n+1))

对y(n)和y(n+1)相减,可以得到另一种计算形式:

5.信号处理(1) --常用信号平滑去噪的方法_数据_02


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)

5.信号处理(1) --常用信号平滑去噪的方法_时域_03


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)

5.信号处理(1) --常用信号平滑去噪的方法_时域_04

中间部分两者完全一致,但是两端有所差别。主要是因为,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),即

5.信号处理(1) --常用信号平滑去噪的方法_数据_05


因此对于filter滤波函数,输入的格式为:

y = filter(b,a,x)

其中b和a的定义为:

5.信号处理(1) --常用信号平滑去噪的方法_频域_06

其中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.信号处理(1) --常用信号平滑去噪的方法_数据_07


5、移动平均的幅频响应

幅频响应可以通过之前4得到的H(z)函数来得到,在单位圆上采样,也就是把z替换为e^iw。

以中心窗口为例,

5.信号处理(1) --常用信号平滑去噪的方法_时域_08

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)
end

5.信号处理(1) --常用信号平滑去噪的方法_频域_09

由于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频率的信号滤波效果最好。所以根据这个特性,一方面我们要好好利用,一方面也要避免其影响。

5.信号处理(1) --常用信号平滑去噪的方法_数据_10

举个应用的例子,比如你的采样频率为10Hz,采样3点滑动平均滤波,但是你的信号在3.3hz左右,你就会发现你的信号被过滤掉了,只留下了噪声。

反之,以美国近期的疫情为例,疫情的采样频率为1天一采样,而且显示出很强的7日一周期的特性(也要过周末

5.信号处理(1) --常用信号平滑去噪的方法_数据_11

)。所以,为了消除这个归一化频率为1/7的噪声影响,采样7点的滑动平均滤波。可以看到所有以7天为一变化的信号分量全部被消除掉了。(下面这个图经常被引用,但是很少有人思考为什么用7天平均方法来平滑数据。)

5.信号处理(1) --常用信号平滑去噪的方法_数据_12

回到原本的幅频特性问题上。当点数较少的时候,比如3个点,高频滤波效果并不是很好。所以当取的点数比较少的时候,需要平滑完一遍之后再平滑一遍,直到满意为止。多次平滑之后,高频的衰减非常明显。这也就是说,即使只有3个点平均,多次平滑之后也可以等效为一个较好的低通滤波器。

所以总结一下,移动平均滤波拥有保低频滤高频的特点,而且对于特点频率的滤波具有良好的效果。但是缺点是所有频率分量的信号都会有不同程度衰减。


6、时域和频域的转换关系

时域上的滤波和频域上的滤波是可以互相转换,且一一对应的。也就是时域上的卷积等于频域上的乘积。

下图为3点移动平均滤波法,时域和频域的转换关系:

5.信号处理(1) --常用信号平滑去噪的方法_数据_13

虽然前面的 movmean()或者conv()等函数都是用时域实现的信号滤波,但是同样也可以完全在频域上实现。采用ifft(fft(x).*fft(F))实现的滤波效果,和完全时域上的滤波效果是等价的。

5.信号处理(1) --常用信号平滑去噪的方法_时域_14

5.信号处理(1) --常用信号平滑去噪的方法_时域_15

5.信号处理(1) --常用信号平滑去噪的方法_数据_16

这也意味着你也可以在频域上操作,实现想要的滤波。比如想要低频通过高频衰减,就把fft后的信号,高频部分强行等于0即可。比如想要消除某个频率的信号(陷波),就令fft后那个信号的频率等于0即可。同理,想要把振幅衰减1/2,就在对应频域上乘以0.5.

标签:采样,函数,--,信号处理,滤波,window,信号,movmean,平滑
From: https://blog.51cto.com/domi/5846529

相关文章

  • 8.图像视觉(1) --鱼眼摄像头标定和图像畸变矫正
    前言:最近项目上研究鱼眼摄像头的画面畸变问题,对比了基于Matlab和PythonOpencv的方法,分别进行了摄像头的标定和图像矫正,实际结果个人认为Opencv的效果为佳,本文分享一下基于M......
  • 元进Python(4)---networkx制作我最喜爱55部美剧天体排名
    利用networkx和matplotlib给我最喜爱的55部美剧来个天体排名【最终呈现效果】【代码实现及注释】importnetworkxasnximportmatplotlib.pyplotaspltUSTVNode=['Everybo......
  • 学Java· 从new说对象实例化
    概念了解想要知道内存中发生了什么,我们先来了解两个内存概念......
  • 计算机网络:随机访问介质访问控制之令牌传递协议
    在轮询访问中,用户不能随机地发送信息,而要通过一个集中控制的监控站,以循环方式轮询每个结点,再决定信道的分配。当某结点使用信道时,其他结点都不能使用信道。典型的轮询访问介......
  • Codeforces 1738 / Codeforces Global Round 22
    目录ContestLinkProblemAGloryAddictsProblemBPrefixSumAddictsProblemCEvenNumberAddictsProblemDPermutationAddictsProblemEBalanceAddictsProblemF......
  • 冒泡排序(数组中的问题)
    问题:使用冒泡排序的方法,将数组中的元素按照升序的方式将其排列。冒泡排序核心思想:两两相邻元素进行比较,满足条件则交换;     ①先确认趟数;     ②写下一趟冒泡......
  • 那些年被误解的指针和数组
    误解1:&运算符返回一个地址解释:  &叫做取址运算符,运算的结果是返回一个指向某个数据类型对象的指针。    inta=1; int*p=&a;       &a不是地址,&a是......
  • 「题解」Codeforces 1342F Make It Ascending
    只会\(\mathcal{O}(3^nn^2)\),打开题解一看怎么还真是这个玩意/jy首先集合之间形成一个sum和pos的二维偏序,那么思路就是对一维扫描线,然后另一维搞个什么东西。具体到......
  • Linux下的网络管理工具—OpenNMS
    OpenNMS的是一个运营商级别的,高度集成的,开放源码的平台,用于构建网络监控解决方案。OpenNMS有两个发行版:MeridianandHorizon。使用Meridian是可取的,对企业提供稳定和长期的......
  • XAMPP phpmyadmin 设置局域网访问
    1.问题描述本地local和127.0.0.1访问都正常,局域网访问失败,显示以下内容:2.检查用户的主机名是否是任意主机(%)3.修改httpd-xampp.conf配置文件找到phpmyadmin这一行,......