基于小波变换的信号降噪处理及仿真研究_信号处理方法推荐--1(转载自用,侵删)...

综述

作者:aresmiki
链接:https://www.zhihu.com/question/23701194/answer/167005497
来源:知乎
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。
非平稳信号处理应该是现在信号处理技术最新的也是最热的研究方向了,信号处理方法从最早的时域统计到傅里叶变换的频域分析,是人们认识信号本质的一次巨大飞跃,给了分析人员换个角度看世界的方法,这个时期傅里叶变换在调和分析,谐波分析等领域得到的巨大发展,当然工程人员将傅里叶变换也迅速运用的工程运用中,得到了巨大发展,特别是在机械故障诊断和 功能性磁共振成像fMRI 中。但是现实却不是理论那么美好,人们发现,傅里叶变换作为一个全局变换,天然的少了另一个维度,如果将时间域信号比作一个平面中的物体的话,那么频域信号也同样是一个平面中的物体,只是给我们换了一个角度而已,而人们总是对三维世界的物体更具有直观了解,平面的东西始终是不具体不形象,信号一样,工程人员总是想知道信号有哪些频率,且这些频率在何时产生,而这个需求就给分析方法提出了一个要求,必须多一个维度,也就是频域变换需要保留时间信息,聪明的Gabor第一次提出了窗口傅里叶变换,这个信号分析带来了时间尺度,也第一次有了时频分析的概念,当然其理论就不介绍了,而其优点当然就是同时给了我们时间和频率的信息。科学的魅力就是,永远没有完美的理论,窗口福利叶变换的窗口如何选择成为了一个难题,选的太大,时间分辨率太大,选的太小,频率分辨率太大,如何能在该大的时候大,该小的时候小呢?伟大的Morlet告诉我们,傅里叶变换本身就是全局的,而我们为什么还非得坚守傅里叶变换这个全局框架下呢?这种变化快的信号本身就是短时局部的信号不断变换,为什么不是一个个局部的去分析信号呢,这也是小波变换的基本思想,小波将时频分析推向了研究高潮,这一时期小波理论不断发展,出现了许多小波,db系列小波,morlet小波,coif系列小波等等,可以看出来,小波研究都是基于小波基函数的研究,它不像傅里叶变换一样,基是固定的,而小波基函数有很多形式。这种灵活性给了小波广泛的运用优势,但是,科学就是这样,一种方法的诞生必然伴随着缺陷的诞生,小波基的这种灵活性却给分析人员带来困扰,如何针对不同信号选择这些基呢?是否有一个通用基来处理信号分解任务呢?还有一个关键问题是,小波变换受到测不准原理的限制,不能无限制的细分时间和频率。当然研究人员根据实际处理任务提出了相对通用的小波基,如图像处理中的曲波变换和脊波变换,这一时期小波研究含有另一个方向,不是一个小波不能完全符合处理任务吗?不是不好选择小波基吗?研究人员提出把多个小波基一起用,必然可以得到更好的结果,这典型的代表是多小波的发展。当然普通小波变换存在时移缺陷和频率混叠的缺陷,而双树复数小波DTCWT得提出基本解决了该问题。也有针对小波基难以构造的缺陷,提出了提升框架的小波变换,但是小波变换就是小波变换,本质还是没有改变,小波具有的问题,后面发展出来的方法依然存在这些问题,小波处理方法没有过时,但是其诸多问题却在工程运用中越来越明显,工程运用中信号千变万化,小波基选择始终是摆在分析人员面前的问题,还有就是楼主提到的消噪,小波消噪软阈值和硬阈值的选择会对信号消噪能力有巨大影响,选择不好要么消除不够,有么太多。更直白的说,你说你消除了噪声,那么万一别人问你,你凭什么说你消除的是噪声而不是有用信号了,这个问题是一个无法回答的问题。当然后小波时代崛起了一批数据驱动的非线性非平稳信号处理方法,这些方法的提出就是摆脱小波基选择困难,测不准原理限制的问题,这类方法典型代表就是EMD,EMD方法根据信号自身特点不断分离开来,噪声和有用信号安照不同频带划分分离开,当然由于EMD方法本身诸多问题,后来发展出来了EEMD,CEEMD,CEEMDAN。CWEEMDAN等方法,但是其数学基础问题依然存疑,可是,EMD的提出给了人们一个启发,信号分解就是要让信号根据自己特点分开而不是人工的参与,一切应该是自驱动自适应才是最好的,受该思想的启发,小波研究者也希望将小波推广到这类方法中,其典型代表SWT和EWT,虽然方法很妙,但是其本质问题却并非数据驱动,有兴趣可以自行查阅文献,当然继承该思想的还有VMD、NSP、ITD方法,不得不提一下,这些方法都是非常适合非平稳非线性信号处理,但是分解中的一些参数设置却十分困难。关于信号分析最前沿的方法和理论,推荐专栏 信号处理与机器学习 - 知乎专栏 里面讲了很多信号处理方法的核心思想和最新方法。。。。

