精华内容
下载资源
问答
  • 上波找点
    千次阅读 多人点赞
    2021-05-04 11:43:36

    【小波分析】一、小波分析入门基础介绍

    引入

    从线性代数谈起

    我们在线性代数当中学过点、向量、矩阵,线性变换等一些概念,我们一起回顾一下,从而引出函数空间,为小波分析的介绍做准备。

    我们知道在一个线性空间下,比如说三维的线性空间,有一个坐标系,我们喜欢叫 x x x 轴、 y y y 轴、 z z z 轴,在这个三维空间中的每一个点,都有一个坐标表达(点的表示),比如说 ( x , y , z ) (x,y,z) (x,y,z),分别对应这这个点在每个坐标系上面的分量,或者说是这个点在分别在三个坐标轴下的投影。

    线性变换 A \mathscr{A} A 对应这一个矩阵 A A A,得到新的向量 ( x ′ , y ′ , z ′ ) (x',y',z') (x,y,z) 它作用于一个点(坐标表达为一个向量),它意味着什么呢?我们可以理解为,在原来的坐标系下,对这个点进行了一定的操作,变到了新的地方。一个更好地理解是,我们找到了一个新的坐标系,三轴为 x ′ x' x 轴、 y ′ y' y 轴、 z ′ z' z 轴,在这个新的坐标系下,原来这个点被表示为了 ( x ′ , y ′ , z ′ ) (x',y',z') (x,y,z)。也就是说,一个线性变换,其实对应的是一个坐标系的变换,线性变换作用于点,前后对应着点在不同坐标系下的表示。

    举个例子,某个点在坐标系 O ( x , y , z ) O(x,y,z) O(x,y,z) 下的表示为 ( x , y , z ) (x,y,z) (x,y,z),那么我们再给定一个新的坐标系 O ( x ′ , y ′ , z ′ ) O(x',y',z') O(x,y,z),不妨假设原点相同,它在新的坐标系下得坐标表示是什么呢?假设原空间的一组基向量为 ( e 1 , e 2 , e 3 ) (e_1,e_2,e_3) (e1,e2,e3),变换后的一组基向量为 ( e 1 ′ , e 2 ′ , e 3 ′ ) (e_1',e_2',e_3') (e1,e2,e3),这里的 e i e_i ei 表示列向量。若把原空间的基向量写为后面那个空间的基向量的线性组合有,
    ( e 1 , e 2 , e 3 ) = ( e 1 ′ , e 2 ′ , e 3 ′ ) A (e_1,e_2,e_3) = (e_1',e_2',e_3') A (e1,e2,e3)=(e1,e2,e3)A
    则有
    ( x y z ) = A ( x ′ y ′ z ′ ) \left( \begin{array}{l} x\\y\\z \end{array}\right) = A\left( \begin{array} {l} x'\\y'\\z' \end{array}\right) xyz=Axyz.

    这里考虑的是两个坐标系原点的相同的情况。若两个若两个坐标系的坐标原点不同,再加个坐标平移量 O − O ′ O-O' OO即可。当然,也可以把坐标写成 4 维的,矩阵也扩充成 4x4 的,然后把平移量融合到加出来的位置里面,在机器人运动学里面,经常这么搞 。

    从上面所述,我们知道了一个线性变换其实对应着一个坐标系的变换。不同矩阵的性质,其实也就对应着变换的不同性质。比如说正交矩阵(复数域下我们成为酉矩阵),其实就是把一个基向量是单位长度的标准正交系变成了另外一个标准正交系。换个角度看,它也对应着向量的旋转。其他矩阵的性质也有类似的集合意义。

    线性代数这门学科本质上回答了两个问题:点的描述和线性变换的刻画。

    到高等数学

    很自然地,我们现在思考,如果把线性空间中的点和线性变换的刻画,往函数空间去推广一下会怎么样。假设现在空间中的每一个点刻画的是一个函数,线性空间的每个基,表示函数空间的基函数,线性变换,把从一组已知的基函数,变换到了另外一组基函数,把一组基函数张成的空间,变成了另外一组基函数张成的空间。这,就自然而然实现了从代数到分析的一个跨越。

    所谓的变换,譬如说,傅里叶变换,小波变换,无非就是说找到一个变换,或者说找到了一个组新的基,把原空间里面的函数,在新的一组基下进行表示,进而,我们来研究同一个函数,在两个空间下得性质有什么样的联系和区别。这一组新的基,不是随随便便找的,它的寻找必然是有利于我们的某种需求,比如说更少维度的表达从而达到存储的压缩,等等。

    历史和发展

    傅里叶分析

    傅里叶是一名伟大的法国科学家。拿破仑当政期间,法国科学院在全欧洲出了一个竞赛题,去描述热在介质中的传播规律。很多人去做了这个热试验,即拿一个管子,在一端加热,然后去测试管子每一点上的温度。傅里叶不仅做了实验,还给出了一个解决方式。他在 1807 年提出,解能写出这样的形式。
    f ( t ) = ∑ k = − ∞ + ∞ α k e i k t f(t)=\sum_{k=-\infty}^{+\infty} \alpha_{k} e^{i k t} f(t)=k=+αkeikt
    可以看得出,这里的基函数,不是别的,是这种三角函数,或者说 e 指数函数的一系列伸缩。简单地说,傅里叶给出了一个基,它能够解决热传导问题。现在很多相关的名词描述什么都非常成熟了,这都是后来的人给加的。原始的 想法可以看傅里叶本人的专著《热的解析理论》。

    傅里叶把他的想法整理了一下,寄到了法国科学院去,很遗憾,没有任何反应。当时有两派意见,一派认为,这玩意儿太简单了吧,几乎全是文字,没什么意思。另一派的数学家认为,这个吧,一个证明没有,净瞎扯。随手一丢,装箱子里去了。

    后来傅里叶当官去了。之后拿破仑政府倒台了,但是傅里叶没受到太大牵连。那个时候,他对他的工作进一步进行了补充,即有了傅里叶变换。

    傅里叶变换,就不要求函数是周期的了。它可以写成,
    f ( t ) = 1 2 π ∫ − ∞ + ∞ F ( ω ) e i ω t d t f(t) = \frac{1}{2 \pi} \int_{-\infty}^{+\infty} F(\omega) e^{i \omega t} d t f(t)=2π1+F(ω)eiωtdt

    上下两个式子何其相似。只不过把离散的 k k k变成了频率 ω \omega ω,离散求和变成了积分。

    这时候,傅里叶又把他这个提交到了法国科学院,法国科学院,还是分成了两派。一派觉得他解决得很好,相当漂亮的变大。另外一派还是觉得没有数学证明,这个不行啊。后来,是确认傅里叶获奖了,但是不给钱。后来,别的学科就把这个公式拿去用,发现 work 得很好,逐渐就流行开了。 很多东西,因为侧重点不一样,所以在很多人看起来,评价就不一样。咱们学数学的,似乎更倾向于要有一个漂亮的解析证明,工程的人,似乎关心它 work 不 work。比如说,神经网络,很多时候 work 得很好,但是可解释性差,所以数学的,比如说学统计的,就非常抵制它,但是工程上的人又非常地喜欢用它用到具体的问题上去算一些东西。

    哈尔基

    这个傅里叶分析,就这样发展了 100 年,直到 20 世纪初期,一个叫哈尔的哥们,他说,傅里叶这东西,也没啥了不起的,无非不就是说可以把一个函数展开为一些基函数的线性组合么,我随便就能找出这样的基函数来。

    在 1908 年的国际数学家会议上,哈尔就提出了这样一个函数,
    h ( t ) = { 1 0 ≤ x < 1 2 − 1 1 2 ≤ x < 1 0 其他  h(t)=\left\{\begin{array}{ll} 1 & 0\leq x<\frac{1}{2} \\ -1& \frac{1}{2} \leq x<1 \\ 0 & \text {其他 } \end{array}\right. h(t)=1100x<2121x<1其他 

    h j , k ( t ) = 2 j 2 h ( 2 j t − k ) h_{j, k}(t)=2^{\frac{j}{2}} h\left(2^{j} t-k\right) hj,k(t)=22jh(2jtk)
    它构成了 L 2 ( R ) L^{2}(\mathbb{R}) L2(R)空间的一组标准正交基。可以看得出,这个 h j , k h_{j,k} hj,k是对于 h h h的伸缩和平移, 2 j 2^j 2j表示缩了 2 j 2^j 2j倍, − k -k k表示往右拉了 k k k个格子。可以举例的。

    于是,哈尔说,任意函数,也能够写成这组基的线性组合。
    f ( t ) = ∑ j ∑ k α j k h j , k ( t ) f(t)=\sum_{j} \sum_{k} \alpha_{j k} h_{j, k}(t) f(t)=jkαjkhj,k(t)
    这个结果,也很遗憾,在数学家会议上没有引起任何注意。

    1910年,这个工作在 siam 上发表了。就是说,傅里叶那个玩意儿,你能找一个基,函数能表示成他们任意的线性组合,我也能。

    加窗傅里叶变换

    再后来,就到了加窗傅里叶变换,窗,意味着窗口。Gabor 在 1946 年他的博士论文中,提到了一组新的基函数,给原来的三角基函数,加了一个窗口限制,即

    g ( t − t i ) e i ω t g\left(t-t_{i}\right) e^{i \omega t} g(tti)eiωt

    更清楚一点,写成频率分量的表达,是
    F ( t 0 , w ) = ∫ − ∞ + ∞ f ( t ) g ( t − t 0 ) e − i ω t d t F\left({t_0}, w\right)=\int_{-\infty}^{+\infty} f(t) g\left(t-t_{0}\right) e^{-i \omega t} d t F(t0,w)=+f(t)g(tt0)eiωtdt
    这个表示一个信号 f ( t ) f(t) f(t),在时间点 t = t 0 t=t_0 t=t0的附近,频率为 ω \omega ω的频率分量的大小。这里的附近的范围,是由窗函数 g g g确定的。

    Gabor 指出,把所有时间点附近的频率分量拼在一起,就是这个信号在这个频率上的分量表示,写成公式就是,
    F ( w ) = ∫ − ∞ + ∞ F ( t 0 , w ) d t 0 F(w)=\int_{-\infty}^{+\infty} F\left(t_{0}, w\right) d t_0 F(w)=+F(t0,w)dt0

    那么,由傅里叶变换,一个信号,就可以表示为,
    f ( t ) = ∫ − ∞ + ∞ F ( ω ) e i ω t d ω = ∫ − ∞ + ∞ ∫ − ∞ + ∞ F ( t 0 , ω ) e i ω t d ω d t 0 f(t)=\int_{-\infty}^{+\infty} F(\omega) e^{i \omega t} d \omega = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} F\left(t_{0}, \omega\right) e^{i\omega t} d \omega d t_{0} f(t)=+F(ω)eiωtdω=++F(t0,ω)eiωtdωdt0

    这个东西一提出,学术界很快就接受了它。后来就出现了上百种窗函数,譬如高斯窗,矩形窗等。

    小波分析

    这个再往后,就到了19世纪80年代左右, Morlet ,一个搞地质学的兄弟,也就是提出小波变换这个概念的这个哥们。 Morlet 在处理一些地质勘探数据的时候,为了自己的方便,就开始琢磨了。

    他说这个加窗的傅里叶变换,有点问题啊,就是说,窗函数宽度和 ω \omega ω 没关系。 e i ω t e^{i \omega t} eiωt 是周期的正余弦函数, ω \omega ω越大,频率越高,抖得越厉害,每个周期长度就越小,这时候,截取的窗口就可以小一些。反之亦然。现有的窗函数,对于固定时间点附近的这个范围,似乎和 ω \omega ω没有关系,这个不太合理啊。所以,心里萌生了以一种 “高频窄窗低频宽窗” 的思想。

    基本想法是,希望能找到一个函数,在显著的时间内,有明显的波形,在这个显著的时间外,当然是希望没有。做不到没有,当然是希望越来越小,直至没有,在外面迅速地衰减。Morlet 提出了,

    1 ∣ a ∣ ψ ( t − b a ) \frac{1}{\sqrt{|a|}} \psi\left(\frac{t-b}{a}\right) a 1ψ(atb)

    这里的 ψ \psi ψ是满足一定要求的函数,后面还会提到。 b b b 是位移参数,左加右减。 a a a 是尺度参数,表示把这个波形放大了多少倍,体现的是快变信息和慢变信息。外面开根号 a a a分之一,是为了能量的守恒,里面除 a a a,外面除的是 a \sqrt a a ,因为能量是模方的表达,平方一下,刚好消掉。

    同样,1978 年 Morlet 在国际会议上宣读自己的想法,也没激起什么浪花。1982 年刊登他的成果,很长时间内无人问津。

    但是,这个时候,小波这个工具,已然悄悄萌芽。

    参考和致谢

    《Nonlinear approximation》 Ronald A. Devore 剑桥大学出版社

    《小波十讲》Ingrid Daubechies 李建平等译 国防工业出版社

    《小波与算子》(第一卷 ) Y.Meyer 尤众译 世界图书出版社

    《小波变换与分数傅里叶变换:理论及应用》 冉启文 哈尔滨工业大学出版社

    “小波理论及应用” 冉启文 【哈尔滨工业大学】

    相关的 Wikipedia 和网络资源

    更多相关内容
  • 去噪基本概念

    千次阅读 2020-12-30 09:26:41
    一、前言在现实生活和工作中,噪声无处不在,在许多领域中,如天文、医学图像和计算机视觉方面收集到的数据常常是含有噪声的。...小分析是近年来发展起来的一种新的信号处理工具,这种方法源于傅立叶分析,小(...

    一、前言

    在现实生活和工作中,噪声无处不在,在许多领域中,如天文、医学图像和计算机视觉方面收集到的数据常常是含有噪声的。噪声可能来自获取数据的过程,也可能来自环境影响。由于种种原因,总会存在噪声,噪声的存在往往会掩盖信号本身所要表现的信息,所以在实际的信号处理中,常常需要对信号进行预处理,而预处理最主要的一个步骤就是降噪。

    小波分析是近年来发展起来的一种新的信号处理工具,这种方法源于傅立叶分析,小波(wavelet),即小区域的波,仅仅在非常有限的一段区间有非零值,而不是像正弦波和余弦波那样无始无终。小波可以沿时间轴前后平移,也可按比例伸展和压缩以获取低频和高频小波,构造好的小波函数可以用于滤波或压缩信号,从而可以提取出已含噪声信号中的有用信号。

    二、小波去噪的原理

    Donoho提出的小波阀值去噪的基本思想是将信号通过小波变换(采用Mallat算法)后,信号产生的小波系数含有信号的重要信息,将信号经小波分解后小波系数较大,噪声的小波系数较小,并且噪声的小波系数要小于信号的小波系数,通过选取一个合适的阀值,大于阀值的小波系数被认为是有信号产生的,应予以保留,小于阀值的则认为是噪声产生的,置为零从而达到去噪的目的。

    从信号学的角度看 ,小波去噪是一个信号滤波的问题。尽管在很大程度上小波去噪可以看成是低通滤波 ,但由于在去噪后 ,还能成功地保留信号特征 ,所以在这一点上又优于传统的低通滤波器。由此可见 ,小波去噪实际上是特征提取和低通滤波的综合 ,其流程图如下所示:

    一个含噪的模型可以表示如下:

    其中 ,f( k)为有用信号,s(k)为含噪声信号,e(k)为噪声,ε为噪声系数的标准偏差。

    假设,e(k)为高斯白噪声,通常情况下有用信号表现为低频部分或是一些比较平稳的信号,而噪声信号则表现为高频的信号,我们对 s(k)信号进行小波分解的时候,则噪声部分通常包含在HL、LH、HH中,如下图所示,只要对HL、LH、HH作相应的小波系数处理,然后对信号进行重构即可以达到消噪的目的。

    我们可以看到,小波去噪的原理是比较简单类,类似以往我们常见的低通滤波器的方法,但是由于小波去找保留了特征提取的部分,所以性能上是优于传统的去噪方法的。

    三、小波去噪的基本方法

    一般来说, 一维信号的降噪过程可以分为 3个步骤

    信号的小波分解。选择一个小波并确定一个小波分解的层次N,然后对信号进行N层小波分解计算。

    小波分解高频系数的阈值量化。对第1层到第N层的每一层高频系数(三个方向), 选择一个阈值进行阈值量化处理.

    这一步是最关键的一步,主要体现在阈值的选择与量化处理的过程,在每层阈值的选择上matlab提供了很多自适应的方法, 这里不一一介绍,量化处理方法主要有硬阈值量化与软阈值量化。下图是二者的区别:

    上面左图是硬阈值量化,右图是软阈值量化。采用两种不同的方法,达到的效果是,硬阈值方法可以很好地保留信号边缘等局部特征,软阈值处理相对要平滑,但会造成边缘模糊等失真现象。

    信号的小波重构。根据小波分解的第 N层的低频系数和经过量化处理后的第1层到第N 层的高频系数,进行信号的小波重构。

    小波阀值去噪的基本问题包括三个方面:小波基的选择,阀值的选择,阀值函数的选择。

    (1)小波基的选择:通常我们希望所选取的小波满足以下条件:正交性、高消失矩、紧支性、对称性或反对称性。但事实上具有上述性质的小波是不可能存在的,因为小波是对称或反对称的只有Haar小波,并且高消失矩与紧支性是一对矛盾,所以在应用的时候一般选取具有紧支的小波以及根据信号的特征来选取较为合适的小波。

    (2)阀值的选择:直接影响去噪效果的一个重要因素就是阀值的选取,不同的阀值选取将有不同的去噪效果。目前主要有通用阀值(VisuShrink)、SureShrink阀值、Minimax阀值、BayesShrink阀值等。

    (3)阀值函数的选择:阀值函数是修正小波系数的规则,不同的反之函数体现了不同的处理小波系数的策略。最常用的阀值函数有两种:一种是硬阀值函数,另一种是软阀值函数。还有一种介于软、硬阀值函数之间的Garrote函数。

    另外,对于去噪效果好坏的评价,常用信号的信噪比(SNR)与估计信号同原始信号的均方根误差(RMSE)来判断。

    参考:

    小波阀值去噪法基础  http://blog.sina.com.cn/s/blog_4d7c97a00101cib3.html

    在python中使用小波分析进行阈值去噪声,使用pywt.threshold函数

    #coding=gbk

    #使用小波分析进行阈值去噪声,使用pywt.threshold

    import pywt

    import numpy as np

    import pandas as pd

    import matplotlib.pyplot as plt

    import math

    data = np.linspace(1, 10, 10)

    print(data)

    # [ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]

    # pywt.threshold(data, value, mode, substitute) mode 模式有4种,soft, hard, greater, less; substitute是替换值

    data_soft = pywt.threshold(data=data, value=6, mode='soft', substitute=12)

    print(data_soft)

    # [12. 12. 12. 12. 12. 0. 1. 2. 3. 4.] 将小于6 的值设置为12, 大于等于6 的值全部减去6

    data_hard = pywt.threshold(data=data, value=6, mode='hard', substitute=12)

    print(data_hard)

    # [12. 12. 12. 12. 12. 6. 7. 8. 9. 10.] 将小于6 的值设置为12, 其余的值不变

    data_greater = pywt.threshold(data, 6, 'greater', 12)

    print(data_greater)

    # [12. 12. 12. 12. 12. 6. 7. 8. 9. 10.] 将小于6 的值设置为12,大于等于阈值的值不变化

    data_less = pywt.threshold(data, 6, 'less', 12)

    print(data_less)

    # [ 1. 2. 3. 4. 5. 6. 12. 12. 12. 12.] 将大于6 的值设置为12, 小于等于阈值的值不变

    1、小波变换的理解

    傅里叶变换——短时傅里叶变换——小波变换。

    参考文献:以下两篇参考资料讲述得十分清楚,有助于理解小波变换。

    但具体的数学角度阐述,请参考其他资料。

    (1)知乎专栏:形象易懂讲解算法I——小波变换

    https://zhuanlan.zhihu.com/p/22450818

    (2)知乎专栏:傅里叶分析之掐死教程。

    https://zhuanlan.zhihu.com/p/19763358

    2、小波包分解

    小波包是为了克服小波分解在高频段的频率分辨率较差,而在低频段的时间分辨率较差的问题的基础上而提出的。

    它是一种更精细的信号分析的方法,提高了信号的时域分辨率。

    下面是两者的对比图:

    3、能量谱

    基于小波包分解提取多尺度空间能量特征的原理是把不同分解尺度上的信号能量求解出来,将这些能量值按尺度顺序排列成特征向量供识别使用。

    20180510补充更新:具体计算公式如下所示,本文中未使用重构后的系数进行能量值计算,直接使用小波包分解后的系数,参考文献《基于小波包能量特征的滚动轴承故障监测方法 》。

    4、Matlab代码

    给出两部分代码,写成两个函数。一个是小波包分解与重构,另一个是能量谱函数。

    下载地址:https://download.csdn.net/download/ckzhb/10030651

    代码名称:wavelet_packetdecomposition_reconstruct

    function wpt= wavelet_packetdecomposition_reconstruct( x,n,wpname )

    %% 对信号进行小波包分解,得到节点的小波包系数。然后对每个节点系数进行重构。

    % Decompose x at depth n with wpname wavelet packets.using Shannon entropy.

    %

    %  x-input signal,列向量。

    %  n-the number of decomposition layers

    %  wpname-a particular wavelet.type:string.

    %

    %Author hubery_zhang

    %Date 20170714

    %%

    wpt=wpdec(x,n,wpname);

    % Plot wavelet packet tree (binary tree)

    plot(wpt)

    %% wavelet packet coefficients.default:use the front 4.

    cfs0=wpcoef(wpt,[n 0]);

    cfs1=wpcoef(wpt,[n 1]);

    cfs2=wpcoef(wpt,[n 2]);

    cfs3=wpcoef(wpt,[n 3]);

    figure;

    subplot(5,1,1);

    plot(x);

    title('原始信号');

    subplot(5,1,2);

    plot(cfs0);

    title(['结点 ',num2str(n) '  1',' 系数'])

    subplot(5,1,3);

    plot(cfs1);

    title(['结点 ',num2str(n) '  2',' 系数'])

    subplot(5,1,4);

    plot(cfs2);

    title(['结点 ',num2str(n) '  3',' 系数'])

    subplot(5,1,5);

    plot(cfs3);

    title(['结点 ',num2str(n) '  4',' 系数'])

    %% reconstruct wavelet packet coefficients.

    rex0=wprcoef(wpt,[n 0]);

    rex1=wprcoef(wpt,[n 1]);

    rex2=wprcoef(wpt,[n 2]);

    rex3=wprcoef(wpt,[n 3]);

    figure;

    subplot(5,1,1);

    plot(x);

    title('原始信号');

    subplot(5,1,2);

    plot(rex0);

    title(['重构结点 ',num2str(n) '  1',' 系数'])

    subplot(5,1,3);

    plot(rex1);

    title(['重构结点 ',num2str(n) '  2',' 系数'])

    subplot(5,1,4);

    plot(rex2);

    title(['重构结点 ',num2str(n) '  3',' 系数'])

    subplot(5,1,5);

    plot(rex3);

    title(['重构结点 ',num2str(n) '  4',' 系数'])

    end

    代码名称:wavelet_energy_spectrum

    function E = wavelet_energy_spectrum( wpt,n )

    %% 计算每一层每一个节点的能量

    %  wpt-wavelet packet tree

    %  n-第n层能量

    %

    % Author hubery_zhang

    % Date  20170714

    %%

    % 求第n层第i个节点的系数

    E(1:2^n )=0;

    for i=1:2^n

    E(i) = norm(wpcoef(wpt,[n,i-1]),2)^2; %20180604更新 原代码:E(i) = norm(wpcoef(wpt,[n,i-1]),2)

    end

    %求每个节点的概率

    E_total=sum(E);

    for i=1:2^n

    p_node(i)= 100*E(i)/E_total;

    end

    % E = wenergy(wpt); only get the last layer

    figure;

    x=1:2^n;

    bar(x,p_node);

    title(['第',num2str(n),'层']);

    axis([0 2^n 0 100]);

    xlabel('结点');

    ylabel('能量百分比/%');

    for j=1:2^n

    text(x(j),p_node(i),num2str(p_node(j),'%0.2f'),...

    'HorizontalAlignment','center',...

    'VerticalAlignment','bottom')

    end

    end

    一维、二维离散卷积的计算方法:

    一、定义

    离散信号f(n),g(n)的定义如下:

    N-----为信号f(n)的长度

    s(n)----为卷积结果序列,长度为len(f(n))+len(g(n))-1

    以3个元素的信号为例:

    f(n) = [1 2 3]; g(n) = [2 3 1];

    s(0) = f(0)g(0-0) + f(1)g(0-1)+f(2)g(0-2) = 1*2 + 2*0 + 3*0 =2

    s(1) = f(0)g(1-0) + f(1)g(1-1) + f(2)g(1-2) = 1*3 + 2*2 + 3*0 = 7

    s(2) = f(0)g(2-0) + f(1)g(2-1) + f(2)g(2-2) =1*1 + 2*3 + 3*2=13

    s(3) = f(0)g(3-0) + f(1)g(3-1) + f(2)g(3-2) =1*0 + 2*1 + 3*3=11

    s(4) = f(0)g(4-0) + f(1)g(4-1) + f(2)g(4-2) =1*0 + 2*0 + 3*1=3

    最终结果为:

    s(n) = [2 7 13 11 3]

    上述计算图示如下:

    在数学里我们知道f(-x)的图像是f(x)对y轴的反转

    g(-m)就是把g(m)的序列反转,g(n-m)的意义是把g(-m)平移的n点:

    如上图g(m)在信号处理中通常叫做滤波器或掩码,卷积相当于掩码g(m)反转后在信号f(n)上平移求和。Matlab计算卷积的函数为conv,

    >> A = [1 2 3];

    B = [2,3,1];

    convD = conv(A,B)

    convD =

    2     7    13    11     3

    相应的二维卷积定义如下:

    展开全文
  • 一篇分享了一些均值滤波相关的算法,均值滤波作为一种线性滤波器,在滤除噪声的同时也会导致边缘模糊问题。而且均值滤波对高斯噪声的效果很好,但是对于椒盐噪声的效果就很一般。但是中值滤波作为一种顺序滤波器,...

    上一篇分享了一些均值滤波相关的算法,均值滤波作为一种线性滤波器,在滤除噪声的同时也会导致边缘模糊问题。而且均值滤波对高斯噪声的效果很好,但是对于椒盐噪声的效果就很一般。但是中值滤波作为一种顺序滤波器,对于椒盐噪声的效果很好,而且保边能力很强,所以这一篇主要讨论一下中值相关的算法。

    中值滤波

    算法原理

    中值滤波很好理解,均值滤波就是在一个小窗口中求均值来取代当前像素值,而中值滤波就是通过求小窗口中的中位值来取代当前位置的方式来滤波。

    在这里插入图片描述

    如图绿色窗口就是当前的滤波窗口,在一个3X3的邻域窗口中进行滤波。那么中值滤波做的就是:

    1. 对这个邻域中的像素值进行排序:[45, 50, 52, 60, 75, 80, 90, 200, 255];
    2. 从排序后的数据中找出中位值:75;
    3. 用中位值取代当前位置的像素值,所以就可以得到右侧的滤波后的数据。

    代码实现

    %% --------------------------------
    %% author:wtzhu
    %% email:wtzhu_13@163.com
    %% date: 20210306
    %% fuction: 中值滤波
    %% --------------------------------
    clear;
    clc;
    close all;
    img = imread('./images/test_pattern_blurring_orig.tif');
    [m, n] = size(img);
    figure;subplot(221);imshow(img);title('original image');
    
    %% 运算的时候需要对边缘进行扩展
    % 需要特殊处理四周最外圈的行和列,本算法中将其向外扩展一圈,用最外圈的值填充
    headRowMat = img(1,:);%取f的第1行
    tailRowMat = img(m,:);%取f的第m行
    % 行扩展后,列扩展时需要注意四个角需要单独扩展进去,不然就成了十字架形的
    headColumnMat = [img(1,1), img(:,1)', img(m,1)];
    tailColumnMat = [img(1,n), img(:,n)', img(m,n)];
    expandImage = [headRowMat; img; tailRowMat];
    expandImage = [headColumnMat; expandImage'; tailColumnMat];
    expandImage = uint8(expandImage');
    subplot(222);imshow(expandImage);title('expand image');
    
    newImg = zeros(m, n);
    for i =2: m+1
        for j =2: n+1
           imgRoi = [expandImage(i-1, j-1) expandImage(i-1, j) expandImage(i-1, j+1) ...
                     expandImage(i  , j-1) expandImage(i  , j) expandImage(i  , j+1) ...
                     expandImage(i+1, j-1) expandImage(i+1, j) expandImage(i+1, j+1)];
           orderedList = sort(imgRoi);
           sizeRoi = size(imgRoi);
           newImg(i-1, j-1) = orderedList((sizeRoi(2)+1)/2);
        end
    end
    newImg = uint8(newImg);
    subplot(223);imshow(newImg);title('new image');
    subplot(224);imshow(uint8(newImg-img));title('newImg-img');
    

    在这里插入图片描述

    从滤波效果图中可以看出去噪能力还可以,百块中的噪点去除了很多,而且边缘信息保留得也很好。

    多级中值滤波

    算法原理

    简单的中值滤波器滤波效果有限,于是就有人提出了将多个中值滤波进行多级级联实现更好的滤波效果。

    在这里插入图片描述

    如图就是一种多级级联的方式,先在窗口中定义一个’+'和’X’形的窗口,然后分别求出这两个窗口的中位值,然后结合当前窗口的中心点就有3个候选值,再从这三个值中求出一个中位值作为滤波后的结果。

    这种方式也可以直接应用到RAW图中做BNR,需要修改的就是窗口设置为5X5,然后在做滤波的时候需要区分G和RB通道,因为这个前面已经讲过,RAW图中的RGB分布是不均匀的,G占50%,R和B各占25%。
    在这里插入图片描述

    左侧就是针对G通道的滤波器,右侧是R和B通道的滤波器,都是定义了一个’+'和’X’形的窗口,不同的只是取的点的位置不同。

    代码实现

    %% --------------------------------
    %% author:wtzhu
    %% email:wtzhu_13@163.com
    %% date: 20211023
    %% fuction: multistage median filters
    %% --------------------------------
    
    close all;
    clear all;
    clc
    img = imread('./images/lena.bmp');
    I = double(img);
    % I = double(imresize(img, [64, 64]));
    figure();
    imshow(uint8(I));
    title('org file');
    
    I_noise = I + 10 * randn(size(I));
    figure();
    imshow(uint8(I_noise));
    title('noise file');
    [m,n] = size(I_noise);
    
    DenoisedImg = zeros(m,n);
    PaddedImg = padarray(I,[1, 1],'symmetric','both');
    
    tic
    for i = 1: m
        for j = 1: n
            roi = PaddedImg(i:i+2, j:j+2);
            % first stage
            median_HV = median([roi(1,2), roi(2,1), roi(2,2), roi(2,3), roi(3,2)]);
            median_diag = median([roi(1,1), roi(1,3), roi(2,2), roi(3,1), roi(3,3)]);
            % second stage
            DenoisedImg(i, j) = median([median_HV, roi(2,2), median_diag]);
        end
    end
    figure();
    imshow(uint8(DenoisedImg));
    title('denoise file');
    toc
    
    b = medfilt2(I_noise,[3,3]);
    figure();
    imshow(uint8(b));
    title('median filter of matlab denoise file');
    
    

    多级中值混合滤波

    算法原理

    前面介绍过中值滤波和均值滤波各自有各自的优点和缺点,所以就可以考虑将两者结合起来,相互弥补实现更好的滤波效果,于是就有人提出了这种多级混合滤波的方式。

    在这里插入图片描述

    算法流程:

    1. 求出竖直方向相邻三个点的均值和水平方向相邻三个点的均值,再结合当前点,用这三个点再求一个中位值;
    2. 求出45°和135°方向上的均值,然后结合当前点求出一个中位值;
    3. 两个中位值结合当前点组成新的数组,最后求一个中位值作为当前点的值完成滤波。

    这样就巧妙的将均值滤波和中位值滤波结合了。然后如果应用在RAW图上,只需要对滤波器稍作改善即可。

    在这里插入图片描述

    代码实现

    %% --------------------------------
    %% author:wtzhu
    %% email:wtzhu_13@163.com
    %% date: 20211023
    %% fuction: multistage median hybird filters
    %% --------------------------------
    
    close all;
    clear all;
    clc
    img = imread('./images/lena.bmp');
    I = double(img);
    figure();
    imshow(uint8(I));
    title('org file');
    
    I_noise = I + 10 * randn(size(I));
    figure();
    imshow(uint8(I_noise));
    title('noise file');
    [m,n] = size(I_noise);
    
    DenoisedImg = zeros(m,n);
    PaddedImg = padarray(I,[1, 1],'symmetric','both');
    
    tic
    for i = 1: m
        for j = 1: n
            roi = PaddedImg(i:i+2, j:j+2);
            % first stage: average and median
            mean_V = mean(roi(:,2));
            mean_H = mean(roi(2,:));
            median_HV = median([mean_V, roi(2, 2)], mean_H);
            
            mean45 = mean([roi(1, 3), roi(2, 2), roi(3, 1)]);
            mean135 = mean([roi(1, 1), roi(2, 2), roi(3, 3)]);
            median_diag = median([mean45, roi(2, 2)], mean135);
            
            % second stage
            DenoisedImg(i, j) = median([median_HV, roi(2,2), median_diag]);
        end
    end
    figure();
    imshow(uint8(DenoisedImg));
    title('denoise file');
    toc
    
    b = medfilt2(I_noise,[3,3]);
    figure();
    imshow(uint8(b));
    title('median filter of matlab denoise file');
    
    

    多级中值有理混合滤波

    算法原理

    在这里插入图片描述

    WMF加权中值滤波

    在这里插入图片描述

    加权中值滤波也很好理解,和加权均值滤波差不多,就是在原始数据的基础上给每个点分别赋予一个权重,然后在加权后的数据中取出中位值作为滤波后的值。

    算法流程

    1. 求出’+'形和’X’形的窗口的中位值;

    2. 对’+'形窗口再利用CWMF求出一个值,CWMF是WMF的一种特殊情况,就是只对中心点进行加权;

    3. 对以上求出的三个参数用一下公式计算出一个新的值作为滤波后的值
      y ( m , n ) = ϕ 2 ( m , n ) + ϕ 1 ( m , n ) − 2 ∗ ϕ 2 ( m , n ) + ϕ 3 ( m , n ) h + k ( ϕ 1 ( m , n ) − ϕ 3 ( m , n ) ) y(m, n)=\phi 2(m, n)+\frac{\phi 1(m, n)-2 * \phi 2(m, n)+\phi 3(m, n)}{h+k(\phi 1(m, n)-\phi 3(m, n))} y(m,n)=ϕ2(m,n)+h+k(ϕ1(m,n)ϕ3(m,n))ϕ1(m,n)2ϕ2(m,n)+ϕ3(m,n)

    一样的,稍作改动该算法就可以用于raw格式图像

    在这里插入图片描述

    代码实现

    %% --------------------------------
    %% author:wtzhu
    %% email:wtzhu_13@163.com
    %% date: 20211023
    %% fuction: median rational hybird file filter
    %% references: MEDIAN-RATIONAL HYBRID FILTERS
    %% --------------------------------
    
    close all;
    clear all;
    clc
    img = imread('./images/lena.bmp');
    I = double(img);
    figure();
    imshow(uint8(I));
    title('org file');
    
    I_noise = I + 10 * randn(size(I));
    figure();
    imshow(uint8(I_noise));
    title('noise file');
    [m,n] = size(I_noise);
    
    DenoisedImg = zeros(m,n);
    PaddedImg = padarray(I,[1, 1],'symmetric','both');
    
    h = 2;
    k = 0.01;
    
    tic
    for i = 1: m
        for j = 1: n
            roi = PaddedImg(i:i+2, j:j+2);
            median_HV = median([roi(1,2), roi(2,1), roi(2,2), roi(2,3), roi(3,2)]);
            median_diag = median([roi(1,1), roi(1,3), roi(2,2), roi(3,1), roi(3,3)]);
            CWMF = median([roi(1,2), roi(2,1), roi(2,2)*3, roi(2,3), roi(3,2)]);
            
            DenoisedImg(i, j) = CWMF + (median_HV + median_diag - 2 * CWMF) / (h + k * (median_HV - median_diag));
        end
    end
    toc
    figure();
    imshow(uint8(DenoisedImg));
    title('denoise file');
    b = medfilt2(I_noise,[3,3]);
    figure();
    imshow(uint8(b));
    title('median filter of matlab denoise file');
    

    相关链接

    展开全文
  • 点云滤波——半径滤波

    千次阅读 2018-08-31 22:29:30
    如图1所示,有助于形象化理解RadiusOutlierRemoval的作用,在点云数据中,用户指定每个的一定范围内周围至少要有足够多的近邻。例如,如果指定至少要有1个邻居,只有黄色的会被删除,如果指定至少要有2个邻居,...

    如图1所示,有助于形象化理解RadiusOutlierRemoval的作用,在点云数据中,用户指定每个的点一定范围内周围至少要有足够多的近邻。例如,如果指定至少要有1个邻居,只有黄色的点会被删除,如果指定至少要有2个邻居,黄色和绿色的点都将被删除。
    半径滤波示意图
    图1 RadiusOutlierRemoval滤波处理示意图

    pcl::RadiusOutlierRemoval<pcl::PointXYZ> outrem;
    outrem.setInputCloud(pointCloud_raw);
    outrem.setRadiusSearch(0.02);//设置在0.02半径的范围内找邻近点
    outrem.setMinNeighborsInRadius(5);//设置查询点的邻近点集数小于5的删除
    outrem.filter(*pointCloud_filter);
    
    展开全文
  • stm32简易示器(标准库)

    万次阅读 多人点赞 2020-02-16 11:59:28
    此项案例是基于正点原子精英板制作的一个简易示器,可以读取信号的频率和幅值,并可以通过按键改变采样频率和控制屏幕的更新暂停。也是通过学习网上大佬的经验加以改编。下面将介绍所用到的基本原理以及相关代码的...
  • 图像卷积与滤波的一些知识

    万次阅读 多人点赞 2015-10-12 21:24:06
    图像卷积与滤波的一些知识zouxy09@qq.comhttp://blog.csdn.net/zouxy09 之前在学习CNN的时候,有对卷积经常一些学习和整理,后来就烂尾了,现在稍微整理下,先放上来,以提醒和交流。一、线性滤波与卷积的基本...
  • matlab与示器连接及电脑连接

    千次阅读 2021-04-18 11:33:14
    标签:最近进行了示器的数据采集,MSO2014,openChoice软件+Tekvisa驱动就可以了,采集的波形可以直接用matlab处理。后面又发现可以直接将示器跟matlab进行连接。...2.左边,test&measurement窗...
  • 分析——1. 初识小分析

    千次阅读 多人点赞 2018-07-17 12:07:25
    事实,自从小分析被提出后,已经在很多地方广泛使用,并且在一定程度取代了傅立叶分析。在这一章节里,我将为大家介绍这种神奇的小分析。 小变化/小分析 (Wavelet Transform) 首先,有请小函数: ...
  • matlab小阈值降噪

    千次阅读 2022-03-15 21:42:06
    一个合适的阈值,对噪声引起的小系数进行萎缩来去除噪声。 计算含噪信号的正交小波变换,选择分解层次N,进行N层小分解。 对小系数进行非线性阈值处理,保留低频系数,对1到N层的高频系数采用硬阈值和软...
  • PCL 几种滤波方法

    万次阅读 多人点赞 2018-07-19 10:51:09
    在获取点云数据时,由于设备精度、操作者经验、环境因素等带来的影响,以及电磁衍射特性、被测物体表面性质变化和数据拼接配准操作过程的影响,点云数据中将不可避免地出现一些噪声,属于随机误差。除此之外,...
  • 相位累加器输出的数据,作为波形存储器的相位采样地址,这样就可以把存储在波形存储器里的波形采样值经查表出,完成相位到幅度的转换。波形存储器的输出数据送到 D/A 转换器,由 D/A 转换器将数字信号转换成模拟...
  • 可是这样的题目并不常见,一次见到是什么时候我已经不记得了,昨天,一道让人忍不住叫好的作业题目摆在我的面前。 先看题目,是某大学大四学生的课程作业: ------- 题目 Heart Beat Period De...
  • 利用小分解后,频率计算问题

    千次阅读 2019-03-29 10:10:02
    matlab中使用小工具箱对信号进行小分解后,得到各频率分量的重构信号,分解后的这些信号的频段具体怎么计算??? 答: 小波变换并不是纯频域的变换,它无法完全脱离时空域,所以小的应用的多数领域并不十分...
  • 电信号实际是以电磁的形式在传输线中传播的。当传输线的尺寸不再远小于电磁波长时,就不得不考虑这个“”的特性了。下图是将一个窄脉冲施加到100m左右的终端短路的网线时,示器在信号源端测量到的图片。...
  • 峰峰值定义_示器峰峰值怎么看

    千次阅读 2021-04-25 14:16:05
    描述峰峰值定义峰峰值是指一个周期内信号最高值和最低值之间差的值,就是最大和最小之间的范围。...示器峰峰值怎么看一、对于模拟示器:用示器测量信号的电压峰-峰值和周期时,把波形峰-峰值高度调到显示屏...
  • 将市场销售的 6 V 直流电源接到收音机,为什么有的声音清晰,有的含有交流声? 对于同样标称输出电压为 6 V 的直流电源,在未接收音机时,为什么测量输出端子的电压,有的为 6 V,而有的为 7~8 V?用后者为收音机...
  • ARS408-21毫米雷达笔记

    千次阅读 多人点赞 2021-06-11 17:36:00
    ARS 408 , 德国大陆的毫米雷达: 型号: ARS408-21(XX, Far Range 250m), ARS408-21SC(特殊的软件版本, 扩展到1200m, 优化目标分类行人检测, 增加多边形配置等), 本篇测的是前者ARS408-21XX 发射频率76…77GHz, ...
  • 基于STM32的示

    万次阅读 多人点赞 2019-07-20 20:14:13
    网上的轮子,很好用。最大频率可以接近90KHZ,再往可能受限于中断的处理速度,因为ADC采样的频率一般可以设为12M,最大可达14M. 绘制波形时为前后连线,效果很好。 19.11.24更新 代码链接: 链接:...
  • 点云处理之点云滤波去噪

    千次阅读 2021-09-09 22:02:56
    目录 点云滤波简介 什么是点云滤波? 为什么要点云滤波?...点云滤波作为常见的点云处理算法,一般是点云处理的...有很多方面也有很多种功能,比如去除噪声、离群、点云平滑以及空洞、数据压缩等 为什么要点云滤
  • 利用STM32F103单片机输出SPWM

    万次阅读 多人点赞 2019-11-29 14:43:31
    最近需要用到单片机输出SPWM功能,在网上了好多资料,发现都不完整,有算法的没有代码,有代码的看不懂算法。于是只好自己摸索,现将方法整理如下。 关于什么是SPWM,为什么要用SPWM,网上的介绍有很多,就...
  • 器参数详解

    千次阅读 2020-05-19 08:32:09
    一、带宽 带宽:示器测得的单频信号幅度与实际幅度值相差小于3dB的频率范围。...经常看到示类似于这样的标注:500MHz,4GSa/s。采样率跟带宽是什么关系呢?经典的采样定理表达了这样一种意思:
  • 讲一个故事: 今年九月,一个新学期的开始,课很少。...一家柜台摆着“低价处理 LCD 模块”的牌子,对于像我这样的穷学生来说,价格往往是考虑的主要因素。我径直走了过去,老板说这些低价屏都...
  • 分析与神经网络 故障诊断

    千次阅读 2021-05-28 18:23:28
    分析能够将任何信号分解到一个由小伸缩而成的基函数族,对信号进行高、低频部分局部细化并保留原信号的时域特征,因而具有良好的时频特性,能对非平稳信号进行有效识别,达到故障诊断的目的。 小波变换与...
  • MATLAB--数字图像处理 中值滤波

    万次阅读 多人点赞 2019-09-14 16:23:46
    中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术,中值滤波的基本原理是把数字图像或数字序列中一点的值用该的一个邻域中各值的中值代替,让周围的像素值接近的真实值,从而消除孤立的噪声...
  • 均值滤波vs中值滤波

    千次阅读 2017-04-08 10:13:16
    中值滤波也是一种很常用的数字滤波器,它通过窗内的所有像素值的中值然后赋给中心像素,然后得到输出图像,这样做的好处是,它不创造新的像素值,只是取周围像素值作为它的输出,这一方法可以有效的消除脉冲噪声...
  • 关于点云滤波去噪的方法

    万次阅读 多人点赞 2019-01-24 23:04:04
    (2) 因为遮挡等问题造成离群需要去除 (3) 大量数据需要下采样 (4) 噪声数据需要去除 点云数据去噪滤波方法: 双边滤波、高斯滤波、分箱去噪、KD-Tree、直通滤波、随机采样一致性滤波等 方法定义以及...
  • (带宽定义:示器带宽的定义没有变,就是输入一个正弦,保持幅度不变,增加信号频率,当示显示的信号是实际信号幅度的70.7%(即3dB衰减)的时候,该对应的频率就等于示器带宽。)使用正弦信号发生器,在...
  • 利用MATLAB命令窗口绘制Simulink仿真示器波形的方法最近写了一篇有关步进电机控制仿真分析的文章,需要将一部分仿真波形图贴到WORD里面去。但贴图时发现,如果直接将simulink中示器的输出波形截图后贴到word文...
  • 锯齿FMCW雷达目标检测原理

    万次阅读 多人点赞 2019-06-18 23:03:45
    FMCW(Frequency Modulated Continuous Wave),即调频的连续信号。目前在安防、自动驾驶等领域的应用极为广泛。本文从最基本的数学原理出发进行推导,对锯齿体制FMCW雷达目标检测进行详细说明。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 88,390
精华内容 35,356
关键字:

上波找点