1、Hilbert解调原理

我们知道信号解调是信号调制的反过程。解调就是从已调制的高频信号中解调出原始调制信号。信号调制包括调幅、调频、调相，因此信号解调的目的就是：根据已有信号，提取出信号中的包络、相位、频率信息。

Hilbert变换的作用是把信号频率分量的相位推迟90度，因此也叫做90度移相器。下面我们看一下通过Hilbert变换是如何实现信号解调的。

首先，我们假设有一个调制后的信号，形式如下所示。其中A(t)是幅度调制信息，fn是载波频率，ψ(t)是相位调制信息。

信号x(t)的希尔伯特表达式，可以用下式来表达。y(t)和x(t)相位差90度，相当于余弦变成正弦。

我们用下式来表示x(t)，y(t)的相位。

用式（一）和式（二）左右两边平方、相加再开根号可得瞬时幅值，如下所示：

用式（二）除以式（一），两端同时取反正切可得瞬时相位，如下所示：

对式（三）两端求导数可得瞬时频率（单位为频率），如下所示：

上述求式（四）、式（五）、式（六）的过程有什么意义呢？在信号解析过程中，我们实际上只知道x(t)信号，A(t)、fn和ψ(t)都可能是未知的。我们可以通过hilbert算法得到y(t) = hilbert(x(t))。那么上面的解调过程可以理解为，已知x(t)和y(t)得出瞬时幅值、瞬时相位和瞬时频率。

2、Hilbert解调仿真

在Matlab软件中，y=hilbert(x)，其中x表示输入信号，输出信号y是一个复数序列，y的实部是原始实数序列x，虚部是Hilbert变换的结果。

本次仿真的目的是构造一个调制信号，用三种方法计算其瞬时幅值、相位、频率。第一种方法是直接使用matlab内置的算法；第二种方法是直接求解方程，解出瞬时幅值、相位、频率。这种方法只能在已知调制信号的场景下使用，用于事后分析。第三种方法是使用本文中用到的数学推导公式来计算瞬时幅值、相位、频率。如果三种方法结论一致，那就是最好的结果。

%基本参数配置
fs = 1000;
T = 1;
N = fs * T;
dt = 1 / fs;
t = (0:N-1) * dt;

%调制信号设置，a是调幅，b是调相，c是载波
a = 1 + 0.5 * cos(2 * pi * 5 * t);
b = sin(2 * pi * 10 * t);
c = cos(2 * pi * 30 * t + b);
x = a .* c;
y = hilbert(x);

figure(1),
hold on
plot(t,real(y),'red');
plot(t,imag(y),'blue');
plot(t,abs(y),'green');
hold off

%方法1，用matlab内置算法得到的瞬时幅值、相位、频率
z1 = abs(y);
z2 = unwrap(angle(y));
z3 = instfreq(x, fs,'Method','hilbert')';

figure(2),
subplot(3, 1, 1); 
plot(t,z1,'red');
subplot(3, 1, 2); 
plot(t,z2,'green');
subplot(3, 1, 3); 
plot(t(1:N-1),z3,'blue');
 
%方法2，求解x(t)，y(t)得到的瞬时幅值、相位、频率
m1 = a;
m2 = 2*pi*30*t + sin(2*pi*10*t);
m3 = 2*pi*30 + 20*pi*cos(2*pi*10*t);

figure(3),
subplot(3, 1, 1); 
plot(t,m1,'red');
subplot(3, 1, 2); 
plot(t,m2,'green');
subplot(3, 1, 3); 
plot(t,m3,'blue');
 
%方法3，根据定义得到的瞬时幅值、相位、频率
p1 = sqrt(real(y).^2 + imag(y).^2);
p2 = atan(imag(y)./real(y));
p3 = diff(p2);

figure(4),
subplot(3, 1, 1); 
plot(t,p1,'red');
subplot(3, 1, 2); 
plot(t,p2,'green');
subplot(3, 1, 3); 
plot(t,p3,'blue');

• 图1表示hilbert变换的实部、虚部、包络线之间的关系。

• 图2使用matlab内置算法画出了调制信号的瞬时幅值、瞬时相位、瞬时频率。

• 图3我们对调制信号直接求导可以算出它的瞬时幅值、瞬时相位、瞬时频率。
• 我们看到图3的瞬时幅值、瞬时相位和图2完全相同，瞬时频率形态相同，但是纵坐标有差异。

• 图4我们用推导出的公式计算信号的瞬时幅值、瞬时相位、瞬时频率。
• 图4的瞬时幅值和图2，图3完全相同，但是瞬时相位、瞬时频率完全不同，差异很大。

这个仿真结果有些不理想，推导公式应该是没有问题，但是画出来的图差异却非常大，这个结果有些令人费解（后续会继续找原因）。

3、Hilbert解调应用1

本节我们用小波律动公司最近刚接收到的一组早期预警数据进行分析。振动传感器安装在一个齿轮箱上，其中包含多个轴承和齿轮，输入轴转速大约为1135r/min，齿轮箱内部参数都已知，所以特征频率都能计算出来，这里不一一列举了。数据的采样频率为2000Hz，能覆盖所有的轴承特征频率，但是不能覆盖齿轮啮合频率，所以这里是存在频率混叠现象的，可以参考之前发的文章《振动频谱分析中的频率混叠现象》。

预警信息内容如下所示：

特征频率:(100.586Hz,0.795)出现了2倍频(201.172Hz,0.146);角接触球轴承162250LB滚动体过外圈频率89.625Hz 出现了边频对(62.5Hz，0.363)、(119.141Hz，0.255);

这条预警信息中包含两个内容。

