希爾伯特變換：給信號做個“全身體檢”的數學黑科技

一、 揭開面紗：什么是希爾伯特變換？

在現實世界中，傳感器采集到的信號（如麥克風錄下的聲音、地震儀記錄的震動）都是實信號。實信號就像是一個人的平面照片，雖然能看出樣貌，但缺乏立體感。數學家大衛·希爾伯特（David Hilbert）提出了一種天才的構想：如果我們給這個實信號強行加上一個“虛幻的影子”，讓它變成一個復數空間里的解析信號，會發生什么？

具體而言，希爾伯特變換的作用，就是把原信號中的所有頻率成分，統一進行90度的相移。原信號作為“實部”，變換后的信號作為“虛部”，兩者結合，就構成了一個完美旋轉的解析信號向量。通過這波操作，我們獲得了三個極為關鍵的“上帝視角”數據：瞬時幅值（包絡）、瞬時相位、瞬時頻率。

希爾伯特變換的數學定義可以從時域卷積、頻域變換以及解析信號三個角度來理解，它們本質上是等價的。

Matlab手動實現希爾伯特變換

% ==============================
% 希爾伯特變換 - 手動公式實現
% ==============================
clear;
clc;
% 1. 測試信號
fs = 48000;
f=1000;
t = 0:1/fs:1-1/fs;
x = sin(2*pi*f*t);
N = length(x);
% 2. FFT
X = fft(x);
% 3. 構造頻域希爾伯特濾波器：-j * sgn(f)
f = (0:N-1)*(fs/N);
H = -1j * sign(f - fs/2);
% 頻域核
% 4. 頻域相乘 → IFFT
X_hilbert = X .* H;
x_hilbert = ifft(X_hilbert);
% 5. 取實部（消除數值虛部）
x_hilbert = real(x_hilbert);
% 畫圖
figure(1), plot(t(1:200), x(1:200), 'b', t(1:200), x_hilbert(1:200), 'r');
legend('原始','希爾伯特變換');
grid on;

四、 工業聽診器：機械故障預測

在工業監測領域，當風力發電機的軸承出現早期微小裂紋時，它會產生一種高頻的沖擊信號，但這種信號通常被巨大的機器轟鳴聲掩蓋。工程師通過希爾伯特變換，對振動信號進行包絡分析（即提取瞬時幅值）。這就像是給信號剝去了一層層厚厚的洋蔥皮，直接暴露出低頻的故障特征頻率，實現設備的預測性維護。

五、音頻魔術師：聲音領域的奇妙應用

在我們日常聆聽音樂和語音交互中，希爾伯特變換同樣扮演著"音頻魔術師"的角色，悄無聲息地提升著我們的聽覺體驗。

1. 包絡提取與動態壓縮
通過希爾伯特變換構造解析信號，我們可以輕松提取出音頻信號的瞬時幅度包絡，這個包絡信息反映了聲音能量的起伏變化，被廣泛用于自動增益控制和動態范圍壓縮，讓輕柔的音符更清晰，爆炸聲也不會過于刺耳，優化整體聽感。

2. 泛音分析與音色處理
音樂之所以動人，在于豐富的泛音層次。希爾伯特變換能夠幫助我們分離出音頻信號中的泛音成分，提取各次諧波的瞬時頻率和相位，這對于樂器識別、音色分析以及電子音樂合成都至關重要。通過調整泛音結構，我們甚至可以讓一把小提琴聽起來像是另一件樂器。

3. 聲反饋抑制（防嘯叫）
在會議室、演唱會現場，經常會出現令人頭疼的嘯叫問題。這是因為揚聲器發出的聲音被麥克風重新采集，形成了正反饋。希爾伯特變換通過對信號進行90°相移，改變了反饋路徑的相位關系，能夠有效破壞嘯叫產生的條件，讓擴聲系統更加穩定。

4. 音頻調制與解調
在一些特殊的音頻傳輸場景中，希爾伯特變換用于完成幅度調制信號的解調，精確恢復出原始的音頻信息，保證了語音和音樂傳輸的保真度。

從HIFI耳機的音質優化，到手機語音助手的語音增強，再到現場演出的嘯叫抑制，希爾伯特變換以其數學之美，讓我們聽到的聲音更加清晰動聽。

六、希爾伯特變換Matlab代碼應用

在前面我們已經看到了希爾伯特變換的手動實現，現在讓我們通過幾個完整的實例，來感受一下Matlab中希爾伯特變換的強大功能。

1. 基礎實例：信號移頻