VMD(变分模态分解)(https://zhuanlan.zhihu.com/p/66898788?utm_source=qq)

这一篇写一下变分模态分解(原始论文:Variational Mode Decomposition),跟原始论文思路思路一致但有一点点不太一样,原始论文写的很好,但我不是通信专业没有学过信号相关课程一开始看起来有点费劲。模态分解认为信号是由不同“模态”的子信号叠加而成的,而变分模态分解则认为信号是由不同频率占优的子信号叠加而成的,其目的是要把信号分解成不同频率的子信号。变分模态分解的分解结果如图所示

79a1390353b59c3f9c12683df1a539aa.png

基础

一开始论文看不懂的原因是缺少相关前置知识,但一旦顺下来就会感觉其实没有那么难,难的是作者的思路很巧妙,先写下我遇到的这些知识盲点

第一点是傅立叶变换的微分性质

equation?tex=f%28t%29 的傅立叶变换为
equation?tex=F%28%5Comega%29 ,其导数
equation?tex=f%27%28t%29 的傅立叶变换为
equation?tex=j%5Comega+F%28%5Comega%29
equation?tex=%5Cbegin%7Balign%7D+f%28t%29+%26%3D+%5Cfrac%7B1%7D%7B%5Csqrt%7B2%5Cpi%7D%7D+%5Cint_%5Cmathbb%7BR%7DF%28%5Comega%29e%5E%7Bj%5Comega+t%7D%5Cmathrm%7Bd%7D%5Comega+%5C%5C+%5Cfrac%7B%5Cpartial+f%28t%29%7D%7B%5Cpartial+t%7D+%26%3D+%5Cfrac%7B1%7D%7B2%5Cpi%7D%5Cint_%5Cmathbb%7BR%7D%5Cfrac%7B%5Cpartial+F%28%5Comega%29e%5E%7Bj%5Comega+t%7D%7D%7B%5Cpartial+t%7D+%5Cmathrm+d%5Comega+%5C%5C+++++%26+%3D+%5Cfrac%7B1%7D%7B2%5Cpi%7D+%5Cint_%5Cmathbb%7BR%7D+j%5Comega+F%28%5Comega%29e%5E%7Bj%5Comega+t%7D%5Cmathrm+d+%5Comega+%5Cend%7Balign%7D+%5C%5C

另一点是解析信号,现实世界只能采集实信号,但实信号有很多不好用的性质,如存在负频率,无法直接得到调制频率后的实信号等。 设原始信号是一个实信号

equation?tex=f%28t%29+%3D+A%28t%29%5Ccos%28%5Cphi+t%29+%5C%5C为了方便表示
equation?tex=%5Cphi 为随
equation?tex=x 变化的函数,相当于瞬时频率,解析信号是一个复信号,可以通过希尔伯特变换得到
equation?tex=f_A%28t%29+%3D+f%28t%29+%2B+j%5Cmathcal%7BH%7Df%28t%29+%3D+A%28t%29e%5E%7Bj%5Cphi+t%7D+%5C%5C 解析信号的实部是原本的实数信号,并且经过调频之后复数信号的实部仍然是调频之后的实信号,如对信号
equation?tex=F_A%28t%29 增加频率
equation?tex=%5Comega_0+ ,只需要乘以
equation?tex=e%5E%7Bj%5Comega_0+t%7D 即可
equation?tex=%5Cbegin%7Balign%7D+f_A%28t%29%5Ccdot+e%5E%7Bj%5Comega_0+t%7D+%26+%3D+A%28t%29e%5E%7Bj%5Cphi+t%7D+%5Ccdot+e%7Bj%5Comega_0+t%7D+%5C%5C+++++%26+%3D+A%28t%29%5Ccdot+e%5E%7Bj%28%5Cphi+%2B+%5Comega_0%29t%7D+%5C%5C+%5CRe+%5CLarge%5C%7B+F_A%28t%29e%5E%7Bj%5Comega_0+t%7D+%5CLarge%5C%7D+%26+%3D+A%28t%29%5Ccdot+%5Ccos%28%28%5Cphi%2B%5Comega_0%29t%29+%5Cend%7Balign%7D+%5C%5C 由此,其实部相当于在原本频率
equation?tex=%5Cphi 的基础上增加了频率
equation?tex=%5Comega_0 ,如下面的
matlab脚本
clear;close all;clc;
t = 1:0.01:10;
%%
f1 = sin(20*t).*(t-5).^2;
subplot(3,1,1);
plot(f1);
ylim([-25 25]);
%%
f2 = sin(50*t).*(t-5).^2;
subplot(3,1,2);
plot(f2);
ylim([-25 25]);
%%
H = hilbert(f1);
f_hat = H.*exp(1i*30.*t);
subplot(3,1,3);
plot(real(f_hat));
ylim([-25 25]);

f713e1f76e51fef95e1c378acad9e171.png

matlabhilbert函数包括希尔伯特变换和解析函数转换两部分,直接得到实信号的解析信号,其中希尔伯特变换

equation?tex=%5Cmathcal%7BH%7D+f+%3D+%5Cfrac%7B1%7D%7B%5Cpi%7D+%2A+f++%5C%5C

正文

接下来我们看看如何一步一步得到变分模态分解的思路

原论文通过一个信号降噪问题进行说明,现需要对采样信号

equation?tex=f_0 进行降噪重构,假设观测信号是由原始信号叠加一个独立的高斯噪音
equation?tex=f_0+%3D+f+%2B+%5Ceta+%5C%5C,需要求
equation?tex=f ,又说该等式是一个不适定问题(ill-posed problem),不满足识定问题的三个条件,所以要用一个正则化的方法
equation?tex=%5Cunderset%7Bf%7D%7B%5Cmin%7D%7B%5CLarge%7B%7D+%5Cleft+%5C%7C+f-f_0+%5Cright+%5C%7C_2%5E2+%2B+%5Calpha+%5Cleft+%5C%7C+%5Cpartial_t+f+%5Cright+%5C%7C_2%5E2+%7B%5CLarge%7D%7D%5C%5C 第一部分是对原信号进行重构,第二部分是为了解决不适定问题的解不唯一,而且不同于机器学习的建模问题一样
equation?tex=f 是一个权值形式可以直接加权值
equation?tex=W 的L1-norm或L2-norm,这里的
equation?tex=f 是一个纯函数的形式,其导数的L2-norm最小化感觉上应该是保证了函数
equation?tex=f 不会产生太大的波动,这里可不可以跟防止过拟合联系起来,接下来说明这个约束会使
equation?tex=f 在频率域产生什么影响。

下面解这个式子,首先来看一下直接最小化泛函

equation?tex=J 能不能解出来
equation?tex=%5Cbegin%7Balign%7D++J%5Bf%5D+%26%3D+%5Cint_%5Cmathbb%7BR%7D%28f%28t%29+-+f_0%28t%29%29%5E2+%2B+%5Calpha+%28f%27%28t%29%29%5E2+%5Cmathrm%7Bd%7Dt+%5C%5C+++++%26+%3D+%5Cint_%5Cmathbb%7BR%7DF%28f%2C+f%27%29%5Cmathrm%7Bd%7Dt+%5Cend%7Balign%7D+%5C%5C
equation?tex=F%28f%2Cf%27%29+%3D+%28f-f_0%29%5E2+%2B+%5Calpha%28f%27%29%5E2 ,然后根据E-L方程的引理
equation?tex=%5Cfrac%7B%5Cpartial+F%7D%7B%5Cpartial+f%27%7D+f%27+-+F+%3D+C%5C%5C
equation?tex=%5Calpha%7Bf%27%7D%5E2+-+%28f-f_0%29%5E2+%3D+C%5C%5C ,往下忘了怎么求了。。。。。想起来怎么求再继续写-_-!,这个直接求的方法与原文没什么关系的

由上面的傅立叶变换的微分性质知道,用傅立叶变换在复数域很方便可以把微分约掉,而且还有一个叫什么Plancherel傅立叶等距映射的东西,时域的L2范数与傅立叶变换到的频率域的L2范数等距,故上面的最小化的项可以直接用傅立叶变换转换到频域

equation?tex=%5Cunderset%7B%5Chat%7Bf%7D%28%5Comega%29%7D%7B%5Cmin%7D%7B%5CLarge%7B%7D+%5Cleft+%5C%7C+%5Chat%7Bf%7D%28%5Comega%29-%5Chat%7Bf_0%7D%28%5Comega%29+%5Cright+%5C%7C2%5E2+%2B+%5Calpha+%5Cleft+%5C%7C+j%5Comega%5Chat%7Bf%7D%28%5Comega%29+%5Cright+%5C%7C_2%5E2+%7B%5CLarge%7D%7D+%5C%5C
,这里需要注意的是利用上面的那个什么等距定理有
equation?tex=%5Chat%7Bf%7D%28%5Comega%29
equation?tex=%5Chat%7Bf_0%7D%28%5Comega%29
equation?tex=%5Comega 是同一个
equation?tex=%5Comega ,这一点很重要,第二个是L2范数大于0,则把其展开为泛函
equation?tex=J%5B%5Chat%7Bf%7D%5D+%3D+%5Cint%5Cmathbb%7BR%7D%28%5Chat%7Bf%7D%28%5Comega%29+-+%5Chat%7Bf_0%7D%28%5Comega%29%29%5E2+%2B+%5Calpha%5Comega%5E2%28%5Chat%7Bf%7D%28%5Comega%29%29%5E2+%5Cmathrm%7Bd%7D%5Comega%5C%5C 求最极值,由于这里面只有
equation?tex=%5Chat%7Bf%7D 直接求偏导即可也不用解微分方程
equation?tex=%5Cbegin%7Balign%7D+%5Cfrac%7B%5Cdelta+J%5B%5Chat%7Bf%7D%5D%7D%7B%5Cdelta+%5Chat%7Bf%7D%7D+%26%3D+2%28%5Chat%7Bf%7D+-+%5Chat%7Bf_0%7D%29+%2B+2%5Calpha%5Comega%5E2%5Chat%7Bf%7D+%5C%5C+++++%26%3D+2%28%5Chat%7Bf%7D%281%2B%5Calpha%5Comega%5E2%29+-+%5Chat%7Bf_0%7D%29+%3D+0+%5C%5C+%5Chat%7Bf%7D+%26%3D+%5Cfrac%7B%5Chat%7Bf_0%7D%7D%7B1+%2B+%5Calpha%5Comega%5E2%7D+%5Cend%7Balign%7D%5C%5C

可以看到得到的

equation?tex=%5Chat%7Bf%7D 相当于对观测信号在
equation?tex=%5Chat%7Bf_0%7D 在频率段进行滤波,过滤掉了高频部分,这说明加了该导数的L1正则约束与上面的直观感觉是一样的,过滤掉了高频部分,减弱
equation?tex=f 的波动

再往下进行,模态分解需要把原始信号分解成多个子信号的和,我们为了和原文对应修改一下符号表示,

equation?tex=f%28t%29 表示观测的采样信号,
equation?tex=u_k%28t%29 表示分解得到的基函数,则上面的约束对象变为
equation?tex=%5Cunderset%7Bu_k%7D%7B%5Cmin%7D+%5Cleft%5C%7B++%7C%7C+%5Csum_k+u_k%28t%29+-+f%28t%29%7C%7C2%5E2+%2B+%5Calpha%5Csum_k%7C%7C+%5Cpartial+t+%5Bu_k%28t%29%5D+%7C%7C2%5E2++%5Cright%5C%7D+%5C%5C 同样先转化为频率域再求极值
equation?tex=%5Cunderset%7B%5Chat%7Bu_k%7D%7D%7B%5Cmin%7D+%5Cleft%5C%7B++%7C%7C+%5Csum_k+%5Chat%7Bu_k%7D%28%5Comega%29+-+%5Chat%7Bf%7D%28%5Comega%29%7C%7C_2%5E2+%2B+%5Calpha%5Csum_k%7C%7C+j%5Comega+%5Chat%7Bu_k%7D%28%5Comega%29+%7C%7C_2%5E2++%5Cright%5C%7D+%5C%5C+ 泛函
equation?tex=J%5B%5Chat%7Bu_1%7D%2C+%5Chat%7Bu_2%7D%2C+%5Ccdots%2C%5Chat%7Bu_K%7D%5D+%3D+%5Cint%5Cmathbb%7BR%7D%28%5Csum_k%5Chat%7Bu_k%7D%28%5Comega%29+-+%5Chat%7Bf%7D%28%5Comega%29%29%5E2+%2B+%5Calpha%5Comega%5E2%5Csum_k%28%5Chat%7Bu_k%7D%28%5Comega%29%29%5E2+%5Cmathrm%7Bd%7D%5Comega+%5C%5C

equation?tex=%5Cbegin%7Balign%7D+%5Cfrac%7B%5Cdelta+J%7D%7B%5Cdelta+u_k%7D+%26%3D+2+%5Cleft%28%5Csum_k+%5Chat%7Bu_k%7D+-++%5Chat%7Bf%7D%5Cright%29+%2B+2+%5Calpha+%5Comega+%5Cleft%28%5Chat%7Bu_k%7D+%5Cright%29+%3D+0+%5C%5C+%5Chat%7Bu_k%7D%28%5Comega%29+%26%3D+%5Cfrac%7B%5Chat%7Bf%7D%28%5Comega%29+-+%5Csum_%7Bi%5Cneq+k%7D%5Chat%7Bu_i%7D%28%5Comega%29%7D%7B1+%2B+%5Calpha+%5Comega%5E2%7D+%5Cend%7Balign%7D+%5C%5C 每个基函数基于其他的基函数更新,相当于每个基函数是原信号剩余部分的低通滤波,每次迭代都是保留剩余信号的低频率部分。

到现在为止我们发现每个基函数都会趋向于每次的剩余信号分量的低频部分,这与我们原始的假设“每个基函数都有不同的频率分量”是相悖的,但根据上面的低通滤波的性质,每个基函数进行特定频率的滤波应该就能解决这个问题了,那么上面的式子就简单的变为

equation?tex=%5Chat%7Bu_k%7D%28%5Comega%29+%3D+%5Cfrac%7B%5Chat%7Bf%7D%28%5Comega%29+-+%5Csum_%7Bi%5Cneq+k%7D%5Chat%7Bu_i%7D%28%5Comega%29%7D%7B1+%2B+%5Calpha+%28%5Comega+-+%5Comega_k%29%5E2%7D+%5C%5C 其中,
equation?tex=%5Comega_k 为每个基函数
equation?tex=u_k 的中心频率,该式就是变分模态分解的基函数的更新公式,我们来看一下这个式子应该如何得到,以便于找到中心频率的更新公式

由上面的一步一步的演化发现,基函数对剩余信号的低通滤波是由导数的L2正则最小化带来的,要得到基函数的中心频率约束也要从这个地方入手,由上面的推导可知每个基函数都会被约束到

equation?tex=0 频率附近,那么我们把基函数的频率增加各自的中心频率
equation?tex=%5Comega_k 得到
equation?tex=%5Chat%7Bu_k%7D%28%5Comega%2B%5Comega_k%29 ,并保证
equation?tex=%5Comega+%5Cgeq+0 即可,则相当于对每个基函数乘以了一个
equation?tex=e%5E%7B-j%5Comega_k+t%7D ,这里需要对频率进行变换,我们沿用文章开头的解析信号的性质,认为
equation?tex=u_k%28t%29 是一个复数的解析信号,同样也需要把观测信号
equation?tex=f 预先转化为解析信号,下文默认都是解析信号,则每一个基函数转换频率后的导数变为

equation?tex=%5Cpartial_t+u_k%28t%29+e%5E%7B-j%5Comega_k+t%7D%5C%5C傅立叶变换得
equation?tex=j%5Comega%5Chat%7Bu_k%7D%28%5Comega%2B%5Comega_k%29%5C%5C 先来看看傅立叶变换为什么会得到这个,而不是
equation?tex=j%28%5Comega+%2B+%5Comega_k%29%5Chat%7Bu_k%7D%28%5Comega%2B%5Comega_k%29 反傅立叶变换
equation?tex=%5Cbegin%7Balign%7D+++++u_k%28t%29+%26%3D+%5Cint_0%5E%5Cinfty%5Chat%7Bu_k%7D%28%5Comega%2B%5Comega_k%29e%5E%7Bj%5Comega+t%7D%5Cmathrm%7Bd%7D%5Comega+%5C%5C+++++%5Cpartial_t+u_k%28t%29+%26+%3D+%5Cint_0%5E%5Cinfty+j%5Comega+%5Chat%7Bu_k%7D%28%5Comega%2B%5Comega_k%29e%5E%7Bj%5Comega+t%7D%5Cmathrm%7Bd%7D%5Comega+%5Cend%7Balign%7D+%5C%5C 感觉应该不是这样证的,我数学不是很好,不怎么会这个变量变换

至此,约束变为

equation?tex=%5Cunderset%7Bu_k%7D%7B%5Cmin%7D+%5Cleft%5C%7B++%7C%7C+%5Csum_k+u_k%28t%29+-+f%28t%29%7C%7C2%5E2+%2B+%5Calpha%5Csum_k%5Cleft%5C%7C+%5Cpartial+t+%5Cleft%5Bu_k%28t%29e%5E%7B-j%5Comega_k+t%7D%5Cright%5D+%5Cright%5C%7C_2%5E2++%5Cright%5C%7D+%5C%5C 傅立叶变换
equation?tex=%5Cunderset%7B%5Chat%7Bu_k%7D%7D%7B%5Cmin%7D+%5Cleft%5C%7B++%7C%7C+%5Csum_k+%5Chat%7Bu_k%7D%28%5Comega%29+-+%5Chat%7Bf%7D%28%5Comega%29%7C%7C_2%5E2+%2B+%5Calpha%5Csum_k%5C%7C+j%5Comega+%5Chat%7Bu_k%7D%28%5Comega+%2B+%5Comega_k%29+%5C%7C_2%5E2++%5Cright%5C%7D%5C%5C 根据论文原文把
equation?tex=%5Comega%5Cleftarrow+%5Comega+-+%5Comega_k
equation?tex=%5Cunderset%7B%5Chat%7Bu_k%7D%7D%7B%5Cmin%7D+%5Cleft%5C%7B++%7C%7C+%5Csum_k+%5Chat%7Bu_k%7D%28%5Comega%29+-+%5Chat%7Bf%7D%28%5Comega%29%7C%7C_2%5E2+%2B+%5Calpha%5Csum_k%5C%7C+j%28%5Comega+-+%5Comega_k%29+%5Chat%7Bu_k%7D%28%5Comega%29+%5C%7C_2%5E2++%5Cright%5C%7D+%5C%5C 然后求最小就可以得到上面
equation?tex=%5Chat%7Bu_k%7D 的更新公式
equation?tex=J+%3D+%5Cint_0%5E%5Cinfty%28%5Csum_k%5Chat%7Bu_k%7D%28%5Comega%29+-+%5Chat%7Bf%7D%28%5Comega%29%29%5E2+%2B+%5Calpha%5Csum_k%28%5Comega+-+%5Comega_k%29%5E2%28%5Chat%7Bu_k%7D%28%5Comega%29%29%5E2+%5Cmathrm%7Bd%7D%5Comega+%5C%5C

equation?tex=%5Cbegin%7Balign%7D+%5Cfrac%7B%5Cdelta+J%7D%7B%5Cdelta+%5Comega_k%7D+%26%3D+%5Cint_0%5E%5Cinfty+2%28%5Comega+-+%5Comega_k%29%7C%5Chat%7Bu_k%7D%28%5Comega%29%7C%5E2%5Cmathrm%7Bd%7D%5Comega++%5C%5C+%5Comega_k+%26%3D+%5Cfrac%7B%5Cint_0%5E%5Cinfty+%5Comega+%7C%5Chat%7Bu_k%7D%28%5Comega%29%7C%5E2%5Cmathrm%7Bd%7D%5Comega%7D%7B%5Cint_0%5E%5Cinfty+%7C%5Chat%7Bu_k%7D%28%5Comega%29%7C%5E2%5Cmathrm%7Bd%7D%5Comega%7D+%5Cend%7Balign%7D%5C%5C

最后,为了保证每个点处的重构信号与原信号尽可能相似,增加了每个点处的重构约束,其实这一项并不是必需的,最终的约束对象为

equation?tex=%5Cbegin%7Bmultline%7D+%5Cunderset%7B%7Bu_k%7D%2C%7B%5Comega_k%7D%7D%7B%5Cmin%7D+%5Cleft%5C%7B+%7C%7C+%5Csum_k+u_k%28t%29+-+f%28t%29%7C%7C2%5E2+%2B+%5Calpha%5Csum_k%5Cleft%5C%7C+%5Cpartial+t+%5Cleft%5Bu_k%28t%29e%5E%7B-j%5Comega_k+t%7D%5Cright%5D+%5Cright%5C%7C_2%5E2++%5Cright+%5C%7D+%5C%5C+s.t.+%5Cquad+f+%3D+%5Csum_k+u_k+%5Cend%7Bmultline%7D+%5C%5C ,然后拉格朗日乘子法带进去,但其需要满足下式才有意义
equation?tex=+%5Cint_%5Cmathbb%7BR%7D+%5Clambda%28t%29%5Cleft%28+f%28t%29+-+%5Csum_k+u_k%28t%29+%5Cright%29%5Cmathrm%7Bd%7D+t+%3D+%5Cint_%5Cmathbb%7BR%7D+%5Chat%7B%5Clambda%7D%28%5Comega%29%5Cleft%28+%5Chat%7Bf%7D%28%5Comega%29+-+%5Csum_k+%5Chat%7Bu_k%7D%28%5Comega%29%5Cright%29%5Cmathrm%7Bd%7D%5Comega++%5C%5C 这个式子还是符合Parseval定理,故整理一下
equation?tex=+%5Cbegin%7Bmultline%7D+J+%3D+%5Cint_0%5E%5Cinfty+%5Calpha%5Csum_k%28%5Comega+-+%5Comega_k%29%5E2%28%5Chat%7Bu_k%7D%28%5Comega%29%29%5E2+%2B+%5Cleft%28%5Csum_k%5Chat%7Bu_k%7D%28%5Comega%29+-+%5Chat%7Bf%7D%28%5Comega%29%5Cright%29%5E2+%5C%5C+%2B+%5Chat%7B%5Clambda%7D%28%5Comega%29%5Cleft%28+%5Chat%7Bf%7D%28%5Comega%29+-+%5Csum_k+%5Chat%7Bu_k%7D%28%5Comega%29%5Cright%29+%5Cmathrm%7Bd%7D%5Comega+%5Cend%7Bmultline%7D+%5C%5C

最后整理一下更新公式:

equation?tex=%5Chat%7Bu_k%7D%28%5Comega%29+%3D+%5Cfrac%7B%5Chat%7Bf%7D%28%5Comega%29+-+%5Csum_%7Bi%5Cneq+k%7D%5Chat%7Bu_i%7D%28%5Comega%29+%2B+%5Cfrac%7B%5Chat%7B%5Clambda%7D%28%5Comega%29%7D%7B2%7D%7D%7B1+%2B+%5Calpha+%28%5Comega+-+%5Comega_k%29%5E2%7D+%5C%5C+%5Comega_k+%3D+%5Cfrac%7B%5Cint_0%5E%5Cinfty+%5Comega%7C%5Chat%7Bu_k%7D%28%5Comega%29%7C%5E2+%5Cmathrm%7Bd%7D%5Comega%7D%7B%5Cint_0%5E%5Cinfty+%7C%5Chat%7Bu_k%7D%28%5Comega%29%7C%5E2+%5Cmathrm%7Bd%7D%5Comega%7D+%5C%5C+%5Chat%7B%5Clambda%7D%5E%7Bn%2B1%7D%28%5Comega%29+%3D+%5Chat%7B%5Clambda%7D%5En+%28%5Comega%29+-+%5Ctau%5Cleft%28+%5Chat%7Bf%7D%28%5Comega%29+-+%5Csum_k+%5Chat%7Bu_k%7D%28%5Comega%29+%5Cright%29+%5C%5C 其中
equation?tex=%5Chat%7B%5Clambda%7D 使用梯度下降更新

代码

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbert
T = 1000
fs = 1./T
t = np.linspace(0, 1, 1000,endpoint=True)
f_1 = 10
f_2 = 50
f_3 = 100
mode_1 = (2 * t) ** 2
mode_2 = np.sin(2 * np.pi * f_1 * t)
mode_3 = np.sin(2 * np.pi * f_2 * t)
mode_4 = np.sin(2 * np.pi * f_3 * t)
f = mode_1 + mode_2 + mode_3 + mode_4 + 0.5 * np.random.randn(1000)plt.figure(figsize=(6,3), dpi=150)
plt.plot(f, linewidth=1)

968d7feab34fba6f100a4213d6a0d505.png
class VMD:def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):""":param K: 模态数:param alpha: 每个模态初始中心约束强度:param tau: 对偶项的梯度下降学习率:param tol: 终止阈值:param maxIters: 最大迭代次数:param eps: eps"""self.K =Kself.alpha = alphaself.tau = tauself.tol = tolself.maxIters = maxItersself.eps = epsdef __call__(self, f):T = f.shape[0]t = np.linspace(1, T, T) / Tomega = t - 1. / T# 转换为解析信号f = hilbert(f)f_hat = np.fft.fft(f)u_hat = np.zeros((self.K, T), dtype=np.complex)omega_K = np.zeros((self.K,))lambda_hat = np.zeros((T,), dtype=np.complex)# 用以判断u_hat_pre = np.zeros((self.K, T), dtype=np.complex)u_D = self.tol + self.eps# 迭代n = 0while n < self.maxIters and u_D > self.tol:for k in range(self.K):# u_hatsum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]res = f_hat - sum_u_hatu_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)# omegau_hat_k_2 = np.abs(u_hat[k, :]) ** 2omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)# lambda_hatsum_u_hat = np.sum(u_hat, axis=0)res = f_hat - sum_u_hatlambda_hat -= self.tau * resn += 1u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)u_hat_pre[::] = u_hat[::]# 重构,反傅立叶之后取实部u = np.real(np.fft.ifft(u_hat, axis=-1))omega_K = omega_K * Tidx = np.argsort(omega_K)omega_K = omega_K[idx]u = u[idx, :]return u, omega_K
K = 4
alpha = 2000
tau = 1e-6
vmd = VMD(K, alpha, tau)
u, omega_K = vmd(f)
omega_K
# array([0.85049797, 10.08516203, 50.0835613, 100.13259275]))
plt.figure(figsize=(5,7), dpi=200)
plt.subplot(4,1,1)
plt.plot(mode_1, linewidth=0.5, linestyle='--')
plt.plot(u[0,:], linewidth=0.2, c='r')plt.subplot(4,1,2)
plt.plot(mode_2, linewidth=0.5, linestyle='--')
plt.plot(u[1,:], linewidth=0.2, c='r')plt.subplot(4,1,3)
plt.plot(mode_3, linewidth=0.5, linestyle='--')
plt.plot(u[2,:], linewidth=0.2, c='r')plt.subplot(4,1,4)
plt.plot(mode_4, linewidth=0.5, linestyle='--')
plt.plot(u[3,:], linewidth=0.2, c='r')
# [<matplotlib.lines.Line2D at 0x7fac532f4dd8>]

