希尔伯特变换
1905年，Hilbert在研究黎曼-希尔伯特问题时提出希尔伯特变换。至于啥是黎曼-希尔伯特问题，不是专业搞数学的，我也不懂哈。1946年，Gabor定义了解析信号 y ( t ) = x ( t ) + j x ^ ( t ) y(t)=x(t)+j{{\hat x}}\left( t \right) y(t)=x(t)+jx^(t)将希尔伯特变换正式引入到信号处理领域中。此后希尔伯特变换被广泛应用到通信、机械等学科中。作为一名刚开始搞通信的机械小白，以前只知道希尔伯特变换常用在轴承故障诊断中求包络谱，在学习通信基础书籍的时候偶然看到了基带复信号的生成，才知道我原来对希尔伯特变换了解太少了，遂学习以记之。
1、从通信领域看一下hilbert变换
1.1 简单回顾一下傅里叶变换
信号是信息的载体，实际现实中存在的都是实信号，而实信号总是具有共轭对称的频谱，这个我们可以简单证明一下：
根据傅里叶变换：
X
(
ω
)
=
∫
−
∞
∞
x
(
t
)
e
−
j
ω
t
d
t
X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j{\omega}t}dt
X(ω)=∫−∞∞x(t)e−jωtdt
而根据欧拉公式：
e
−
j
ω
t
=
c
o
s
(
ω
t
)
−
j
s
i
n
(
ω
t
)
e^{-j{\omega}t}=cos({\omega}t)-jsin({\omega}t)
e−jωt=cos(ωt)−jsin(ωt)
傅里叶公式可以表示为：
X
(
ω
)
=
∫
−
∞
∞
x
(
t
)
e
−
j
ω
t
d
t
=
∫
−
∞
∞
x
(
t
)
c
o
s
(
ω
t
)
d
t
−
j
∫
−
∞
∞
x
(
t
)
s
i
n
(
ω
t
)
d
t
X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j{\omega}t}dt=\int_{-\infty}^{\infty}x(t)cos({\omega}t)dt-j\int_{-\infty}^{\infty}x(t)sin({\omega}t)dt
X(ω)=∫−∞∞x(t)e−jωtdt=∫−∞∞x(t)cos(ωt)dt−j∫−∞∞x(t)sin(ωt)dt
其实部和虚部分别为：
{
ℜ
(
X
(
ω
)
)
=
∫
−
∞
∞
x
(
t
)
cos
(
ω
t
)
d
t
ℑ
(
X
(
ω
)
)
=
∫
−
∞
∞
x
(
t
)
sin
(
ω
t
)
d
t
\left\{ \begin{array}{l} \Re \left( {X(\omega )} \right) = \int_{ - \infty }^\infty x (t)\cos (\omega t)dt\\ \Im \left( {X(\omega )} \right) = \int_{ - \infty }^\infty x (t)\sin (\omega t)dt \end{array} \right.
{ℜ(X(ω))=∫−∞∞x(t)cos(ωt)dtℑ(X(ω))=∫−∞∞x(t)sin(ωt)dt
根据上式可以看出：
{
ℜ
(
X
(
ω
)
)
=
ℜ
(
X
(
−
ω
)
)
ℑ
(
X
(
ω
)
)
=
−
ℑ
(
X
(
−
ω
)
)
\left\{ \begin{array}{l} \Re \left( {X(\omega )} \right) = \Re \left( {X( - \omega )} \right)\\ \Im \left( {X(\omega )} \right) = - \Im \left( {X( - \omega )} \right) \end{array} \right.
{ℜ(X(ω))=ℜ(X(−ω))ℑ(X(ω))=−ℑ(X(−ω))
上式说明了实函数的傅里叶变换实部是偶函数，虚部是奇函数。
对应的幅值和相位为：
{
∣
X
(
ω
)
∣
=
ℜ
2
(
X
(
ω
)
)
+
ℑ
2
(
X
(
ω
)
)
φ
(
ω
)
=
arctan
ℑ
(
X
(
ω
)
)
ℜ
(
X
(
ω
)
)
\left\{ \begin{array}{l} \left| {X(\omega )} \right| = \sqrt {{\Re ^2}\left( {X(\omega )} \right) + {\Im ^2}\left( {X(\omega )} \right)} \\ \varphi (\omega ) = \arctan \frac{{\Im \left( {X(\omega )} \right)}}{{\Re \left( {X(\omega )} \right)}} \end{array} \right.
{∣X(ω)∣=ℜ2(X(ω))+ℑ2(X(ω))φ(ω)=arctanℜ(X(ω))ℑ(X(ω))
不难得到实函数的傅里叶变换幅值是偶对称的，相位是奇对称的。
以三角函数 x = c o s ( 2 π f t + φ ) x = cos(2{\pi}ft+\varphi) x=cos(2πft+φ)为例，其幅频特性与相频特性如下图所示：
从信息的角度来看，双边频谱的信息是冗余的，可以只保留其一半的频率部分即可（一般保留正频率部分）。
1.2 通信中单边频谱的必要性
1.2.1、 提高频谱利用率
前面我们已经知道，现实中的信号都是实信号，所以总是双边频谱，但是此时信息是冗余的。在通信中，假设前面的基带信号被调制到频点fc上，即： r ( t ) = x ( t ) ∗ e x p ( j 2 π f c t ) r(t)=x(t)*exp(j2{\pi}f_ct) r(t)=x(t)∗exp(j2πfct)，此时其频谱为
可以看到原来的负频率真实占用了频谱资源。随着无线通信的发展，频谱资源愈发珍贵，这种浪费是可耻的。
1.2.2、抗混叠
在接受端，需要把调制信号从载波里提取出来，通常的做法是 ：将载频变频到零（通称为零中频）。很显然， 将接收到的实信号直接变到零中频是不行的 ，因为随着载频的下移，正、负相互接近，到中频小于信号频带一半时，两部分谱就会发生混叠，当中频为零时混叠最严重，使原信号无法恢复，具体如下所示：
2、公式推导
要想获得单边谱非常简单，只需要将实信号频谱乘以如下阶跃函数的图像，这样正频率部分保持原样输出，而负频率部分乘以0得到消除。
阶跃函数可以表示为：
ε
(
ω
)
=
{
0
,
ω
<
0
0.5
,
ω
=
0
1
,
ω
>
0
\varepsilon \left( \omega \right) = \left\{ \begin{array}{l} 0,{\rm{ }}&\omega < 0\\ 0.5,{\rm{ }}&\omega = 0\\ 1,{\rm{ }}&\omega > 0 \end{array} \right.
ε(ω)=⎩
⎨
⎧0,0.5,1,ω<0ω=0ω>0
或
ε ( ω ) = 1 2 ( 1 + s g n ( ω ) ) \varepsilon \left( \omega \right) = \frac{1}{2}\left( {1 + {\mathop{\rm sgn}} \left( \omega \right)} \right) ε(ω)=21(1+sgn(ω))
其中
s
g
n
(
ω
)
sgn(\omega)
sgn(ω)为符号函数，其定义为：
s
g
n
(
ω
)
=
{
1
ω
>
0
0
ω
=
0
−
1
ω
<
0
{\mathop{\rm sgn}} \left( \omega \right) = \left\{ \begin{array}{l} 1{\rm{ }}&\omega > 0\\ 0{\rm{ }}&\omega = 0\\ {- 1}{\rm{ }}&\omega < 0 \end{array} \right.
sgn(ω)=⎩
⎨
⎧10−1ω>0ω=0ω<0
根据傅里叶变换的性质，频域的乘积对应于时域的卷积，所以为了方便在时域进行处理，我们需要求出
ε
(
ω
)
\varepsilon \left( \omega \right)
ε(ω)或
s
g
n
(
ω
)
{\mathop{\rm sgn}} \left( \omega \right)
sgn(ω)的逆傅里叶变换。为了后面方便重复使用结果，此处以
s
g
n
(
ω
)
{\mathop{\rm sgn}} \left( \omega \right)
sgn(ω) 为例进行证明。信号与系统中常用的证明方式有微分法、指数衰减法等，为了不攀扯太多的傅里叶变换性质，我们用指数衰减法进行证明。如果单纯求其逆变换是很难求出来的：
F
−
1
(
s
g
n
(
ω
)
)
=
1
2
π
∫
−
∞
+
∞
s
g
n
(
ω
)
e
j
ω
t
d
ω
=
1
2
π
(
∫
0
+
∞
e
j
ω
t
d
ω
+
∫
−
∞
0
−
e
j
ω
t
d
ω
)
=
1
2
π
[
∫
0
+
∞
(
e
j
ω
t
−
e
−
j
ω
t
)
d
ω
]
=
j
π
[
∫
0
+
∞
sin
(
ω
t
)
d
ω
]
=
−
j
π
t
cos
(
ω
t
)
∣
0
+
∞
\begin{split} {F^{ - 1}}\left( sgn(\omega) \right) &= \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {{\mathop{\rm sgn}} \left( \omega \right)} {e^{j\omega t}}d\omega \\ &= \frac{1}{{2\pi }}\left( {\int_0^{ + \infty } {{e^{j\omega t}}} d\omega + \int_{ - \infty }^0 { - {e^{j\omega t}}} d\omega } \right)\\ &= \frac{1}{{2\pi }}\left[ {\int_0^{ + \infty } {\left( {{e^{j\omega t}} - {e^{ - j\omega t}}} \right)} d\omega } \right]\\ &= \frac{j}{\pi }\left[ {\int_0^{ + \infty } {\sin \left( {\omega t} \right)} d\omega } \right]\\ &= \frac{{ - j}}{{\pi t}}\cos \left( {\omega t} \right)\left| {_0^{ + \infty }} \right. \end{split}
F−1(sgn(ω))=2π1∫−∞+∞sgn(ω)ejωtdω=2π1(∫0+∞ejωtdω+∫−∞0−ejωtdω)=2π1[∫0+∞(ejωt−e−jωt)dω]=πj[∫0+∞sin(ωt)dω]=πt−jcos(ωt)∣
∣0+∞
这里由于无法确定
c
o
s
(
∞
t
)
cos({\infty}t)
cos(∞t)等于多少，所以求解无法进一步进行。这种信号由于无法保证能量有限，所以很难直接求出其傅里叶变换式，往往需要增加一个指数衰减项
e
−
a
ω
(
a
>
0
)
e^{-a{\omega}}(a>0)
e−aω(a>0)来保证能量有限，才能继续计算：
F
−
1
(
s
g
n
(
ω
)
)
=
1
2
π
[
∫
0
+
∞
e
−
a
ω
(
e
j
ω
t
−
e
−
j
ω
t
)
d
ω
]
=
1
2
π
∫
0
+
∞
e
(
−
a
+
j
t
)
ω
−
e
(
−
a
−
j
t
)
ω
d
ω
=
1
2
π
[
(
1
−
a
+
j
t
)
e
(
−
a
+
j
t
)
ω
−
(
1
−
a
−
j
t
)
e
(
−
a
−
j
t
)
ω
]
∣
0
+
∞
\begin{split} {F^{ - 1}}\left( sgn(\omega) \right) &= \frac{1}{{2\pi }}\left[ {\int_0^{ + \infty } e^{-a{\omega}}{\left( {{e^{j\omega t}} - {e^{ - j\omega t}}} \right)} d\omega } \right]\\ &= \frac{1}{{2\pi }}{\int_0^{ + \infty } e^{(-a+jt)\omega}-e^{(-a-jt)\omega} d\omega } \\ &=\frac{1}{{2\pi }}[(\frac{1}{-a+jt})e^{(-a+jt)\omega}-(\frac{1}{-a-jt})e^{(-a-jt)\omega}]|_0^{+\infty} \end{split}
F−1(sgn(ω))=2π1[∫0+∞e−aω(ejωt−e−jωt)dω]=2π1∫0+∞e(−a+jt)ω−e(−a−jt)ωdω=2π1[(−a+jt1)e(−a+jt)ω−(−a−jt1)e(−a−jt)ω]∣0+∞
此时，由于
lim
ω
→
∞
e
−
a
ω
=
0
\mathop {\lim }\limits_{\omega \to \infty } {e^{ - a\omega }} = 0
ω→∞lime−aω=0，并且假设
a
a
a是一个趋于0的数，上式可以进一步化简为：
F
−
1
(
s
g
n
(
ω
)
)
=
1
2
π
lim
a
→
∞
[
0
−
(
1
−
a
+
j
t
−
1
−
a
−
j
t
)
]
=
1
2
π
−
2
j
t
=
j
π
t
\begin{split} {F^{ - 1}}\left( sgn(\omega) \right) &= \frac{1}{{2\pi }}\mathop {\lim }\limits_{a \to \infty } \left[ {0 - \left( {\frac{1}{{ - a + jt}} - \frac{1}{{ - a - jt}}} \right)} \right]\\ &= \frac{1}{{2\pi }}\frac{{ - 2}}{{jt}}\\ &= \frac{j}{{\pi t}} \end{split}
F−1(sgn(ω))=2π1a→∞lim[0−(−a+jt1−−a−jt1)]=2π1jt−2=πtj
同样的
e
−
a
ω
(
a
>
0
)
e^{-a{\omega}}(a>0)
e−aω(a>0)技巧，也可以直接得到
ε
(
ω
)
\varepsilon \left( \omega \right)
ε(ω)的逆傅里叶变换(这里不再证明了，大家自己尝试证一下)或者利用性质：
F
−
1
(
2
ε
(
ω
)
)
=
F
−
1
(
1
+
s
g
n
(
ω
)
)
=
F
−
1
(
1
)
+
F
−
1
(
s
g
n
(
ω
)
)
=
δ
(
t
)
+
j
π
t
\begin{split} {F^{ - 1}}\left( {2\varepsilon \left( \omega \right)} \right) &= {F^{ - 1}}\left( {1 + {\rm{sgn}}\left( \omega \right)} \right)\\ &= {F^{ - 1}}\left( 1 \right) + {F^{ - 1}}\left( {{\rm{sgn}}\left( \omega \right)} \right)\\ &= \delta \left( t \right) + \frac{j}{{\pi t}} \end{split}
F−1(2ε(ω))=F−1(1+sgn(ω))=F−1(1)+F−1(sgn(ω))=δ(t)+πtj
由此我们可以直接时域卷积得到所需要的单边带信号：
x
+
(
t
)
=
x
(
t
)
∗
[
δ
(
t
)
+
j
π
t
]
=
x
(
t
)
+
j
x
(
t
)
∗
1
π
t
\begin{array}{c} {x_ + }\left( t \right) = x\left( t \right) * \left[ {\delta \left( t \right) + \frac{j}{{\pi t}}} \right]\\ = x\left( t \right) + jx\left( t \right) * \frac{1}{{\pi t}} \end{array}
x+(t)=x(t)∗[δ(t)+πtj]=x(t)+jx(t)∗πt1
这里
x
(
t
)
∗
1
π
t
x\left( t \right) * \frac{1}{{\pi}t}
x(t)∗πt1就被定义为希尔伯特变换，根据卷积的定义，
x
(
t
)
x(t)
x(t)的希尔伯特变换
x
^
(
t
)
\hat x(t)
x^(t)为：
x
^
(
t
)
=
x
(
t
)
∗
1
π
t
=
1
π
∫
−
∞
+
∞
x
(
τ
)
t
−
τ
d
τ
\hat x\left( t \right) = x\left( t \right) * \frac{1}{{\pi t}} = \frac{1}{\pi }\int_{ - \infty }^{ + \infty } {\frac{{x\left( \tau \right)}}{{t - \tau }}} d\tau
x^(t)=x(t)∗πt1=π1∫−∞+∞t−τx(τ)dτ
其中，
h
(
t
)
=
1
π
t
h\left( t \right) = \frac{1}{{\pi t}}
h(t)=πt1是该系统的单位冲击响应。中学我们应该学过，这个冲激响应类似一个反比例函数，在t=0处是没有定义的，一般认为h(0)=0。
常用的希尔伯特变换对如下表所示：
3、希尔伯特变换的性质
性质比较多，我们只在下面列举几个比较常用的性质。
3.1 90°移相器
在以前机械的时候，老师总是讲希尔伯特变换就是一个90度移相器，为啥这么说呢？我们可以求一下h(t)的傅里叶变换就一目了然了：
根据前面已有的证明，我们知道傅里叶变换对：
x
(
t
)
=
j
π
t
↔
X
(
ω
)
=
s
g
n
(
ω
)
x\left( t \right) = \frac{j}{{\pi t}} \leftrightarrow X\left( \omega \right) = {\mathop{\rm sgn}} \left( \omega \right)
x(t)=πtj↔X(ω)=sgn(ω)
则:
h
(
t
)
=
1
π
t
↔
H
(
ω
)
=
−
j
s
g
n
(
ω
)
=
{
−
j
(
ω
>
0
)
0
(
ω
=
0
)
j
(
ω
<
0
)
h\left( t \right) = \frac{1}{{\pi t}} \leftrightarrow H\left( \omega \right) = - j{\rm{sgn}}\left( \omega \right) = \left\{ \begin{array}{l} {- j}\left( {\omega > 0} \right)\\ 0\left( {\omega = 0} \right)\\ j\left( {\omega < 0} \right) \end{array} \right.
h(t)=πt1↔H(ω)=−jsgn(ω)=⎩
⎨
⎧−j(ω>0)0(ω=0)j(ω<0)
可以看到除直流外，模值
H
(
ω
)
=
1
H(\omega)=1
H(ω)=1，相位对正频率偏移-π/2，对负频率偏移π/2。这就是90度移相器的由来。
3.2 多次希尔伯特变换
下图演示了多次进行希尔伯特的变化，正频率部分每次偏移-π/2，负频率偏移π/2。
理论上来说，多次H变换时，有：
H
n
(
x
(
t
)
)
↔
(
−
j
s
g
n
(
ω
)
)
n
F
(
ω
)
H^n(x(t)) \leftrightarrow(-jsgn(\omega))^nF(\omega)
Hn(x(t))↔(−jsgn(ω))nF(ω)
两次变换相当于对原始取反：
H 2 ( ω ) = { − 1 ( ω ≠ 0 ) 0 ( ω = 0 ) H^2\left( \omega \right) = \left\{ \begin{array}{l}\ { - 1}\left( {\omega}{\ne} 0 \right)\\ 0\left( {\omega = 0} \right) \end{array} \right. H2(ω)={ −1(ω=0)0(ω=0)
四次变换等于本身:
H 4 ( ω ) = { 1 ( ω ≠ 0 ) 0 ( ω = 0 ) H^4\left( \omega \right) = \left\{ \begin{array}{l} 1\left( {\omega}{\ne} 0 \right)\\ 0\left( {\omega = 0} \right) \end{array} \right. H4(ω)={1(ω=0)0(ω=0)
三次变换时：
H
3
(
x
(
t
)
)
↔
(
−
j
s
g
n
(
ω
)
)
3
F
(
ω
)
=
−
j
F
(
ω
)
{H^3}(x(t)) \leftrightarrow {( - jsgn(\omega ))^3}F(\omega ) = - jF(\omega)
H3(x(t))↔(−jsgn(ω))3F(ω)=−jF(ω)
4、如何实现希尔伯特变换
4.1 FFT变换法
对于一离散函数
u
[
n
]
u[n]
u[n] ，以及其离散傅里叶变换函数
U
(
ω
)
U(\omega)
U(ω) ，可推得其希尔伯特变换为:
H
(
u
)
[
n
]
=
D
T
F
T
−
1
{
U
(
ω
)
⋅
σ
H
(
ω
)
}
H(u)[n]=D T F T^{-1}\left\{U(\omega) \cdot \sigma_H(\omega)\right\}
H(u)[n]=DTFT−1{U(ω)⋅σH(ω)}
其中
σ
H
(
ω
)
=
{
j
,
−
π
<
ω
<
0
−
j
,
0
<
ω
<
π
0
,
ω
=
−
π
,
0
,
π
\sigma_H(\omega) {=} \begin{cases}j, & -\pi<\omega<0 \\ -j, & 0<\omega<\pi \\ 0, & \omega=-\pi, 0, \pi\end{cases}
σH(ω)=⎩
⎨
⎧j,−j,0,−π<ω<00<ω<πω=−π,0,π
所以实现希尔伯特变换可以从频域入手：
- 对输入向量x进行n点FFT变换;
- 创建如下的向量 h(i):
- 0 当
i= 1,(n/2)+1- -j 当
i= 2, 3, … ,(n/2)- j 当
i=(n/2)+2, … ,n- 频域进行乘积FFT(x)*h；
- 计算FFT(x)*h的IFFT，即为最终的结果.
附上一段MATLAB实现的代码：
function x = myhilbert(xr,n)
if nargin<2, n=[]; end
if ~isreal(xr)
warning(message('signal:hilbert:Ignore'))
xr = real(xr);
end
% Work along the first nonsingleton dimension
[xr,nshifts] = shiftdim(xr);
if isempty(n)
n = size(xr,1);
end
x = fft(xr,n,1); % n-point FFT over columns.
h = 1j*ones(n,~isempty(x),'like',x); % nx1 for nonempty. 0x0 for empty.
if n > 0 && 2*fix(n/2) == n
% even and nonempty
h([1 n/2+1]) = 0;
h(2:n/2) = -1j;
elseif n>0
% odd and nonempty
h(1) = 0;
h(2:(n+1)/2) = -1j;
end
x = ifft(x.*h(:,ones(1,size(x,2))),[],1);
% Convert back to the original shape.
x = shiftdim(x,-nshifts);
end
这个代码中有两点需要解释一下:
1、为啥 h ( i ) h(i) h(i)和 σ H ( ω ) \sigma_H(\omega) σH(ω)不对应？
这个是MATLAB计算FFT的原因，经过FFT之后，MATLAB输出的频率范围是[0,fs]，fs是采样率，前面 σ H ( ω ) \sigma_H(\omega) σH(ω)的对应范围是[-fs/2,fs/2]，也就是零频在中间，所以排列0,-j,j对应移动完成这个过程，本质上还是对应的。
2、为啥最终输出结果和MATLAB中hilbert函数输出结果不一样？
MATLAB中输出的解析信号 y ( t ) = x ( t ) + x ^ ( t ) y(t)=x(t)+{{\hat x}}\left( t \right) y(t)=x(t)+x^(t)，其构建的 h ( i ) h(i) h(i)为
创建如下的向量 h(i):
- 1 当
i= 1,(n/2)+1- 2 当
i= 2, 3, … ,(n/2)- 0 当
i=(n/2)+2, … ,n
也就是对应我们前面推导的 ε ( ω ) \varepsilon \left( \omega \right) ε(ω)。
4.2 希尔伯特滤波器
此外，根据卷积定律，另一个相等的方程式为:
H
(
u
)
[
n
]
=
u
[
n
]
∗
h
[
n
]
H(u)[n]=u[n] * h[n]
H(u)[n]=u[n]∗h[n]
其中
h
[
n
]
=
D
T
F
T
−
1
{
σ
H
(
ω
)
}
=
{
0
,
n为偶数
2
π
n
n为奇数
h[n] \stackrel{}{=} D T F T^{-1}\left\{\sigma_H(\omega)\right\}= \begin{cases}0, & \text n为偶数 \\ \frac{2}{\pi n} & \text n为奇数\end{cases}
h[n]=DTFT−1{σH(ω)}={0,πn2n为偶数n为奇数
所以可以通过
h
(
n
)
h(n)
h(n)的定义直接生成一个滤波器系数，不过由于实际滤波器阶数有限，直接生成的
h
(
n
)
h(n)
h(n)并没有MTATLAB采用最下化加权积分平方误差的效果好。下图给出了采用MATLAB中使用firls和直接根据公式生成的幅频特性，可以看出firls生成的通带波动更小一点。
——未完待续——
欢迎关注我的微信公众号
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
