Hilbert變換是一個(gè)很有用的變換，用它來做包絡(luò)分析更是一種有效的數(shù)據(jù)處理方法?，F(xiàn)用代碼測(cè)試其變換效果

第一個(gè)程序效果如下

% Hilbert變換測(cè)試
clc
clear all
close all

ts = 0.001;
fs = 1/ts;
N = 200;
f = 50;
k = 0:N-1;
t = k*ts;

% 信號(hào)變換
% 結(jié)論：sin信號(hào)Hilbert變換后為cos信號(hào)
y = sin(2*pi*f*t);
yh = hilbert(y);    % matlab函數(shù)得到信號(hào)是合成的復(fù)信號(hào)
yi = imag(yh);      % 虛部為書上定義的Hilbert變換

figure
subplot(211)
plot(t, y)
title('原始sin信號(hào)')
subplot(212)
plot(t, yi)
title('Hilbert變換信號(hào)')

% 檢驗(yàn)兩次Hilbert變換的結(jié)果（理論上為原信號(hào)的負(fù)值）
% 結(jié)論：兩次Hilbert變換的結(jié)果為原信號(hào)的負(fù)值
yih = hilbert(yi);
yii = imag(yih);
max(y + yii)

% 信號(hào)與其Hilbert變換的正交性
% 結(jié)論：Hilbert變換后的信號(hào)與原信號(hào)正交
sum(y.*yi)

% 譜分析
% 結(jié)論：Hilbert變換后合成的復(fù)信號(hào)的譜沒有大于奈氏頻率的頻譜，即其譜為單邊的
NFFT = 2^nextpow2(N);
f = fs*linspace(0,1,NFFT);
Y = fft(y, NFFT)/N;
YH = fft(yh, NFFT)/N;

figure
subplot(211)
plot(f,abs(Y))
title('原信號(hào)的雙邊譜')
xlabel('頻率f (Hz)')
ylabel('|Y(f)|')
subplot(212)
plot(f,abs(YH))
title('信號(hào)Hilbert變換后組成的復(fù)信號(hào)的雙邊譜')
xlabel('頻率f (Hz)')
ylabel('|YH(f)|')

第二個(gè)效果如下

第一個(gè)包絡(luò)測(cè)試

可以看到，此包絡(luò)分析得到的包絡(luò)信號(hào)頻率為20Hz，包絡(luò)信號(hào)的波形為余弦信號(hào)的絕對(duì)值信號(hào)，這是因?yàn)橛?jì)算包絡(luò)時(shí)是取絕對(duì)值得到的，從而使信號(hào)頻率加倍。解決方法是把包絡(luò)提升，遠(yuǎn)離0，如下第二個(gè)包絡(luò)。

第二個(gè)包絡(luò)測(cè)試

可以看到Hilbert包絡(luò)分析可以有效提取包絡(luò)和調(diào)制信號(hào)頻率，和檢波有一樣的效果，而且更實(shí)用。

第三個(gè)包絡(luò)測(cè)試

這是嘗試一個(gè)任意形狀的包絡(luò)，可以看到除在邊緣處有誤差外，整體效果很好。

% 包絡(luò)分析(高中心頻率的窄帶信號(hào)分析)
% 基于：兩個(gè)信號(hào)乘積的Hilbert變換取決于高頻信號(hào)的Hilbert變換
clc
clear all
close all

ts = 0.001;
fs = 1/ts;
N = 200;
k = 0:N-1;
t = k*ts;

% 原始信號(hào)
f1 = 10;
f2 = 70;
% a = cos(2*pi*f1*t);       % 包絡(luò)1
a = 2 + cos(2*pi*f1*t);     % 包絡(luò)2
% a = 1./(1+t.^2*50);       % 包絡(luò)3
m = sin(2*pi*f2*t);         % 調(diào)制信號(hào)
y = a.*m;  % 信號(hào)調(diào)制

figure
subplot(241)
plot(t, a)
title('包絡(luò)')
subplot(242)
plot(t, m)
title('調(diào)制信號(hào)')
subplot(243)
plot(t, y)
title('調(diào)制結(jié)果')

% 包絡(luò)分析
% 結(jié)論：Hilbert變換可以有效提取包絡(luò)、高頻調(diào)制信號(hào)的頻率等
yh = hilbert(y);
aabs = abs(yh);                 % 包絡(luò)的絕對(duì)值
aangle = unwrap(angle(yh));     % 包絡(luò)的相位
af = diff(aangle)/2/pi;         % 包絡(luò)的瞬時(shí)頻率，差分代替微分計(jì)算

% NFFT = 2^nextpow2(N);
NFFT = 2^nextpow2(1024*4);      % 改善柵欄效應(yīng)
f = fs*linspace(0,1,NFFT);

YH = fft(yh, NFFT)/N;           % Hilbert變換復(fù)信號(hào)的頻譜
A = fft(aabs, NFFT)/N;          % 包絡(luò)的頻譜

subplot(245)
plot(t, aabs, t, a, '.')
title('包絡(luò)的絕對(duì)值')
legend('包絡(luò)分析結(jié)果', '真實(shí)包絡(luò)')
subplot(246)
plot(t, aangle)
title('調(diào)制信號(hào)的相位')
subplot(247)
plot(t(1:end-1), af*fs)
title('調(diào)制信號(hào)的瞬時(shí)頻率')
subplot(244)
plot(f,abs(YH))
title('原始信號(hào)的Hilbert譜')
xlabel('頻率f (Hz)')
ylabel('|YH(f)|')
subplot(248)
plot(f,abs(A))
title('包絡(luò)的頻譜')
xlabel('頻率f (Hz)')
ylabel('|A(f)|')