083da026f92a27bc5ecff48f38cb3e33.png

可以看到有比较强的端点效应,端点处会有重叠,文章原始论文中采用的方法是对称拼接的方法,把原信号复制一份然后拼成两半,一半对称放前面,一般对称放后面

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbertT = 1000
fs = 1./T
t = np.linspace(0, 1, 1000,endpoint=True)f_1 = 10
f_2 = 50
f_3 = 100
mode_1 = np.sin(2 * np.pi * f_1 * t)
mode_2 = np.sin(2 * np.pi * f_2 * t)
mode_3 = np.sin(2 * np.pi * f_3 * t)
f = np.concatenate((mode_1[:301], mode_2[301:701], mode_3[701:])) + 0.1 * np.random.randn(1000)plt.figure(figsize=(6,3), dpi=150)
plt.plot(f, linewidth=0.5)
# [<matplotlib.lines.Line2D at 0x7fc2134b8630>]

5e8f700484bb00e5fec7ead75c6038b3.png
class VMD:def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):""":param K: 模态数:param alpha: 每个模态初始中心约束强度:param tau: 对偶项的梯度下降学习率:param tol: 终止阈值:param maxIters: 最大迭代次数:param eps: eps"""self.K =Kself.alpha = alphaself.tau = tauself.tol = tolself.maxIters = maxItersself.eps = epsdef __call__(self, f):N = f.shape[0]# 对称拼接f = np.concatenate((f[:N//2][::-1], f, f[N//2:][::-1]))T = f.shape[0]t = np.linspace(1, T, T) / Tomega = t - 1. / T# 转换为解析信号f = hilbert(f)f_hat = np.fft.fft(f)u_hat = np.zeros((self.K, T), dtype=np.complex)omega_K = np.zeros((self.K,))lambda_hat = np.zeros((T,), dtype=np.complex)# 用以判断u_hat_pre = np.zeros((self.K, T), dtype=np.complex)u_D = self.tol + self.eps# 迭代n = 0while n < self.maxIters and u_D > self.tol:for k in range(self.K):# u_hatsum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]res = f_hat - sum_u_hatu_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)# omegau_hat_k_2 = np.abs(u_hat[k, :]) ** 2omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)# lambda_hatsum_u_hat = np.sum(u_hat, axis=0)res = f_hat - sum_u_hatlambda_hat -= self.tau * resn += 1u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)u_hat_pre[::] = u_hat[::]# 重构,反傅立叶之后取实部u = np.real(np.fft.ifft(u_hat, axis=-1))u = u[:, N//2 : N//2 + N]omega_K = omega_K * T / 2idx = np.argsort(omega_K)omega_K = omega_K[idx]u = u[idx, :]return u, omega_K
K = 3
alpha = 2000
tau = 1e-6
vmd = VMD(K, alpha, tau)u, omega_K = vmd(f)
omega_K
# array([  9.68579292,  50.05232833, 100.12321047])plt.figure(figsize=(5,7), dpi=200)
plt.subplot(3,1,1)
plt.plot(mode_1, linewidth=0.5, linestyle='--')
plt.plot(u[0,:], linewidth=0.2, c='r')plt.subplot(3,1,2)
plt.plot(mode_2, linewidth=0.5, linestyle='--')
plt.plot(u[1,:], linewidth=0.2, c='r')plt.subplot(3,1,3)
plt.plot(mode_3, linewidth=0.5, linestyle='--')
plt.plot(u[2,:], linewidth=0.2, c='r')
# [<matplotlib.lines.Line2D at 0x7fc2134075c0>]