% ==============================
% 希爾伯特變換 + 1000Hz 移頻（正確版）
% ==============================
clear; clc;
%% 1. 生成信號
fs = 48000;
f1 = 1000; % 原始 1000Hz
t = 0:1/fs:1-1/fs;
x = sin(2*pi*f1*t);
N = length(x);
%% 2. 手動實現希爾伯特變換（正確版）
X = fft(x);
H = zeros(1, N);
H(1) = 0; % 直流
H(2:N/2) = -1j; % 正頻率移相 -90°
H(N/2+1) = 0; % 奈奎斯特頻率
H(N/2+2:end) = 1j; % 負頻率移相 +90°
x_hilbert = ifft(X .* H);
x_hilbert = real(x_hilbert); % 希爾伯特結果
%% 3. 構造解析信號（關鍵！）
analytic = x + 1j * x_hilbert;
%% 4. 移頻 1000Hz
shift_hz = 1000;
shift_factor = exp(1j * 2*pi*shift_hz*t);
%shift_factor = cos(2*pi*shift_hz*t) + 1j*sin(2*pi*shift_hz*t);
y_analytic = analytic .* shift_factor; % 正確移頻
y = real(y_analytic); % 輸出實信號
%% 5. 畫圖看結果（前2000點看波形）
figure;
subplot(2,1,1); plot(t(1:2000),x(1:2000)); title('原始 1000Hz');
subplot(2,1,2); plot(t(1:2000),y(1:2000)); title('移頻后 1100Hz');

2. 經典應用：希爾伯特變換提取瞬時幅度、瞬時相位、瞬時頻率

% ==============================
% 希爾伯特變換 → 提取 瞬時幅度、相位、頻率
% ==============================
clear; clc;
%% 1. 生成測試信號（1000Hz 正弦波）
fs = 48000; % 采樣率
f0 = 1000; % 信號頻率
t = 0:1/fs:1-1/fs; % 時間軸
x = sin(2*pi*f0*t); % 原始信號
%% 2. 希爾伯特變換 → 解析信號
% z = hilbert(x); % MATLAB 內置（最穩定）
%% 2. 手動實現希爾伯特變換（正確版）
N = length(x);
X = fft(x);
H = zeros(1, N);
H(1) = 0; % 直流
H(2:N/2) = -1j; % 正頻率移相 -90°
H(N/2+1) = 0; % 奈奎斯特頻率
H(N/2+2:end) = 1j; % 負頻率移相 +90°
x_hilbert = ifft(X .* H);
x_hilbert_real = real(x_hilbert); % 希爾伯特結果
x_hilbert_imag = imag(x_hilbert); % 希爾伯特結果
x_hilbert = real(x_hilbert); % 希爾伯特結果
z = x + 1j * x_hilbert;
%% 3. 提取 3 個核心參數
amp = abs(z); % 瞬時幅度
phi = angle(z); % 瞬時相位（-π ~ π）
phi = unwrap(phi); % 相位解卷（關鍵！否則會跳變）
freq = diff(phi)/(2*pi)*fs; % 瞬時頻率（差分求導）
% 讓頻率和原信號長度一致（補最后一個點）
freq = [freq, freq(end)];
%% 4. 畫圖展示
figure('Position',[100,100,800,600])
subplot(3,1,1);
plot(t, amp);
title('瞬時幅度');
ylim([0,1.5]); grid on;
subplot(3,1,2);
plot(t, phi);
title('瞬時相位（已解卷）');
grid on;
subplot(3,1,3);
plot(t, freq);
title('瞬時頻率 (Hz)');
ylim([900,1100]); grid on;

3. 高階應用：變聲

% ==============================
% 希爾伯特變換 變調（改變音調）
% 原理：修改瞬時相位，不改變時長
% ==============================
clear; clc;
%% 1. 讀取音頻（必須是 wav 文件）
[input, fs] = audioread('input.wav');
input = input(:,1); % 轉單聲道
%% 2. 變調參數（自己調）
pitch_ratio = 1.3; % >1 聲音變尖（蘿莉）
% pitch_ratio = 0.8; % <1 聲音變粗（大叔）
%% 3. 希爾伯特 → 解析信號
z = hilbert(input); % 解析信號
%% 4. 提取幅度 + 相位（核心）
amp = abs(z); % 瞬時幅度（不變）
phase = unwrap(angle(z)); % 瞬時相位（解卷，關鍵！）
%% 5. 變調核心：縮放相位
new_phase = phase * pitch_ratio;
%% 6. 重構信號
output = amp .* cos(new_phase);
%% 7. 歸一化防止爆音
output = output / max(abs(output));
%% 8. 保存 & 播放
audiowrite('voice_hilbert_pitch.wav', output, fs);
sound(output, fs);
disp('變調完成！');