本帖最后由 会笑的星星 于 2021-7-23 14:23 编辑
最后两篇**写希尔伯特变换与小波变换，这两个常见的变换用的也很普遍，他们处理的信号与fft不同，fft的频谱图只是展现时域信号中具有的频率成分，但我们需要知道什么时间点出现了一些什么频率时，fft有效性就大打折扣了。而希尔伯特变换与小波变换就是为此而生，他不仅仅知道时域信号的频谱，还是知道什么时间点出现了哪一个频率。坦白说，我没有系统的学习过这希尔伯特变换与小波变换，只知道他们的基本原理。但好在因为有python，有高人已经写好了现成的代码，可以让我们用来仿真分析，从而能对这些变换有一定的认识，等到需要时，知道有些工具可以拿来用，如果真的有用，到时候再去学习它也不迟。本篇**先介绍一下用python做希尔伯特变换以及希尔伯特黄变换。如果你对这两个变换还没有认识的话，建议你先看看[1][2]。
由于希尔伯特-黄变换是建立在希尔伯特变换之上，所以，我们先说希尔伯特变换。希尔伯特变换实际上就是一个滤波器，他的作用是把输入信号的相位推迟π/2，要注意，输入信号与输出信号的频率依然是一样的。这有什么用呢？假设x（t）是输入信号，y（t）是x（t）的希尔伯特变换，构造一个信号z（t） = x（t） + i*y（t），显然，z（t）与x（t）、y（t）的频率是一样的（假设x（t）、y（t）的周期为T，那么z(t+T)=z(t)）。通过z（t），我们可以知道该信号的瞬时相位为ψ（t）= arctan（y（t）/x（t））。又根据频率定义，我们知道相位对时间的变化率就是频率，即 f = (1/2π)*dψ（t）/dt。由此，我们便知道了这个连续信号的瞬时频率。我们又知道，从连续域的信号转换为离散域信号，其离散信号频域会扩展T倍（T=1/fs,fs为采样率）。由此，对于离散信号而言，其瞬时频率f = T *（1/2π）*（dψ（t）/dt） = (1/(2π*fs))*（dψ（t）/dt）。有了以上分析，我们就可以看下面的python代码了。
# -*- coding: utf-8 -*-
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
duration = 1.0 #持续时间，这个是1s的意思
fs = 400.0
samples = int(fs*duration)
t = np.arange(samples) / fs
#产生一个啁啾信号（单频且频率随时间的增加而增加）
#t为时间，20为起始频率，后面两个参数t[-1]、100表示信号信号最后点对应的频率是100
x = signal.chirp(t, 20.0, t[-1], 100.0)
#signal.hilbert()对信号x做希尔伯特变换变换，函数原型见[3]
ht_s = signal.hilbert(x)
#np.angle()求瞬时相位
#np.unwrap()函数是为了防止相位卷绕[4],如果不理解也没关系,加上他是为了防止
#相位计算错误
instantaneous_phase = np.unwrap(np.angle(ht_s))
#np.diff()求解相位的导数
instantaneous_frequency = (np.diff(instantaneous_phase) /
(2.0*np.pi) * fs)
plt.subplot(2,1,1)
plt.ylabel('chirp signal')
plt.plot(t,x)
plt.subplot(2,1,2)
plt.ylabel('instantaneous frequency')
plt.plot(t[0:-1],instantaneous_frequency)
plt.show()
可以看到上述频谱图中，随着时间增加，频率也在增加，这是符合实际信号特征的。要注意的是，在时间的两端（0s、1s）附近，得到的频谱是不稳定的，而在0~1s之间的频率则是准确且稳定的。因此，通过希尔伯特变换，我们是可以找出时域信号中的频率与时间关系的。
上述我们做希尔伯特变换的信号是啁啾信号，如果我把信号换成由两个频率合成的信号，如x = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*100*t)，他的希尔伯特变换后的频谱图又是怎么样的呢？由于代码与上面所述代码雷同，这里不再重复，直接给出他的频谱图如下。
可以看到，对于合成信号，希尔伯特变换构成的新信号的频谱并不能反应真实的信号频谱。如上图所示，频率已经完全变化了。从这里也可以看出，企图直接用希尔伯特变换找出复合信号的频谱图是不可行的。那怎么办呢？希尔伯特变-黄变换就呼之欲出了。他可以把合成信号分解成若干个单独的能做希尔伯特变换的信号（学名叫做经验模态分解，简称EMD，每个被分解的信号简称imf），从而就知道了任何时刻信号中有哪些频率，见下面代码：
# -*- coding: utf-8 -*-
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import pyhht as hht
duration = 1.0
fs = 1000
samples = int(fs*duration)
t = np.arange(samples) / fs
x = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*50*t)
#这里用到了一个新模块pyhht，该模块专门用于做希尔伯特-黄变换，详见[5]
#hht.EMD()实例化一个对象
decomposer = hht.EMD(x)
#分解这个信号，得到若干组能做希尔伯特变换的信号
imfs = decomposer.decompose();
#绘制各组分解信号
hht.plot_imfs(x,imfs)
上图中，signal是原始信号，imf1是被分解的一组信号，res表示残差，也就是原始信号中剥离imf1后的剩余信号，在这里，singla = imf1 + res。
然后，我们对decomposer.decomposer（）返回的两组信号做希尔伯特变换，完整代码如下：
# -*- coding: utf-8 -*-
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import pyhht as hht
duration = 1.0
fs = 1000
samples = int(fs*duration)
t = np.arange(samples) / fs
x = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*50*t)
decomposer = hht.EMD(x)
imfs = decomposer.decompose();
#hht.plot_imfs(x,imfs)
analytic_signal = signal.hilbert(imfs)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
freq = (np.diff(instantaneous_phase) / (2.0*np.pi) * fs)
plt.ylim(0,500)
plt.plot(t[0:-1], freq[0])
plt.plot(t[0:-1], freq[1])
plt.show()
得到的频谱图如下：
可以看出，希尔伯特-黄变换通过分解信号，然后对每个信号做希尔伯特变换，得到了合成信号的频谱，这样就解决了希尔伯特变换无法得到合成信号的频谱问题。
针对有些时域信号，使用的上述的EMD算法可能无法很好的分解信号[7]，但有改进型的EMD算法可供使用，比如EEMD、CEEMDAN[6]，你可以根据情况使用或者测试不同的算法。
[1]https://www.zhihu.com/question/30372795
[2]https://blog.csdn.net/qq_42688495/article/details/106961315
[3]https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.signal.hilbert.html#scipy.signal.hilbert
[4]https://blog.csdn.net/lishuhuakai/article/details/78812540
[5]https://pyhht.readthedocs.io/en/latest/apiref/pyhht.html
[6]https://pyemd.readthedocs.io/en/latest/eemd.html
[7]https://blog.csdn.net/liu_xiao_cheng/article/details/83897034
|
楼主
标题置顶
标题高亮
