Main Content

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

ヒルベルト変換および瞬時周波数

ヒルベルト変換で推定されるのは、単一成分の信号の瞬時周波数のみです。単一成分の信号は、単一の "リッジ" によって時間-周波数平面に表されます。一連の単一成分信号には単一の正弦波とチャープのような信号が含まれます。

1 kHz で 2 秒間サンプリングしたチャープを生成します。初期周波数が 100 Hz で、1 秒後には 200 Hz に増加するようにチャープを指定します。

fs = 1000;
t = 0:1/fs:2-1/fs;
y = chirp(t,100,1,200);

関数 pspectrum に実装された短時間フーリエ変換を使用して、チャープのスペクトログラムを推定します。信号は、各時点の単一のピーク周波数によって明確に表されます。

pspectrum(y,fs,'spectrogram')

瞬時周波数を測定するために解析信号を計算して位相を微分します。導関数をスケーリングすることで有効な推定が得られます。

z = hilbert(y);
instfrq = fs/(2*pi)*diff(unwrap(angle(z)));

clf
plot(t(2:end),instfrq)
ylim([0 fs/2])

関数 instfreq は、1 ステップで瞬時周波数を計算し、表示します。

instfreq(y,fs,'Method','hilbert')

信号が単一成分ではない場合は、メソッドが失敗します。

1023 Hz で 2 秒間サンプリングされた周波数 60 Hz と 90 Hz の 2 つの正弦波の和を生成します。スペクトログラムを計算してプロットします。各時点で 2 つの成分の存在が示されます。

fs = 1023;
t = 0:1/fs:2-1/fs;
x = sin(2*pi*60*t)+sin(2*pi*90*t);

pspectrum(x,fs,'spectrogram')
yticks([60 90])

解析信号を計算し、その位相を微分します。正弦波の周波数を囲む領域でズームインします。解析信号により、正弦波周波数の平均である瞬時周波数が予測されます。

z = hilbert(x);
instfrq = fs/(2*pi)*diff(unwrap(angle(z)));

plot(t(2:end),instfrq)
ylim([60 90])
xlabel('Time (s)')
ylabel('Frequency (Hz)')

関数 instfreq は、平均の推定も行います。

instfreq(x,fs,'Method','hilbert')

時間の関数として両方の周波数を推定するには、spectrogram を使用して、パワー スペクトル密度と tfridge を探し、2 つのリッジを追跡します。tfridge で、周波数の変更に対してペナルティを 0.1 に設定します。

[s,f,tt] = pspectrum(x,fs,'spectrogram');

numcomp = 2;
[fridge,~,lr] = tfridge(s,f,0.1,'NumRidges',numcomp);

pspectrum(x,fs,'spectrogram')
hold on
plot3(tt,fridge,abs(s(lr)),'LineWidth',4)
hold off
yticks([60 90])

参考

|

関連するトピック