前言
希尔伯特变换的意义本文不提，本文的目标是举例说明到底离散傅里叶变换是如何实现的，并编写程序与matlab自带的hilbert函数结果对比，验证我们的实现过程是否正确。
离散希尔伯特变换
原始信号为，其离散希尔伯特变换的定义公式为：
说明：先让原始信号与信号做卷积，然后一起合并成一个复数信号
。
问题：看上去很简单，但这里的卷积不是一般意义上的卷积的操作！
所以：实际中得到的方法是通过借助离散傅里叶变换DFT来实现的！
因此：本文就用离散傅里叶变换来实现离散希尔伯特变换。
手动实现
设原始信号为，总长度一般
(下标n是从1开始到N)，总体实现步骤可分为3步：
- 对
做DFT，得到
；
- 对
：
- 对
，然后取其虚部的数得
；
- 最后结果：
; 是一个复数信号。
相应的matlab程序：
% 说明:
% 定义: H = h(t) * 1/(pi·t) h(t)与1/(pi·t)做卷积, 但是这不是普通意义上的卷积操作
% 实际: 离散希尔伯特变换是借助离散傅里叶变换实现的!! matlab自带的函数也是用fft的
clear ; clc ;
% x = [2 6 8 -1 12 -7 2 3 2 1 8 3 8 0 3 1];
x = [1 5 7 8 11 8 4 1]; % 还是最好用2^n个
N = length(x);
% 先对原离散信号做离散傅里叶变换
X = fft(x);
% 先是一波重新组合:
Z = [];
Z(1) = X(1);
Z(2:N/2) = 2*X(2:N/2);
Z(N/2+1:N) = 0;
% 然后对重组的结果求IDFT即可:
% 希尔伯特的复数部分: 记得去imag
z = imag(ifft(Z));
% 原信号x做hilbert后的结果: hx
hx = x + z*j
% matlab自带的函数计算结果对比:
hx_zd = hilbert(x)
结果：
hx =
1 至 6 列
1.0000 - 2.4142i 5.0000 - 4.3462i 7.0000 - 2.5355i 8.0000 - 2.7249i 11.0000 + 0.4142i 8.0000 + 4.8462i
7 至 8 列
4.0000 + 4.5355i 1.0000 + 2.2249i
hx_zd =
1 至 6 列
1.0000 - 2.4142i 5.0000 - 4.3462i 7.0000 - 2.5355i 8.0000 - 2.7249i 11.0000 + 0.4142i 8.0000 + 4.8462i
7 至 8 列
4.0000 + 4.5355i 1.0000 + 2.2249i
手动实现正确！
希尔伯特变换的意义
希尔伯特变换的结果是给原始信号提供了一个幅值、频率不变，但相位平移90°的信号
。
所以，希尔伯特变换是从"时域"到"时域"的变换！只改变了相位，所以又叫90°移相滤波器；
所以，原始信号与它的希尔伯特变换
构成正交副。
现在，我们换回到最初的记法：原始信号和它对应的希尔伯特变换信号分别用 和
表示，那么对应的"解析信号"就可以用这两个东西组成：
对于这个解析信号，我们可以得它的3瞬属性：瞬时振幅、瞬时相位、瞬时频率。
瞬时振幅：
瞬时相位：
瞬时频率：
可以看出，3瞬属性是相互关联的！
瞬时振幅、瞬时相位可以直接求，有意义；
但是瞬时频率直接根据解析信号这么按公式求，是没有物理意义的！
并且离散信号，它的瞬时频率求导只能按照"差分"来近似，即：
给出matlab的实现程序：
% 根据希尔伯特变换得到3瞬属性: 瞬时振幅、瞬时相位、瞬时频率
% 注意：这里直接得到的瞬时频率无物理意义！！实际中不这么算。
clear ; clc;
x = xlsread('shuju.xlsx');
x = x(1001:1001+1023)'; % 有效数据长必须是2^n，所以我去1024，最后10几个点是0
N = length(x);
fs = 100; % 采样频率 = 1/采样间隔
tt = (0:N-1)/fs; % 时间刻度
figure(1);
subplot(2,1,1);
plot(tt,x);
axis([min(tt) max(tt) -inf inf]);
xlabel('时间'); ylabel('振幅'); title('原始信号');
% 这里直接用自带的变换函数了hilbert
h = hilbert(x);
imagh = imag(h);
subplot(2,1,2);
plot(tt,imagh);
axis([min(tt) max(tt) -inf inf]);
title('希尔伯特谱(相位滞后90°)');
figure(2);
% 瞬时振幅/包络:
dA = abs(h);
subplot(3,1,1);
plot(tt,dA,'r'); % 包络包着原信号
hold on;
plot(tt,x); % 原信号
axis([min(tt) max(tt) -inf inf]);
xlabel('时间'); ylabel('瞬时振幅'); title('3瞬-瞬时振幅');
% 瞬时相位:
dfai = atan(imagh./x);
subplot(3,1,2);
plot(tt,dfai);
axis([min(tt) max(tt) -inf inf]);
xlabel('时间'); ylabel('瞬时相位'); title('3瞬-瞬时相位');
% 瞬时频率: 离散数据只能用"差分/梯度"来近似代替
df = gradient(dfai)./gradient(tt); % 梯度代替
% df = diff(dfai)./diff(tt); % 差分代替
subplot(3,1,3);
plot(tt(1:length(df)),df);
axis([min(tt) max(tt) -inf inf]);
xlabel('时间'); ylabel('瞬时频率'); title('3瞬-瞬时频率');
效果：
图1：原始信号与希尔伯特变换后结果
图2：3瞬属性
3瞬属性中的瞬时频率，很明显可以看出它有很多的"负频率"！这很明显是错误的。
所以，直接根据"解析信号"算瞬时频率是无意义的！
所以，真正做3瞬属性的分析，做原信号的"时频谱"分析，我们用的：
—— 希尔伯特-黄变换。