ce15fee789d11fc592e2a67177188ed1.png

好像结果要好一点,VMD的一个缺点是K的值对结果有很大影响,但这个迭代过程,怎么开始怎么迭代怎么结束都可以自己控制,感觉可以按照自己的需求来定制怎么动态决定K

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/290276.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

js温故而知新11(AJAX)——学习廖雪峰的js教程

Asynchronous JavaScript and XML&#xff0c;意思就是用JavaScript执行异步网络请求。 如果仔细观察一个Form的提交&#xff0c;你就会发现&#xff0c;一旦用户点击“Submit”按钮&#xff0c;表单开始提交&#xff0c;浏览器就会刷新页面&#xff0c;然后在新页面里告诉你操…

linux下重启mysql php nginx

# /etc/init.d/mysql restart # /etc/init.d/php-fpm restart # /etc/init.d/nginx restart转载于:https://www.cnblogs.com/yiluxiuxing/p/4140627.html

C和指针之字符串编程练习9(在参数1中查找匹配参数2额任意字符)

1、问题 函数应该在第一个参数中进行查找&#xff0c;并返回匹配第二个参数所包含的字符的数目 2、代码实现 #include <stdio.h> #include <string.h>//函数应该在第一个参数中进行查找&#xff0c;并返回匹配第二个参数所包含的字符的数目 int count_chars(char…

