精华内容
下载资源
问答
  • 几种常用信号平滑去噪的方法(附Matlab代码)

    万次阅读 多人点赞 2020-07-31 20:44:36
    几种常用信号平滑去噪的方法(附Matlab代码)1 滑动平均法1.0 移动平均法的方法原理1.1 matlab内自带函数实现移动平均法1.2 利用卷积函数conv()实现移动平均法1.3 利用filter滤波函数实现移动平均法1.4 移动平均的...

    2020年8月更新:增加一个时域和频域的转换关系图,增加相应小节1.5

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

    还是惯例声明:本人没有相关的工程应用经验,只是纯粹对相关算法感兴趣才写此博客。所以如果有错误,欢迎在评论区指正,不胜感激。本文主要关注于算法的实现,对于实际应用等问题本人没有任何经验,所以也不再涉及。

    1 滑动平均法

    本章参考目录
    1《数字信号处理》-胡广书
    2《Digital Signal Processing - A Practical Guide For Engineers and Scientists》 - Steven W.Smith

    1.0 移动平均法的方法原理

    作为开篇第一个方法,会夹带一些数字信号处理的基本方法,可能会导致篇幅比较啰嗦,之后几章我会尽量挑重点的讲解。

    滑动平均法(moving average)也叫做移动平均法、平均法、移动平均值滤波法等等,是一种时间域思想上的信号光滑方法。算法思路为,将该点附近的采样点做算数平均,作为这个点光滑后的值。
    在这里插入图片描述
    一般窗口为对称窗口,防止出现相位偏差。窗口一般为奇数。
    以3点平均(窗口长度为3)公式为例,原数据为x,平滑后的数据为y:
    y(n)=1/3(x(n1)+x(n)+x(n+1)) y(n)=1/3*( x(n-1)+x(n)+x(n+1) )
    对y(n)和y(n+1)相减,可以得到另一种计算形式:
    y(n+1)=y(n)13x(n1)+13x(n+2) y(n+1)=y(n)-\frac{1}{3}x(n-1)+\frac{1}{3}x(n+2)
    当然这两者都是等价的。

    1.1 matlab内自带函数实现移动平均法

    matlab有两个函数实现滑动平均法,一个是smoothdata()函数,一个是movmean()函数。
    以窗口长度为5为例,smoothdata()函数调用方法为:

    y = smoothdata( x , 'movmean' , 5 );
    

    但是这个smoothdata函数实际上是调用了movmean()函数。所以如果直接使用的话,直接用movmean()会更快。
    movmean()函数的调用方法为:

    y = movmean( x , 5 );
    

    下面以一个加噪声的正弦信号为例:

    %移动平均滤波
    clear
    clc
    close all
    
    N_window = 5;%窗口长度(最好为奇数)
    t = 0:0.1:10;
    A = cos(2*pi*0.5*t)+0.3*rand(size(t));
    B1 = movmean(A,N_window);
    figure(1)
    plot(t,A,t,B1)
    

    在这里插入图片描述

    1.2 利用卷积函数conv()实现移动平均法

    根据之前的公式,我们可以看到
    y(n)=1/3(x(n1)+x(n)+x(n+1)) y(n)=1/3*( x(n-1)+x(n)+x(n+1) )
    就相当于一个对x和向量[1/3 1/3 1/3]做卷积。
    可以验证:

    N_window = 5;%窗口长度(最好为奇数)
    t = 0:0.25:10;%时间
    A = cos(2*pi*0.5*t)+0.3*rand(size(t));
    B1 = movmean(A,N_window);
    
    F_average = 1/N_window*ones(1,N_window);%构建卷积核
    B2 = conv(A,F_average,'same');%利用卷积的方法计算
    figure(2)
    plot(t,B1,t,B2)
    

    在这里插入图片描述
    中间部分两者完全一致,但是两端有所差别。主要是因为,movmean()函数在处理边缘时,采用减小窗口的方式,而conv()相当于在两端补零。所以如何处理边缘也是值得注意的。

    1.3 利用filter滤波函数实现移动平均法

    首先介绍一下Z变换。以向前的滑动平均为例(这里中间值不是n而是n+1,所以相位会移动)。
    y(n)=1/3(x(n)+x(n+1)+x(n+2)) y(n)=1/3*( x(n)+x(n+1)+x(n+2) )
    它的Z变换可以简单的理解为,把x(n+k)替换为z^(-k),即
    H(z)=13(z0+z1+z2) H(z)=\frac{1}{3}(z^{0}+z^{-1}+z^{-2})

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

    y = filter(b,a,x)
    

    其中b和a的定义为:
    H(z)=b1+b2z1+b3z2++bnzn+11+a2z1+a3z2++amzm+1 H(z)=\frac{b_1+b_2\cdot z^{-1}+b_3\cdot z^{-2}+\cdots +b_{n}\cdot z^{-n+1}}{1+a_2\cdot z^{-1}+a_3\cdot z^{-2}+\cdots +a_{m}\cdot z^{-m+1}}
    其中a1必须为1。所以对应的移动平均滤波可以表示为:

    y = filter( [1/5,1/5,1/5,1/5,1/5] , [1] , x );
    

    它和下面代码的是等价的(在边缘上的处理方式有所不同)

    y = movmean( x , [4,0] );
    

    代表有偏移的滑动平均滤波,4是向后4个点的意思,0是向前0个点的意思。

    因为 filter滤波器使用有偏移的向后滤波。滤波后,相位会发生改变。所以通常采用零相位滤波器进行滤波,matlab内的函数为filtfilt()。原理从函数名字上就可以体现出来,就是先正常滤波,之后再将信号从后向前再次滤波。正滤一遍反滤一遍,使得相位偏移等于0。这个之后再IIR滤波器会讲一下。

    1.4 移动平均的幅频响应

    幅频响应可以通过之前1.3得到的H(z)函数来得到,在单位圆上采样,也就是把z替换为e^iw。
    以中心窗口为例,
    y(n)=13(x(n1)+x(n)+x(n+1)) y(n)=\frac{1}{3} ( x(n-1)+x(n)+x(n+1) )
    H(iw)=13(exp(iw)1+exp(iw)0+exp(iw)1) H(iw)=\frac{1}{3}(\exp(iw)^{1}+\exp(iw)^{0}+\exp(iw)^{-1})
    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
    H_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个点,高频滤波效果并不是很好。所以当取的点数比较少的时候,需要平滑完一遍之后再平滑一遍,直到满意为止。这个原理也可以通过幅频特性来解释,多次平滑相当于乘了多个H(z)。对于平滑2次,它的图像也就是abs(sum(H_iw,1).*sum(H_iw,1));对于平滑4次,它的图像相当于乘以四个sum(H_iw,1)。(注:因为时域上的卷积等于频域上的乘积,多次卷积对应着多次乘积。)
    在这里插入图片描述
    可以看到,多次平滑之后,高频的衰减非常明显。这也就是说,即使只有3个点平均,多次平滑之后也可以等效为一个较好的低通滤波器。

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

    1.5 时域和频域的转换关系

    额外补充一部分小内容,可能前面有些概念加入的太突然。很多人可能觉得之前时域上的平均法非常好理解,为什么突然加入幅频特性图,又是Z变换又是fft的。

    其实时域上的滤波和频域上的滤波是可以互相转换,且一一对应的。也就是时域上的卷积等于频域上的乘积。
    下图为3点移动平均滤波法,时域和频域的转换关系:
    在这里插入图片描述
    同样的滤波操作,可以用时域公式: y(n)=13(x(n1)+x(n)+x(n+1))y(n)=\frac{1}{3} ( x(n-1)+x(n)+x(n+1) ),进行描述。也可以用频域上,滤波器的幅频特性进行描述。

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

    下面是展示了窗口长度为3的平滑滤波,从时域上和频域上对信号进行滤波的对比:

    %实验,检验频域和时域的一致性
    %以3点滤波为例
    clear
    clc
    close all
    
    N_window = 3;%窗口长度(最好为奇数)
    t = 0:0.1:10;
    A = cos(2*pi*0.3*t)+0.1*cos(2*pi*5*t)+0.2*randn(size(t));
    
    F_average = 1/N_window*ones(1,N_window);%创建滤波器
    B2 = conv(A,F_average,'same');%利用卷积的方法计算
    
    figure(1)
    plot(t,A,'k','linewidth',0.8)
    
    %计算原信号的fft
    A_fft=fft(A);
    %构建频域上的滤波器
    F_average2=zeros(size(t));%长度与x相同,为了后面.*运算
    F_average2(1:(N_window-1)/2+1) = 1/N_window;
    F_average2(end-(N_window-1)/2+1:end) = 1/N_window;%前后设置对称,使得相位不变
    F_Fft = fft(F_average2);
    
    figure(2)
    subplot(1,2,1)
    plot(linspace(0,1,length(t)),abs(A_fft),'linewidth',1);
    subplot(1,2,2)
    plot(linspace(0,1,length(t)),abs(F_Fft),'linewidth',1);
    
    %进行反逆变换
    B3=ifft(A_fft.*F_Fft);
    
    figure()
    plot( t,B2,t,B3 )%对比时域操作和频域操作的差异
    

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

    2 Savitzky-Golay法

    本章参考目录:
    1《数字信号处理》-胡广书
    2 Savitzky-Golay平滑滤波器的最小二乘拟合原理综述 https://wenku.baidu.com/view/b63017eed0d233d4b04e690e.html?fr=search

    2.1 Savitzky-Golay法的方法原理

    Savitzky-Golay法,又叫做平滑滤波器,最著名的就是5点3次滤波器。这是一种基于时间域上的多项式拟合,来消除噪声的方法。

    以简单的5点2次构造方法为例,介绍一下基本的求解系数的方法:
    首先选取5个点x[-2],x[-1],x[0],x[1],x[2],根据这5个点,构造一条2次抛物线:
    f(i)=a0+a1i+a2i2 f(i)=a_{0}+a_{1} \cdot i+a_{2} \cdot i^2
    这里i=-2,-1,0,1,2。因此,我们需要寻找最优的a0,a1,a2,使得最小二乘拟合最小。最小二乘拟合的函数为:
    E=(f(i)x(2))2=(f(2)x(2))2+(f(1)x(1))2++(f(0)x(0))2+(f(1)x(1))2+(f(2)x(2))2 E=\sum{(f(i)-x(-2))^2}\\ =(f(-2)-x(-2))^2+(f(-1)-x(-1))^2+\\ +(f(0)-x(0))^2+(f(1)-x(1))^2+(f(2)-x(2))^2
    之后希望最小二乘E最小,我们使得其导数等于0,也就是
    Eap=0 \frac{\partial E}{\partial a_p}=0
    上面式子中,a0、a1和a2三个偏导数,可以得到3个方程。联立,即可求得a0、a1和a2。
    对于无相位差的滤波,我们希望窗口是对称的。也就是说,我们希望用5个点,去估计f[0]的值。因此我们需要的只有a0,因为:
    f(0)=a0+a10+a20=a0 f(0)=a_{0}+a_{1} \cdot 0+a_{2} \cdot 0=a_{0}
    这里借助Mathematica的符号运算进行编程求解,程序如下:

    Clear["Global`*"]
    M = 5;(*M个点*)
    P = 2;(*P次*)
    
    Mk = (M - 1)/2;
    f[i_] = Sum[a[k]*i^k, {k, 0, P}]
    En = Sum[(f[k] - x[k])^2, {k, -Mk, Mk}];
    DEn = D[En, {Array[a, P + 1, {0, P}]}];
    s = Solve[DEn == 0, Array[a, P + 1, {0, P}]];
    s[[1]][[1]]
    

    可以得到结果:

    a[0] -> 1/35 (-3 x[-2] + 12 x[-1] + 17 x[0] + 12 x[1] - 3 x[2])
    

    对于经典的5点3次平滑结果,只需要把上面代码P改为3即可,恰好结果相同:

    a[0] -> 1/35 (-3 x[-2] + 12 x[-1] + 17 x[0] + 12 x[1] - 3 x[2])
    

    Savitzky-Golay法可以看做是移动平均法的推广。因为当次数P等于1时,Savitzky-Golay法退化为求解平均值的形态。当然,它也可以视为一种加权平均。其它加权平均比如高斯加权等等,可以在smoothdata函数里的帮助文档看到,就不细说了。

    除此之外,也可以直接用matlab中的sgolay()函数,直接得到系数,还是以5点3次为例:

    M=5;%5点
    P=3;%3次
    b = sgolay( P , M );
    a0 = b((M+1)/2,:)
    

    2.2 Savitzky-Golay法的matlab实现

    matlab中实现方法和第一章中滑动滤波的方法相似。
    利用matlab自带的smoothdata(A,‘sgolay’)函数就可以实现Savitzky-Golay法滤波。但是,该函数只支持N点2次的滤波(2.1阶证明了5点2次和5点3次的a0滤波系数相同,但是其他点数则不相同,这一点一定注意)。

    如果想得到不同点数不同次数的公式,参考2.1中的sgolay()函数,可以得到不同的系数。下面代码演示了一个利用matlab自带的5点2阶插值和自己编程实现的5点3次插值,其中自己编程的部分采用了1.2节部分的卷积运算,和2.1节介绍的matlab自带函数sgolay()。

    %sgolay滤波
    clear
    clc
    close all
    
    N_window = 5;%窗口长度(最好为奇数)
    
    t = 0:0.1:10;
    A = cos(2*pi*0.5*t)+0.4*rand(size(t));
    
    %matlab自带的5点2次插值
    B1 = smoothdata(A,'sgolay' ,N_window);
    figure(1)
    plot(t,A,t,B1)
    
    %5点3次插值
    B2 = smooth_SG_hyh(A,5,3,1);
    figure(2)
    plot(t,B1,t,B2)
    
    
    function y_new = smooth_SG_hyh(y,M,P,n)
    %M点P次插值
    %M为窗口长度
    %P为拟合阶数
    %n为光滑n次
    
    m = length(y);
    N_window = M;%窗口长度
    b = sgolay( P , N_window );
    F_SG = b((N_window+1)/2,:);%5点3次插值
    y_new = y;
    for k=1:n
        y_new2 = wextend(1,'sym',y_new,(N_window-1)/2);%镜像延拓
        y_new2 = conv(y_new2,F_SG,'same');%利用卷积的方法计算
        y_new = wkeep1(y_new2,m);%中间截断
    end
    
    end
    

    在这里插入图片描述

    2.3 Savitzky-Golay法的幅频响应

    利用之前1.4的内容,可以求出Savitzky-Golay法的幅频响应。

    %频率响应对比
    [w,A_W_5_3] = A_W_SG_hyh(5,3);
    [w,A_W_5_1] = A_W_SG_hyh(5,1);
    figure()
    plot(w,A_W_5_1,w,A_W_5_3)
    
    %频率响应对比
    [w,A_W_7_3] = A_W_SG_hyh(7,3);
    [w,A_W_7_1] = A_W_SG_hyh(7,1);
    figure()
    plot(w,A_W_7_1,w,A_W_7_3)
    ylim([0,1])
    
    function [w_2p,H_iw_Sum] = A_W_SG_hyh(M,P)
    N_window=M;
    b = sgolay( P , N_window );
    F_SG = b((N_window+1)/2,:);%5点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,:) =F_SG(n)* exp(w.*1i).^(-k);%公式(e^iw)^(-k)
    end
    H_iw_Sum = abs(sum(H_iw,1));%相加
    w_2p = w/2/pi;
    end
    

    在这里插入图片描述
    可以看到,当拟合的次数增加,该方法对低频的滤波效果变差。所以和平均法相比,5点3次方法的频率特性与3点移动平均比较相似。因此,5点3次方法适用于噪声频率比较高,信号频率比较低(说的都是归一化频率,也就是和采样频率之比)的情况。相比较相同频率效果的平均法,Savitzky-Golay法用了更多的点,在时域上保留了更好信息。

    3 处理离群值(粗大误差)的方法

    离群值是指在测量值中,出现了一些反常的值,这个反常值与附近的其它正常值差别非常大。

    常用的方法有利用3σ3 \sigma、中位值法、以及基于中位值的MAD方法等等。

    本章参考目录:
    1 离群值处理方法 https://blog.csdn.net/zjz199303/article/details/79135530
    2 常用去除离群值的算法 https://blog.csdn.net/dulingwen/article/details/97006884

    3.1 中位值法

    中位值法,也叫移动中位数法、中值滤波法等。其思想是将窗口内的中位数作为输出结果,如下图所示:
    在这里插入图片描述
    优点是,在数据采样点密集,且比较平滑的情况下,中位数法可以很好地剔除离群值。缺点是不适用于噪声较大的情况。而且平滑之后,数据光滑度不足。经过中位值法处理之后,极值点会丢失。

    matlab中自带的movmedian()函数可以实现中位值法滤波。

    clear
    clc
    close all
    
    %离群值的删除
    t = 0:0.05:10;
    A = sin(t);
    Ri = randi(length(t),10,1);
    A(Ri) = A(Ri)*2;
    
    %移动中位数,窗口设置为7
    B1 = movmedian(A,7);
    figure(1)
    plot(t,A,t,B1)
    

    下图可以看到中位值法对于离群值处理的结果。
    在这里插入图片描述
    虽然中值方法能够将所有离群值筛除,但是在8s处的波峰位置,极值点被不属于中位数,所以被消除了。在2s处,由于噪声比较复杂,所以出现了信号不光滑的现象(这里如果加上高频噪声,不光滑现象会更明显)。

    3.2 标准差法和MAD法

    标准差法又叫做3σ3 \sigma法。目的是规定一个数据波动阈值,当数据超过这个阈值的时候,便认为该数据离群。这个方法阈值的选取方法,采用窗口数据的3倍标准差。

    MAD法也是定义了一个阈值,这个阈值叫做中位数绝对偏差MAD。如果超过了3倍的MAD,则认为该数据离群。

    两者在Matlab里,可以用filloutliers()函数进行实现。下面代码对比了标准差法和MAD法在消除离群值的效果。

    clear
    clc
    close all
    
    %离群值的删除
    t = 0:0.06:10;
    A = sin(t)+0.1*rand(size(t));
    
    Ri = randi(length(t),4,1);
    A2 = A;
    A2(Ri) = A(Ri)*3;
    
    figure(3)
    B2 = filloutliers(A2,'linear','movmean',11);%标准差法
    B3 = filloutliers(A2,'linear','movmedian',11);%MAD法
    
    subplot(3,1,1)
    plot(t,A2)
    YL = ylim;
    subplot(3,1,2)
    plot(t,B2)
    ylim(YL)
    subplot(3,1,3)
    plot(t,B3)
    ylim(YL)
    

    下图可以看到,原本4个离群值,标准差法只找到了1个,但是MAD方法能够做到全部找到,说明MAD方法比标准差法的效果更好,一般更推荐用MAD方法。
    在这里插入图片描述

    3.3 Matlab中其它离群值消除方法

    matlab自带的函数smoothdata还有两种方法,可以兼顾去噪和去除离群噪声,一个是’rlowess’ 方法,一个是’rloess’ 方法。下面以’rloess’ 方法为例:

    t = 0:0.06:10;
    A = sin(t)+0.2*rand(size(t));
    
    Ri = randi(length(t),4,1);
    A2=A;
    A2(Ri) = A(Ri)*3;
    figure(4)
    B4 = smoothdata(A2,'rloess',8) ;
    plot(t,A2,t,B4)
    legend('原函数','光滑')
    

    下图可以看到,数据在光滑的同时,顺便也把离群值去掉了。
    在这里插入图片描述

    4 其它一些FIR滤波器实现光滑去噪

    本章参考:
    Matlab数字滤波入门 https://zhuanlan.zhihu.com/p/65483011

    4.1 FIR和IIR的区别

    FIR滤波器也叫做 有限冲击响应滤波器。和它相对的是IIR 无限冲击响应滤波器。在matlab里,有一张图可以比较直观。
    在这里插入图片描述
    对于FIR滤波器,用到的数据只有输入项b,输出项为1。但是对于IIR滤波器,不仅有输入项b(分子),还有输出项a(分母)。

    对于FIR滤波器,由于只有分子b,所以可以用卷积conv()函数去进行滤波。但是对于IIR,就不能用卷积函数去计算。一般认为,IIR函数在相同阶数下,滤波效果要比FIR函数要好,但是有可能出现发散问题。

    比如前面的移动平均滤波中的1.3章所介绍的,b=[1/3,1/3,1/3],a=1。就属于一个典型的FIR滤波器。其中,conv(x,b)卷积方法相当于无相位偏移的中心滤波,
    y(n)=13(x(n1)+x(n)+x(n+1)) y(n)=\frac{1}{3} ( x(n-1)+x(n)+x(n+1) )
    而filter(b,a,x)函数相当于有相位偏移的向后滤波
    y(n)=13(x(n)+x(n+1)+x(n+2)) y(n)=\frac{1}{3} ( x(n)+x(n+1)+x(n+2) )

    事实上,通过下面的等式,我们还可以构建一个等价的IIR滤波器:
    H(z)=13(z0+z1+z2)=131z31z1 H(z)=\frac{1}{3}(z^{0}+z^{-1}+z^{-2})=\frac{1}{3}\frac{1-z^{-3}}{1-z^{-1}}
    此时b=[1,0,0,0,-1],a=[1,-1]。

    Fs=20;
    t = 0:1/Fs:10;
    A = cos(2*pi*0.3*t)+1*sin(2*pi*1.2*t)+0.8*rand(size(t));
    
    b=[1/3,1/3,1/3];a=1;
    B1 = filter(b,a,A);
    
    b=1/3*[1,0,0,0,-1];a=[1,-1];
    B2 = filter(b,a,A);
    
    plot(t,B1,t,B2)
    

    对比这两个滤波器,可以看到滤波效果是基本等价的。
    在这里插入图片描述

    4.2 利用Matlab构建FIR滤波器

    回到FIR滤波器上,之前提到的时域角度提出的均值滤波器和Savitzky-Golay滤波器,都属于FIR滤波器。他们在频域上的表现并不能随意的控制。

    如果我们希望构建一个滤波器,让它在频域上满足自己的要求,则可以用fir1()函数去构建。我这里只演示一个简单地,更专业的就不在这里献丑了。

    比如,我的采样频率为10hz,我想提取的信号在1hz附近,而噪声则大于等于2hz。因此,需要构建一个低频通过,高频减弱的滤波器,使得在1.5hz以下的信号能够保留,1.5hz以上的信号删除。

    于是用下面代码去构建,滤波器长度为31,保留0-0.15之内的信号。(注1:matlab的归一化信号不是0-0.5,而是0-1,所以对应的0.15需要乘以2倍,变成0.3)(注2:最小频率不能是0,必须大于0,所以用0.001代替)(注3:滤波器长度为31是为了演示,实际上不需要这么大的长度,滤波器长度大了的话计算量也会增大)

    b = fir1(31,[0.001,0.3]);
    

    构建的滤波器形状和频率特性为:
    在这里插入图片描述
    滤波器形状很像sin(x)/x函数曲线。因为如果频率特性为理想的阶梯形状,则滤波器时域形态从数学上求解就是sin(x)/x函数,也被简写为sinc(x)函数。这种滤波器叫做sinc滤波器,由于其频域像一个矩形,也叫作矩形滤波器。特点就是低频完全通过,高频全部衰减为0,是一种理想的滤波器。

    频率特性满足预期的设计。滤波效果如下:
    在这里插入图片描述
    附上上面的matlab代码

    clear
    clc
    close all
    
    Fs=10;
    t = 0:1/Fs:10;
    A = cos(2*pi*1*t)+0.8*sin(2*pi*2*t+0.5*randn(size(t)))+0.5*randn(size(t));
    
    b = fir1(31,[0.001,0.3]);
    
    figure()
    subplot(1,2,1)
    N=length(b);
    stem(-(N-1)/2:1:(N-1)/2,b)%绘制滤波器形状
    subplot(1,2,2)
    [h,w]=freqz(b,1,512);%输出频率特性
    plot(w/2/pi,abs(h))
    xlim([0,0.5]);ylim([0,1])
    
    
    B2 = conv(A,b,'same');%利用卷积的方法计算
    %因为FIR滤波可以利用卷积方法代替,这两个等效的
    
    figure()
    plot(t,A,t,B2)
    legend('原函数','滤波')
    

    5 IIR滤波器实现光滑去噪

    比较著名有巴特沃斯滤波器butter()和切比雪夫滤波器cheby1()。这里感觉还有地方没看懂,就不详细写了。matlab里可以通过designfilt()来自定义的构建滤波器。

    以butter滤波器为例,采用零相位滤波,和第4.2节一样,希望保留0-0.15之内的信号。那么构建滤波器格式为:

    [b,a] = butter(6,0.15*2,'low' );
    

    b为分子,a为分母。6为6阶,代表滤波器的长度为6+1个点。'low’是低通滤波器的意思,0.15*2代表保留0-0.15之内的信号。得到的滤波器幅频响应曲线如下,可以看到能够满足要求。
    在这里插入图片描述
    利用butter滤波器进行光滑去噪的matlab代码如下:

    %butter滤波器
    clear
    clc
    close all
    
    
    Fs = 10;
    t = 0:1/Fs:10;
    A = cos(2*pi*1*t)+0.8*sin(2*pi*2*t+0.5*randn(size(t)))+0.5*randn(size(t));
    [b,a] = butter(6,0.15*2,'low' );
    B = filtfilt(b,a,A);
    
    figure()
    [h,w]=freqz(b,a,512);
    plot(w/2/pi,abs(h))
    xlim([0,0.5]);ylim([0,1])
    
    
    figure()
    plot(t,A,t,B)
    legend('原函数','滤波')
    

    平滑后的曲线效果如下图。
    在这里插入图片描述

    展开全文
  • 数字图像处理的课程设计,是关于图像的平滑去噪,包括高斯滤波,中值滤波以及箱式滤波,使用的是MFC界面
  • PCNN平滑去噪滤波程序

    2013-06-11 21:07:32
    这是一个基于PCNN的平滑滤波程序,可以有效对椒盐噪声造成干扰的图像平滑去噪
  • 图像平滑去噪

    2016-08-06 16:02:59
    图像平滑去噪,形态学处理
  • bregman matlab兼顾图像平滑去噪与边缘保留的自适应全变分模型
  • bregman matlab兼顾图像平滑去噪与边缘保留的自适应全变分模型matlab
  • 记录自己用python加opencv实现的图像处理的入门操作,各种平滑去噪滤波器的实现。 包括有:产生的椒盐噪声、高斯噪声等等,以及使用的中值滤波、平均滤波、高斯滤波等等。 分成了两部分来实现:一是自编写函数来实现...

    记录自己用python加opencv实现的图像处理的入门操作,各种平滑去噪滤波器的实现。
    包括有:产生的椒盐噪声、高斯噪声等等,以及使用的中值滤波、平均滤波、高斯滤波等等。
    分成了两部分来实现:一是自编写函数来实现,二是调用opencv中的相应函数,对比效果。

    噪声的产生:分别是椒盐噪声和高斯噪声,原理的话可以参考别人的博客或我之后再补充,噪声就是在原来的图像上以一定的特殊规律给图像增添一些像素,使图像变得模糊等等。

    1.噪声的自编写

    椒盐噪声:
    输入图像和自定义的噪声阈值,输出处理后的图像

    # 向图片中添加椒盐噪声
    def salt_pepper_noise(image, prob):  # prob:盐噪声阈值,由用户自己决定
        output = np.zeros(image.shape, np.uint8)
        thres = 1 - prob  # 胡椒噪声阈值
        for i in range(image.shape[0]):  # 遍历整个图片的灰度级
            for j in range(image.shape[1]):
                randomnum = random.random()  # 生成一个随机0-1之间的随机数
                if randomnum < prob:  # 如果随机数大于盐噪声阈值0.1,则将此位置灰度级的值设为0,即添加盐噪声
                    output[i][j] = 0
                elif randomnum > thres:  # 如果随机数大于胡椒噪声阈值1-0.1,则将此位置灰度级的输出设为255,即添加胡椒噪声
                    output[i][j] = 255
                else:  # 如果随机数处于两者之间,则此位置的灰度级的值等于原图的灰度级值
                    output[i][j] = image[i][j]
        return output
    

    高斯噪声:
    输入图像,自定义的参数:均值、方差,输出的是处理后的图像

    # 向图片中添加高斯噪声
    def gasuss_noise(image, mean=0, var=0.001):         # mean : 均值,var : 方差
        image = np.array(image/255, dtype=float)
        noise = np.random.normal(mean, var ** 0.5, image.shape)  # 使用numpy库中的函数生成正态分布矩阵,对应数据分别为概率均值,概率标准差,图像的大小
        output = image + noise  # 输出结果为原图灰度级概率与噪声概率相加
        output_handle = np.array([[[0]*3 for i in range(output.shape[1])] for i in range(output.shape[0])], dtype=float)
        # 处理后最终输出矩阵将齐大小设置为与原图一样
        if output.min() < 0:  # 确定一个比较中间值
            low_clip = -1.
        else:
            low_clip = 0.
        for i in range (output.shape[0]):  # 遍历整个三位矩阵
            for j in range (output.shape[1]):
                for k in range (output.shape[2]):
                    if output[i][j][k] < low_clip:  # 将输出的概率矩阵内的值限定在(-1,1)范围内
                        output_handle[i][j][k] = low_clip   # 使其之后*255变为灰度级时不会超出[0-255]的范围
                    elif output[i][j][k] > 1.0:
                        output_handle[i][j][k] = 1.0
                    else:
                        output_handle[i][j][k] = output[i][j][k]    # 在最大值和最小值之间的不变
        output = np.uint8(output_handle*255)   # 将处理后的灰度级转化为[0-255]的整数级
        return output
    

    加性噪声:

    # 向图片中添加加性噪声
    def addrandom_noise(image,prob=0.1):
        output = image   # 将原始图像数据拷贝至输出矩阵
        n = random.randint(1, 1000) + int(prob*20000)
        for k in range(n-500):
            a = random.randint(0, 50)
            b = random.randint(0, 50)
            c = random.randint(0, 50)
            i = random.randint(0, image.shape[0]-1)
            j = random.randint(0, image.shape[1]-1)
            output[i][j][0] = 255-a
            output[i][j][1] = 255-b
            output[i][j][2] = 255-c
        for k in range(n):
            a = random.randint(0, 50)
            b = random.randint(0, 50)
            c = random.randint(0, 50)
            i = random.randint(0, image.shape[0]-1)
            j = random.randint(0, image.shape[1]-1)
            output[i][j][0] = a
            output[i][j][1] = b
            output[i][j][2] = c
        return output
    

    2.平滑去噪滤波器的自编写

    设计了均值滤波、中值滤波、高斯滤波,可分别对应处理以上的三种噪声,效果较好。
    具体的对应情况可自己比较一下哦。

    中值滤波:

    # 中值滤波 a为要处理的图像  windowsize为采用的模版大小
    def medianfliter(a, windowsize):
        output = a
        if windowsize == 3 :
            output1 = np.zeros(a.shape, np.uint8)
            for i in range(1, output.shape[0]-1):  # 求齐周围9个方格与模版进行冒泡排序
                for j in range(1, output.shape[1]-1):
                    value1 = [output[i-1][j-1], output[i-1][j], output[i-1][j+1], output[i][j-1], output[i][j], output[i][j+1], output[i+1][j-1], output[i+1][j], +output[i+1][j+1]]
                    np.sort(value1)   # 对这九个数进行排序
                    value = value1[4]    # 中值为排序后中间这个数的正中间
                    output1[i-1][j-1] = value
        elif windowsize == 5:
            output1 = np.zeros(a.shape, np.uint8)
            for i in range(2, output.shape[0]-2):       # 求齐周围25个方格与模版进行卷积
                for j in range(2, output.shape[1]-2):
                    value1 = [output[i-2][j-2],output[i-2][j-1],output[i-2][j],output[i-2][j+1],output[i-2][j+2],output[i-1][j-2],output[i-1][j-1],output[i-1][j],output[i-1][j+1],\
                                output[i-1][j+2],output[i][j-2],output[i][j-1],output[i][j],output[i][j+1],output[i][j+2],output[i+1][j-2],output[i+1][j-1],output[i+1][j],output[i+1][j+1],\
                                output[i+1][j+2],output[i+2][j-2],output[i+2][j-1],output[i+2][j],output[i+2][j+1],output[i+2][j+2]]
                    value1.sort()   # 对这九个数进行排序
                    value = value1[12]    # 中值为排序后中间这个数的正中间
                    output1[i-2][j-2] = value   # 将计算结果填入原本位置
        else :
            print('模版大小输入错误,请输入3或5,分别代表3*3或5*5模版!')
        return output1
    

    均值滤波:

    # 均值滤波 a为要处理的图像  windowsize为采用的模版大小
    def meanflite(a, windowsize):
        output = a
        if windowsize == 3:
            window = np.ones((3, 3)) / 3 ** 2  # 生成3*3模版
            output1 = np.zeros(a.shape, np.uint8)
            for i in range(1, output.shape[0] - 1):  # 求齐周围9个方格与模版进行卷积
                for j in range(1, output.shape[1] - 1):
                    value = (output[i - 1][j - 1] * window[0][0] + output[i - 1][j] * window[0][1] + output[i - 1][j + 1] *
                             window[0][2] + \
                             output[i][j - 1] * window[1][0] + output[i][j] * window[1][1] + output[i][j + 1] * window[1][
                                 2] +\
                             output[i + 1][j - 1] * window[2][0] + output[i + 1][j] * window[2][1] + output[i + 1][j + 1] *
                             window[2][2])
                    output1[i - 1][j - 1] = value  # 将计算结果填入原本位置
        elif windowsize == 5:
            window = np.ones((5, 5)) / 5 ** 2  # 生成5*5模版
            output1 = np.zeros(a.shape, np.uint8)
            for i in range(2, output.shape[0] - 2):  # 求齐周围25个方格与模版进行卷积
                for j in range(2, output.shape[1] - 2):
                    value = (output[i - 2][j - 2] * window[0][0] + output[i - 2][j - 1] * window[0][1] + output[i - 2][j] *
                             window[0][2] + output[i - 2][j + 1] * window[0][3] + output[i - 2][j + 2] * window[0][4] + \
                             output[i - 1][j - 2] * window[1][0] + output[i - 1][j - 1] * window[1][1] + output[i - 1][j] *
                             window[1][2] + output[i - 1][j + 1] * window[1][3] + output[i - 1][j + 2] * window[1][4] + \
                             output[i][j - 2] * window[2][0] + output[i][j - 1] * window[2][1] + output[i][j] * window[2][
                                 2] + output[i][j + 1] * window[2][3] + output[i][j + 2] * window[2][4] + \
                             output[i + 1][j - 2] * window[3][0] + output[i + 1][j - 1] * window[3][1] + output[i + 1][j] *
                             window[3][2] + output[i + 1][j + 1] * window[3][3] + output[i + 1][j + 2] * window[3][4] + \
                             output[i + 2][j - 2] * window[4][0] + output[i + 2][j - 1] * window[4][1] + output[i + 2][j] *
                             window[4][2] + output[i + 2][j + 1] * window[4][3] + output[i + 2][j + 2] * window[4][4])
                    output1[i - 2][j - 2] = value  # 将计算结果填入原本位置
        else:
            print('模版大小输入错误,请输入3或5,分别代表3*3或5*5模版!')
    
        return output1
    

    高斯滤波:

    ef gaussian(im):
        im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
        b = np.array([[2, 4,  5,  2,  2],
                   [4, 9,  12, 9,  4],
                   [5, 12, 15, 12, 5],
                   [4, 9,  12, 9,  4],
                   [2, 4,  5,  4,  2]]) / 156
        kernel = np.zeros(im.shape)
        kernel[:b.shape[0], :b.shape[1]] = b
    
        fim = np.fft.fft2(im)
        fkernel = np.fft.fft2(kernel)
        fil_im = np.fft.ifft2(fim * fkernel)
    
        return abs(fil_im).astype(int)
    

    3.自编写函数的主程序

    我用到的包有:numpy、matplotlib.pyplot、cv2、random
    主程序很好理解,读入原图,再对原图调用各种函数进行处理,再以子图的形式展示出来。

    if __name__ == "__main__":
        image = cv2.imread('whl.jpg')
    
        plt.subplot(4, 2, 1)
        plt.imshow(image)
        plt.axis('off')
        plt.title('Original')
    
        salt = salt_pepper_noise(image, 0.05)
    
        plt.subplot(4, 2, 2)
        plt.imshow(salt)
        plt.axis('off')
        plt.title('salt')
    
        gauss = gasuss_noise(image)
    
        plt.subplot(4, 2, 3)
        plt.imshow(gauss)
        plt.axis('off')
        plt.title('gauss')
    
        random = addrandom_noise(image)
    
        plt.subplot(4, 2, 4)
        plt.imshow(random)
        plt.axis('off')
        plt.title('random')
    
        median = medianfliter(salt, 3)
    
        plt.subplot(4, 2, 5)
        plt.imshow(median)
        plt.axis('off')
        plt.title('median')
    
        mean = meanflite(random, 3)
    
        plt.subplot(4, 2, 6)
        plt.imshow(mean)
        plt.axis('off')
        plt.title('mean')
    
        gaussout = gaussian(gauss)
    
        plt.subplot(4, 2, 7)
        plt.imshow(gaussout)
        plt.axis('off')
        plt.title('gaussout')
    
        plt.show()
    
    

    效果展示为:
    在这里插入图片描述

    4.opencv实现噪声+滤波器

    以opencv中的函数实现各种噪声的产生和滤波器的编写更为简单,基本一行代码就可以实现调用,以下是我的整个程序:

    # 图像去噪平滑滤波
    # 使用opencv的自带函数实现,与自编写作比较
    # 产生椒盐噪声,高斯噪声等
    # 使用中值滤波,平均滤波,高斯滤波,方框滤波
    
    import numpy as np
    import cv2
    import matplotlib.pyplot as plt
    
    
    # 加噪声
    def noise(img):
        out = img
        rows, cols, chn = img.shape
        for i in range(5000):
            x = np.random.randint(0, rows)
            y = np.random.randint(0, cols)
            out[x, y, :] = 255
        return out
    
    
    if __name__ == "__main__":
        image = cv2.imread('whl.jpg')
    
        plt.subplot(3, 2, 1)
        plt.imshow(image)
        plt.axis('off')
        plt.title('Original')
    
        noise_img = noise(image)
    
        plt.subplot(3, 2, 2)
        plt.imshow(noise_img)
        plt.axis('off')
        plt.title('noise')
    
        # 均值滤波
        result1 = cv2.blur(noise_img, (5, 5))
    
        plt.subplot(3, 2, 3)
        plt.imshow(result1)
        plt.axis('off')
        plt.title('mean')
    
        # 方框滤波
        result2 = cv2.boxFilter(noise_img, -1, (5, 5), normalize=1)
    
        plt.subplot(3, 2, 4)
        plt.imshow(result2)
        plt.axis('off')
        plt.title('box')
    
        # 高斯滤波
        result3 = cv2.GaussianBlur(noise_img, (3, 3), 0)
    
        plt.subplot(3, 2, 5)
        plt.imshow(result3)
        plt.axis('off')
        plt.title('gaussian')
    
        # 中值滤波
        result4 = cv2.medianBlur(noise_img, 3)
    
        plt.subplot(3, 2, 6)
        plt.imshow(result4)
        plt.axis('off')
        plt.title('median')
    
        plt.show()
    
    

    效果展示为:
    在这里插入图片描述

    展开全文
  • 输入一幅灰度图像,分别采用邻域平均法、均值滤波法、中值滤波法对图像进行平滑去噪,在同一个窗口输出显示原始图像和3种结果图像。
  • 几种经典的图像平滑去噪算法程序

    热门讨论 2011-01-04 15:37:26
    几种经典的图像平滑去噪算法程序,其中包括小波软阈值和硬阈值处理。还有三个图像质量评价程序
  • 滤波 平滑 去噪

    千次阅读 2019-01-06 22:23:03
    3. 去噪:一种具体应用,用有噪音污染的信号猜测干净的信号。 他们的联系? 一句话给串起来: 要降噪,常用滤波来实现,平且常用低通滤波,副作用是往往使得信号平滑化。 滤波:是从信号...

    1.滤波:信号处理的一种工具,并且是最常用的工具。

    2. 平滑:形容词,顾名思义,描述信号的特性smooth;    转换到频域,能量大概率集中在低频。

      做动词 - 使平滑化:使信号变得具有平滑的属性

    3. 去噪:一种具体应用,用有噪音污染的信号猜测干净的信号。

    他们的联系?

    一句话给串起来:

    要降噪,常用滤波来实现,平且常用低通滤波,副作用是往往使得信号平滑化。

    滤波:是从信号中获取需要的频率成分或者去除掉不需要的频率成分。

    平滑:在图像处理中,平滑会使得图像模糊。对于高斯噪声有较好的滤除效果,但对于椒盐噪声则无能为力。平滑是一种时域手段(图像处理时域也叫空域,因为是二维的),在频域下就是低通滤波了。(以上可以参考冈萨雷斯的数字图像处理,讲解的比较详细。)

    噪声:影响通信质量的信号是噪声。滤波的一类问题就是滤除噪声。噪声也分很多种,研究不同噪声的特性,就是为了找出对应滤除的方法。

    联系:平滑是时域方式,用频域滤波(低通滤波)的方法也可以实现平滑所达到的效果,也就是去除高斯噪声。

     

    最近在学图像处理,说说我的看法

    我认为的滤波,他是一种处理手段,我理解的滤波包含去噪平滑,但是同时也远不止这些,滤波的含义是让待处理信号通过滤波器(系统)然后得出一个我们需要的目标信号,这个过程叫做滤波,是个很宽泛的概念。

    至于平滑,在图像处理领域,平滑是对空域提出的一种要求,而非频域,它的目的是为了使图像的边缘看上去不那么锐利,与锐化的结果相反,为了实现平滑,我们就要对图像进行滤波,可以通过空域技术进行滤波,比如采用领域均值模板运算,也可以采用频域技术,比如通过低通滤波器。

    至于去噪,它的目的是为了去除信号中不要的成分,可以是图像中的噪点,可以是通信信号中的白噪声,而去噪也需要通过滤波。

    总而言之,滤波是信号处理的最基本步骤,任何对信号的操作都可以称之为滤波,而平滑和降噪算是滤波的一种手段

    以上是我个人的理解

     

     

    作者:SuperMHP

    链接:https://www.zhihu.com/question/35073143/answer/149190299

    来源:知乎

    著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

     

    为了好好解释这个问题,只能疯狂偷老板的图了。这个问题下的各位对滤波有点误解,不要被中文名字迷惑了。

    滤波器不是为了滤除什么,而是为了提取我们感兴趣的信号分量,如果没有这个概念,就很难理解除了去噪/去干扰之外的其他应用。

    滤波和平滑是滤波器的两种状态,去噪是滤波器的应用。

    1.滤波器分类

    滤波器是想通过某种方式获取感兴趣的信号分量,传统的滤波器(比如低通,高通等)是通过某些确定性规则/门限来获取,自适应滤波(比如维纳,卡尔曼滤波)是在通过统计规则来获取。传统滤波器我们不讨论,没有那么多操作。

     

     

     

    2.自适应滤波器的三种状态: 滤波,预测,平滑

    平滑只是滤波器的一种工作状态,如果感兴趣的信号部分在过去,现在,未来,那么滤波器是在分别做平滑,滤波,预测。

     

    根据时间序号 K和信号样本集N ,通过统计规则实现的自适应滤波器有三种效果:

     

    当 N<k ,这是在预测未来的状态,我们在做预测估计(predicted estimate),也就是预测

    当 N=k ,这是在用所有过去的观测值和当前的观测值,在估计当前的状态,这是在做滤波估计(filtered estimate),也就是滤波

    当N>k ,这是在估计过去的状态,这是在做平滑估计(smoothed estimate),也就是平滑

    偷2张我老板在讲维纳滤波时的ppt:

    需要注意的是,预测和滤波可以在real time 内完成,平滑不可能在real time 内完成。

     

    3. 滤波器的应用

     

    滤波器的输出是由我们感兴趣的信号部分所决定的,比如在

    里,我们感兴趣的是信号主成分

    ,那么它就是在去噪。

     

    如果维纳滤波用于时序信号去噪,那么如下图,它通过对噪声的补偿,来从被噪声污染的信号v(n)中复现信号s(n)。

    如果它应用于空域滤波去噪,假设有两天线信号被噪声污染,是这样的,

    或者,比如在MIMO信道模型Y=HX+W 里,我们感兴趣的部分是信道 H ,那么它就是在做信道估计。当然它也可以做去卷积和空间滤波。

     

    4. 时域和空域

    因果性(Causality)是real time 信号处理最重要的问题之一,然而如果滤波器输入是空域参数(图像处理,阵列响应),而不是时域参数(通信信号),因果性可以不必考虑。

    展开全文
  • 我写的一些关于平滑去噪、锐化、灰度、伪彩色操作的C#代码。大家可以参考一下
  • 前言:最近研究汽车碰撞的加速度信号,在信号的采集过程中难免遇到噪音,导致信号偏差,为了更好的反映系统情况,故常需要信号去噪,本文分享一些常用信号平滑去噪的方法。关键字:信号;去噪;Matlab信号在实际测量...

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

    关键字:信号;去噪;Matlab


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


    1、移动平均法

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

    b06caf9b042a4633a60e8e470faa9e4f.png

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

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

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

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

    5a319881f41a3a1b769066e9a254b52d.png


    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)

    8feffd7e7700feeee5f7f6cb78d370a9.png


    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)

    65df897f55cb57a617a5de49b86c1e6f.png

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

    c6b4fd4713048b3c28e9d598ab1e3237.png

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

    y = filter(b,a,x)

    其中b和a的定义为:

    7d24015323f59ad1ef8eceace359e0a1.png

    其中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。

    3d93dddd5e854cbdbde0741af41ebc00.png


    5、移动平均的幅频响应

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

    以中心窗口为例,

    b58d117156ceb317032a8a3e7f19a550.png

    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)

    9f6b5a3b67916816fb6a92ce0cb5a6bc.png

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

    e6bc9f37ed4a66ced69ebc3ff209b335.png

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

    反之,以美国近期的疫情为例,疫情的采样频率为1天一采样,而且显示出很强的7日一周期的特性(病毒也要过周末bd27916862ff985cae2d96a9ae0766b3.png)。所以,为了消除这个归一化频率为1/7的噪声影响,采样7点的滑动平均滤波。可以看到所有以7天为一变化的信号分量全部被消除掉了。(下面这个图经常被引用,但是很少有人思考为什么用7天平均方法来平滑数据。)

    9ffafe01d44867b084f1816bf3b458ea.png

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

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


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

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

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

    526f5d670807c53730c22227ba56752a.png

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

    59c695d0085f02338cae4c6d2b34389c.png468c2ac8625c2f7f807ccc984fb74a3f.png22f76fa7aabfc78d27582ff3e1ba97f1.png

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

    展开全文
  • Savitzky-Golay平滑去噪

    千次阅读 2019-12-04 17:36:27
    推导及使用如下博客讲的很清楚 参考如下博客: https://blog.csdn.net/qq_20823641/article/details/51537461 https://blog.csdn.net/schuffel/article/details/85239770
  • MATLAB对图像进行直方图统计、灰度变换、直方图均衡、平滑去噪 本文基于matlab实现直方图统计以及均衡化,并做简单的灰度变换和去噪。所有的实现基于对像素点的操作,并不调用库函数。 线性拉伸代码 a9=imread('D:a2...
  • [原创]Java调用MATLAB实现三点平滑去噪并绘图 http://www.matlabsky.com/thread-24805-1-1.html 【题外话】周末的晚上,窗外下着雨,感谢那个她,使我觉得屏幕上的一行行代码是那么地优美、动人...... ...
  • C# 多种方法对图像进行平滑去噪

    热门讨论 2011-04-30 20:27:19
    中值和中值滤波 灰度形态学滤波 小波变换去噪 高斯低通滤波 统计滤波
  • 实验五 图像增强 平滑去噪 分别在灰度图像中加入一定量的高斯噪声和椒盐噪声,然后采用3×3的均值滤波器和3×3中值滤波器分别对噪声图像进行处理,给出两种处理方法的峰值信噪比(PSNR);仿效“中值滤波”的方法,对...
  • CUDA:并行计算实现图像平滑去噪

    千次阅读 2018-10-18 22:24:05
         

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 723
精华内容 289
关键字:

平滑去噪