信号特征之希尔伯特变换（Python、C++、MATLAB实现）
Vistyxio 2024-07-17 12:35:06 阅读 65
希尔伯特变换
1 特征描述
希尔伯特变换广泛使用于信号处理应用中，以获得信号的解析表示。其可以计算瞬时频率和相位，相位被定义为原始信号和信号的希尔伯特变换之间的角度。对于信号
x
(
t
)
x(t)
x(t)，其希尔伯特变换定义如下。
x
^
(
t
)
=
H
[
x
(
t
)
]
=
1
π
∫
−
∞
+
∞
x
(
τ
)
t
−
τ
d
τ
\hat{x}(t)=H[x(t)]=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{t-\tau}\,\text{d}\tau
x^(t)=H[x(t)]=π1∫−∞+∞t−τx(τ)dτ
由上式可知，随着变换的结果，自变量不变，因此输出
x
^
(
t
)
\hat{x}(t)
x^(t)也是与时间有关的函数。此外，
x
^
(
t
)
\hat{x}(t)
x^(t)是
x
(
t
)
x(t)
x(t)的线性函数。它是由
(
π
t
)
−
1
{({\pi}t)}^{-1}
(πt)−1与
x
(
t
)
x(t)
x(t)卷积获得的，如下关系所示。
x
^
(
t
)
=
1
π
t
∗
x
(
t
)
\hat{x}(t)=\frac{1}{ {\pi}t}\,\ast\,x(t)
x^(t)=πt1∗x(t)
由上式可以得到
x
(
t
)
x(t)
x(t)的希尔伯特变换的傅立叶变换如下所示。
F
[
x
^
(
t
)
]
=
F
[
1
π
t
]
F
[
x
(
t
)
]
=
−
j
s
g
n
(
f
)
F
[
x
(
t
)
]
\begin{split} F[\hat{x}(t)]&=F[\frac{1}{ {\pi}t}]F[x(t)] \\ &=-j\,sgn(f)F[x(t)] \end{split}
F[x^(t)]=F[πt1]F[x(t)]=−jsgn(f)F[x(t)]
其中，
s
g
n
(
f
)
sgn(f)
sgn(f)如下所示。
s
g
n
(
f
)
=
{
+
1
,
f
>
0
0
,
f
=
0
−
1
,
f
<
0
sgn(f)= \begin{cases} \begin{split} +1,&{\quad}f>0\\ 0,&{\quad}f=0\\ -1,&{\quad}f<0 \end{split} \end{cases}
sgn(f)=⎩
⎨
⎧+1,0,−1,f>0f=0f<0
希尔伯特变换对实际数据作90度相移，正弦变为余弦，反之亦然。希尔伯特变换可用于从实际信号中产生解析信号，解析信号在通信领域中很有用，尤其是在带通信号处理中。解析信号如下所示。
s
(
t
)
=
x
(
t
)
+
j
x
^
(
t
)
s(t)=x(t)+j\hat{x}(t)
s(t)=x(t)+jx^(t)
解析信号的实部是原始的信号数据，虚部是实际的希尔伯特变换结果。
2 数据来源
提供一串信号的数据:
这组数据存在三个频率分量，分别为100Hz，200Hz和300Hz, 用MATLAB生成。
<code># 定义三个频率分量的参数
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='')code>
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绘制的希尔伯特变换图像
声明
本文内容仅代表作者观点，或转载于其他网站，本站不以此文作为商业用途
如有涉及侵权，请联系本站进行删除
转载本站原创文章，请注明来源及作者。