最流行的 .NET 开源项目合集

Github 上有很多优秀的 .NET 开源项目&#xff0c;它们很适合用来学习 .NET 、设计模式和架构。nopCommerce https://github.com/nopSolutions/nopCommercestar: 7k ⭐nopCommerce 是最受欢迎最好的开源电子商务购物车解决方案之一&#xff0c;它基于 ASP.NET Core&#xff…

绿米开关如何重置_三种墙壁开关:绿米、调光、凌动/米家之间的异同

三种墙壁开关&#xff1a;绿米、调光、凌动/米家之间的异同2020-06-18 18:17:179点赞28收藏6评论三种墙壁开关&#xff1a;绿米、调光、凌动/米家之间的异同详细讲解绿米Aqara智能墙壁开关、Yeelight调光开关(86盒版)、Yeelight凌动开关/米家墙壁开关(米家墙壁开关也是Yeelight…

NES模拟器开发-PPU笔记

20151008 占坑,暂时没弄清楚PPU数据如何初始化,绘制顺序等. 转载于:https://www.cnblogs.com/Red_angelX/p/4860497.html

phpRedisAdmin 安装

phpRedisAdmin是一个web端管理redis的工具&#xff0c;每次都是命令行操作&#xff0c;今天安装了下&#xff0c;做个笔记[rootlocalhost linshi]# git clone https://github.com/ErikDubbelboer/phpRedisAdmin.git Initialized empty Git repository in /root/linshi/phpRedis…

