信号特征之希尔伯特变换（Python、C++、MATLAB实现）

希尔伯特变换

1 特征描述

希尔伯特变换广泛使用于信号处理应用中，以获得信号的解析表示。其可以计算瞬时频率和相位，相位被定义为原始信号和信号的希尔伯特变换之间的角度。对于信号$x(t)$，其希尔伯特变换定义如下。

$$\hat{x}(t)=H[x(t)]=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{x(\tau )}{t-\tau }\text{d}\tau$$

由上式可知，随着变换的结果，自变量不变，因此输出$\hat{x}(t)$也是与时间有关的函数。此外，$\hat{x}(t)$是$x(t)$的线性函数。它是由${(\pi t)}^{-1}$与$x(t)$卷积获得的，如下关系所示。

$$\hat{x}(t)=\frac{1}{\pi t}\ast x(t)$$

由上式可以得到$x(t)$的希尔伯特变换的傅立叶变换如下所示。

$$\begin{array}{rl}F[\hat{x}(t)]& =F[\frac{1}{\pi t}]F[x(t)]\\ & =-jsgn(f)F[x(t)]\end{array}$$

其中，$sgn(f)$如下所示。

$$sgn(f)=\{\begin{array}{l}\begin{array}{rl}+1,& f>0\\ 0,& f=0\\ -1,& f<0\end{array}\end{array}$$

希尔伯特变换对实际数据作90度相移，正弦变为余弦，反之亦然。希尔伯特变换可用于从实际信号中产生解析信号，解析信号在通信领域中很有用，尤其是在带通信号处理中。解析信号如下所示。

$$s(t)=x(t)+j\hat{x}(t)$$

解析信号的实部是原始的信号数据，虚部是实际的希尔伯特变换结果。

2 数据来源

提供一串信号的数据:
这组数据存在三个频率分量，分别为100Hz，200Hz和300Hz, 用MATLAB生成。

# 定义三个频率分量的参数
t = 0:0.0001:2/100;

freq1 = 100; # 第一个频率分量的频率(Hz)
freq2 = 200; # 第二个频率分量的频率(Hz)
freq3 = 300; # 第三个频率分量的频率(Hz)

# 生成三个频率分量的正弦波信号
signal1 = sin(2 * pi * freq1 * t);
signal2 = sin(2 * pi * freq2 * t);
signal3 = sin(2 * pi * freq3 * t);

# 将三个频率分量的信号相加
result = signal1 + signal2 + signal3;

3 函数的Python代码

3.1 Python库要求

from scipy import signal

3.2 hilbert希尔伯特变换函数

其中，inputS为输入的信号序列，函数的返回值outputS为解析信号序列。解析信号的实部是原始的信号数据，虚部是实际的希尔伯特变换结果。

outputS = signal.hilbert(inputS)

3.3 希尔伯特变换函数hilbert验证

from scipy import signal

inputS = input().split()

outputS = signal.hilbert(inputS)

for i in outputS:
    print("{} ".format(i.imag), end='')

4 函数的C++代码

4.1 C++库要求

#include<iostream>
#include<complex>
#include<vector>
#include<algorithm>

4.2 void HilbertTransform(vector<double> &cinS, vector<complex<double>> &outputS)希尔伯特变换函数

其中，cinS为输入的信号序列，outputS为得到的解析信号序列。解析信号的实部是原始的信号数据，虚部是实际的希尔伯特变换结果。

void HilbertTransform(vector<double> &cinS, vector<complex<double>> &outputS) {

    //flag=-1时为正变换, flag=1时为反变换
    auto dft = [](vector<complex<double>> &Data, int flag) {

        int length = Data.size();
        complex<double> wn;
        vector<complex<double>> temp(length);
        
        for (size_t i = 0; i < length; i++) {
            temp[i] = 0.;
            for (size_t j = 0; j < length; j++) {
                wn = complex<double>(cos(2. * acos(-1) / length * i * j), sin(flag * 2. * acos(-1) / length * i * j));
                temp[i] = temp[i] + Data[j] * wn;
            }
        }

        if (flag == -1) {
            Data = temp;
        } else if (flag == 1) {
            for (size_t i = 0; i < length; i++) {
                Data[i] = temp[i] / double(length);
            }
        }
        return;
    };
    int length = cinS.size();
    outputS.resize(length);

    for (size_t i = 0; i < length; i++) {
        outputS[i].real(cinS[i]);
        outputS[i].imag(0.);
    }
    dft(outputS, -1); //DFT
    int half = ((length % 2) == 0) ? (length / 2) : ((length + 1) / 2);
    for (size_t i = half; i < length; i++) {
        outputS[i] = 0.;
    }
    dft(outputS, 1); //IDFT
    for (size_t i = 0; i < length; i++) {
        outputS[i] = 2. * outputS[i];
    }
    return;
}

4.3 希尔伯特变换函数void HilbertTransform(vector<double> &cinS, vector<complex<double>> &outputS)验证

int main() {

    double temp;
    vector<double> cinS;
    vector<complex<double>> outputS;

    int cinNum;
    cin >> cinNum; //201

    for (size_t i = 0; i < cinNum; i++) {
        cin >> temp;
        cinS.emplace_back(temp);
    }

    HilbertTransform(cinS, outputS);
    
    cout << endl << endl;
    for (size_t i = 0; i < cinNum; i++) {
        cout << imag(outputS[i]) << " ";
    }

    return 0;
}

5 函数的MATLAB代码

5.1 hilbert希尔伯特变换函数

其中，result为输入的信号序列，函数的返回值HilbertR为解析信号序列。解析信号的实部是原始的信号数据，虚部是实际的希尔伯特变换结果。

HilbertR = hilbert(result);

5.2 希尔伯特变换函数hilbert验证

HilbertR = hilbert(result);
plot(t, real(HilbertR))
hold on
plot(t, imag(HilbertR))
hold off
legend('Real Part','Imaginary Part')

5.3 hilbert希尔伯特变换函数原理解析

MATLAB中的希尔伯特变换并没有直接给出信号的变换结果，而是得到了解析信号。解析信号的实部是原始的信号数据，虚部是实际的希尔伯特变换结果。

事实上，MATLAB是先通过fft得到原始信号的快速傅立叶变换，再将fft后的后半部分信号置零，最后通过ifft得到了解析信号。

5.4 希尔伯特变换函数hilbert原理验证

下面对上述信号处理过程进行验证。

fftR = fft(result);
if mod(length(fftR), 2) == 0
    half = length(fftR) / 2;
else
    half = (length(fftR) + 1) / 2;
end
fftR(half + 1:length(fftR)) = 0;
ifftfftR = 2 * ifft(fftR);

通过下面的绘制可以发现上述处理过程与直接使用hilbert函数进行希尔伯特变换得到的结果相同。

plot(t, real(HilbertR))
hold on
plot(t, real(ifftfftR))
hold off
figure
plot(t, imag(HilbertR))
hold on
plot(t, imag(ifftfftR))
hold off

上文中的C++代码就是根据此处理过程编写的，不过使用的是DFT而非FFT。

6 结果

6.1 Python绘制的希尔伯特变换图像

6.2 C++绘制的希尔伯特变换图像

6.3 MATLAB绘制的希尔伯特变换图像