1、(100.586Hz,0.795)出现了2倍频(201.172Hz,0.146)，但是这个特征频率没有匹配到，可能是某个高频频率混叠后变成了100.586Hz。这个问题我们今天不分析。

2、角接触球轴承162250LB滚动体过外圈频率89.625Hz 出现了边频对(62.5Hz，0.363)、(119.141Hz，0.255)。这个是我们今天要分析的重点，我们用python代码来进行分析。

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft
from scipy.signal import hilbert

#载入数据文件，生成数值数组
f = open("/Users/tom/Desktop/123456.txt", encoding = "utf-8")
fstr = f.read()
f.close()
flist = fstr.split()
x = []
for key in flist:
	x.append(float(key))

#基本信息设置
fs = 2000
N = len(x)
t = np.arange(0, N/fs, 1/fs)
t1 = np.arange(0, fs/2, fs/N)

#fft变换 + hilbert变换
f1 = fft(x)
yabs = abs(f1)
y = hilbert(x)
h = abs(y)
f2 = fft(h)
fabs2 = abs(f2)
fabs2[0] = 0;#去掉初始相位的影响

#画原始波形图
plt.figure()
plt.subplot(311)
plt.plot(t, x)

#画原始波形的频谱图，频谱图上画出一些点，这些点需要额外计算
plt.subplot(312)
plt.plot(t1, yabs[0 : int(N/2)]*2/N)

xd1 = [62.5, 90.82, 100.6, 119.1, 201.2]
yd1 = [0.36, 0.38, 0.79, 0.25, 0.15]

for (a, b) in zip(xd1, yd1):
	plt.plot(a, b,'r.')
	plt.text(a, b, (a, b), ha='center', va='bottom', fontsize=10)

#画包络谱图，包络谱上画出一些点，这些点需要额外计算
plt.subplot(313)
plt.plot(t1, fabs2[0 : int(N/2)]*2/N)

xd2 = [28.32]
yd2 = [0.2]

for (a, b) in zip(xd2, yd2):
	plt.plot(a, b,'r.')
	plt.text(a, b, (a, b), ha='center', va='bottom', fontsize=10)

plt.show()

• 图1是振动波形数据。
• 图2是波形数据的fft变换，我们为了显示点坐标，把fft图像做了放大处理。python最大的缺点就是画出的图不能用鼠标打点，导致我们想画的点必须提前计算好，再用程序画出来。
• 图3是对数据进行hilbert变换后再进行fft变换得到的。

从图2，我们可以算出边频的间隔大约是28.32Hz，这个频率恰好就是图3中的频率，这个频率的幅值在低频部分是峰值最大的一个。经过和设备频率表对比，我们发现这个频率和轴承对应的转轴频率非常接近。这说明轴上出现了故障隐患，所以轴转频28.32Hz对滚动体过外圈频率90.82hz产生了调制现象（和报警频率有差异是因为匹配产生了误差）。这个案例比较简单，证明了hilbert解调算法在工业信号解调中的有效性。

4、Hilbert解调应用2

本节我们分析一组轴承试验台数据，电机为恒转速，其转速 R=1496r/min，轴承的大径 D=80mm，小径 d=35mm，滚动体个数为Z=8，接触角为0度。轴承的点蚀实验中，轴承各元件上的点蚀均为单点点蚀，点蚀缺陷的大小均为直径2mm,深0.1mm的小凹坑。轴承由滚动体、内圈、外圈、保持架四类元件组成。通过这些参数，我们可以算出轴承的特征频率表。我们这里只列出了仿真故障轴承的相关频率，其他设备部件的特征频率没有列出。本案例我们将会使用第二代小波分析技术辅助hilbert解调得出故障诊断结果。

我们用6000Hz采样频率进行采样，其波形和频谱图如下图a，b所示。

从图(b)我们很难直接分析出故障频率，但是我们看到在 2000Hz （某个高频特征频率）附近有明显调制现象产生，位置恰好为小波包 1 层分解第 2 个频带附近。因此对轴承点蚀信号进行 1 层小波分解，然后对(1，2)频带进行单支重构。对重构信号进行 Hilbert 解调，再进行频谱分析， 我们得到了内圈特征频率及其多倍频（很接近），如下图(d)所示。这样我们就实现了对轴承内圈故障频率的识别。

这个案例我们还是用到了Hilbert 解调的原理，但是因为信号整体比较复杂，如果直接用Hilbert 解调会导致解调后的信号依旧很复杂，故障频率不是主要成分，所以我们用了第二代小波分析技术，对部分频带的信号进行单支重构，这样就去掉了无效信号，只留下了最重要的信号部分，使得最终解调后的信号特征明显。

5、结论

本文讲到了Hilbert解调分析，对原理设计到的公式进行了推导。之前查了一些资料对里面不理解的部分进行了仿真实验，结果还是有些出乎意料，没有完全达到预期，对于这些疑惑后面如果找到答案，我会再陆续更新。

Hilbert解调是作者近期在学习的一个信号处理方法，主要目的还是想应用到故障诊断系统中，所以学到哪里博客就写到哪里吧。小波律动公司的在线检测故障诊断系统一直在不断完善中，Hilbert解调等算法也会支持。工业软件不同于仿真实验，仿真实验我们可以根据人工干预来达到想要的实验效果，相当于事后分析。但是工业软件不能像人一样做事后分析，都是事前分析，设定好了算法，直接分析数据，要得出仿真实验的结果还是有难度的。我们的目的是希望构造一个专家系统能尽量模仿人的分析行为，最终代替人，或者半代替人得出故障诊断结论。