GFS(Genetic Fuzzy Systems)—— 基于专家先验规则库和遗传算法相结合的智能体决策算法

文章目录1. FIS 系统&#xff08;Fuzzy Inference Systems&#xff09;1.1 什么是 FIS 系统&#xff1f;1.2 使用 FIS 算法的几个步骤2. GFS 系统&#xff08;GA FIS&#xff09;2.1 什么是基因遗传算法&#xff08;GA&#xff09;?2.2 使用GA算法进化FIS规则库在大规模的多智…

3-unit1 IPv6网络的管理

##########IPv6网络的管理#############学习目标了解IPv6管理IPv6##)IPv6简介Internet Protocol Version 6IPv6是IETF&#xff08;和互联网工程任务组&#xff09;设计的用与替代现行版本IP协议的下一代IP协议。IPv6采用128位2进制数码表示##IPv6示意图:##IPv6表示方式为方便操…

C和指针之字符串编程练习8实现char *my_strnchr(char const *str, int ch, int which)

1、问题 编写函数类似strchr函数,但是它有3个参数,第三个参数是1, 这个函数的功能就和strchr完全一样, 如果第三个参数是2,这个函数就返回一个指向ch字符在str字符串第二次出现的位置的指针,以此类推 2、代码实现 1 #include <stdio.h>2 #include <string.h>3…

