|
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)|')
|