Xamarin效果第一篇之时间轴

一直都想找个时间玩玩移动端,中间也去各种的调研快速的方式去开发;过程中还是太浮躁木有沉下心去学习;好早尝试过Flutter,二点让我不爽:1、配置环境好费劲(VS把我惯坏了)&#xff1b;2、套娃的写法真是恶心;来看看酸爽不:因为一直都是C#开发,所以这次再次摸索Xamarin去开发;前面…

bootstrap tabale 点击_jquery+bootstrap实现tab切换, 每次切换时都请求数据, 点击提交分别向不同的地址提交数据...

今天一个朋友叫帮做一个tab切换, 每一个tab内容区域都是从后台取出的数据, 这些数据要用表格的形式显示处理, 并且表格的内容区域可以修改, 如下所示:例子查看请演示查看.截图如图所示:实现步骤:几个需要注意的点:1. tab部分加一个data-id, 当中的id与下面要显示的具体内容的ta…

AJAX+json+jquery实现预加载瀑布流布局

AJAXjsonjquery实现预加载瀑布流布局转载于:https://www.cnblogs.com/zhujiabin/p/4860954.html

Lync 小技巧-42-动态-IP-统一沟通-环境-IP-变更后-操作

1. 查看-你的-公网IPhttp://www.ip138.com/2. 修改-你的-公网A记录https://www.godaddy.com/3. 修改-你的-拓朴-For-边缘服务器3.1.远程-前端服务器3.2.下载-拓朴3.3.选择-边缘服务器3.4.选择-边缘服务器3.5.修改-公网IP116.230.255.783.5.发布-拓朴3.6.导出-拓朴3.7.复制-拓朴…

C和指针之字符串memcpy、memmove、memset使用总结

1、介绍memcpy、memmove、memset 1) void *memcpy(void *dest, const void *src, size_t n); 从源src所指的内存地址的起始位置开始拷贝n个字节到目标dest所指的内存地址的起始位置中 2) void *memmove( void* dest, const void* src, size_t count ); 从src拷贝count个字节…

Blazor University (1)介绍 - 什么是 Blazor?

原文链接&#xff1a;https://blazor-university.com/overview/what-is-blazor/什么是 Blazor&#xff1f;Blazor 是一个单页应用程序开发框架。Blazor 这个名称是单词 Browser 和 Razor&#xff08;.NET HTML 视图生成引擎&#xff09;的组合/变异。这意味着 Blazor 不必在服务…

Git在版本2.13中继续改进了安全性和UI

Git的最新版本做了很多改进&#xff0c;旨在改进其用户界面&#xff0c;同时也修复了两个重要的漏洞。\\众所周知&#xff0c;Git用于唯一识别对象的SHA-1哈希算法最近被证明容易受到碰撞攻击。虽然Git团队准备过渡到一个新的更安全的散列算法&#xff0c;但它们已经实现了一种…

Mysql日期和时间函数

转载自:http://fanqiang.chinaunix.net/a2/b1/20010705/150000802.html 对于每个类型拥有的值范围以及并且指定日期何时间值的有效格式的描述见7.3.6 日期和时间类型。 这里是一个使用日期函数的例子。下面的查询选择了所有记录&#xff0c;其date_col的值是在最后30天以内&am…

jquery特效(1)—点击展示与隐藏全文

下班了~~~我把今天整理的一个jquery小特效发一下&#xff0c;个人觉得比较简单&#xff0c;嗖嗖的就写出来了~~~ 下面先来看最终的动态效果&#xff1a; 一、来看一下主体框架程序&#xff1a; <!DOCTYPE html> <html><head><meta charset"utf-8&quo…

python中return和printf的区别_Python格式化输出:%s和format()用法比较

1、python格式化输出历史起源python2.5之前&#xff0c;我们使用的是老式格式化输出&#xff1a;%s。从python3.0开始起(python2.6同期发布)&#xff0c;同时支持两个版本的格式化&#xff0c;多出来的一个新版本就是利用format()函数&#xff0c;进行格式化输出。2、为什么要学…