ドキュメンテーション

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

hilbert

ヒルベルト変換を使用した離散時間解析信号

構文

x = hilbert(xr)
x = hilbert(xr,n)

説明

x = hilbert(xr) 実数データ シーケンスから複素螺旋列 ("解析信号" とも呼ばれる) を返します。解析信号 x = xr + i*xi は、実数部 xr と虚数部 xi から成り、実数部は元のデータで、虚数部はそのヒルベルト変換です。虚数部は、元の実数シーケンスを 90°位相シフトしたものです。したがって、正弦値は余弦値に変換され、余弦値は正弦値に変換されます。ヒルベルト変換された級数は、元のシーケンスと同じ振幅と周波数成分をもち、元の位相に対応した位相情報を含みます。

xr が行列の場合、x = hilbert(xr) では、行列の列方向に動作して各列に対応する解析信号を検出します。

x = hilbert(xr,n) では、n 点の FFT を使用してヒルベルト変換が計算されます。入力データ xr は、長さ n になるように、ゼロが付加されるか切り取られます。

ヒルベルト変換は、振幅や周波数などの時系列の瞬間的な特性を計算するのに役立ちます。瞬間的な振幅は、複素ヒルベルト変換の振幅であり、瞬間的な周波数は、瞬間的な位相角の時間変化率です。純粋な正弦波の場合、瞬間的な振幅と瞬間的な周波数は一定値です。ただし、瞬間的な位相は、局所的な位相角が単一の周期で線形に変化する様子を反映して、ノコギリ形です。混合正弦波の場合、属性は多くても 2~3 点の短期的すなわち局所的な平均です。例については、「ヒルベルト変換および瞬間周波数」を参照してください。

参考文献 [1] には、時系列のスペクトル密度の対数のヒルベルト変換を実行して、最小位相を復元する Kolmogorov の方法が記述されています。このツールボックスの関数 rceps はこの復元を行います。

離散時間解析信号 x では、fft(x) の後半部はゼロになり、fft(x) の最初の成分 (DC) と中央の要素 (ナイキスト) には虚数部はありません。

すべて折りたたむ

シーケンスを定義し、hilbert を使用してその解析信号を計算します。

xr = [1 2 3 4];
x = hilbert(xr)
x =

   1.0000 + 1.0000i   2.0000 - 1.0000i   3.0000 - 1.0000i   4.0000 + 1.0000i

x の虚数部は、xr のヒルベルト変換であり、実数部は xr 自体です。

imx = imag(x)
rex = real(x)
imx =

     1    -1    -1     1


rex =

     1     2     3     4

x の DFT の最後の半分はゼロです。(この例では、変換の最後の半分がちょうどその最後の要素です)。fft(x) の DC 要素とナイキスト要素は純粋な実数です。

dft = fft(x)
dft =

  10.0000 + 0.0000i  -4.0000 + 4.0000i  -2.0000 + 0.0000i   0.0000 + 0.0000i

関数 hilbert は、有限データ ブロックに正確に一致する解析信号を検出します。また、FIR ヒルベルト変換器フィルターを使用して解析信号を生成し、虚数部に対する近似を計算することもできます。

周波数が 203 Hz、721 Hz および 1001 Hz である 3 つの正弦波で構成されるシーケンスを生成します。シーケンスは 10 kHz で約 1 秒間サンプリングされています。関数 hilbert を使用して解析信号を計算します。これを 0.01 秒と 0.03 秒の間でプロットします。

fs = 1e4;
t = 0:1/fs:1;

x = 2.5+cos(2*pi*203*t)+sin(2*pi*721*t)+cos(2*pi*1001*t);

y = hilbert(x);

plot(t,real(y),t,imag(y))
xlim([0.01 0.03])
legend('real','imaginary')
title('hilbert Function')

元のシーケンスおよび解析信号のパワー スペクトル密度のウェルチ推定を計算します。ハミング ウィンドウを適用した長さ 256 のオーバーラップのないセクションにシーケンスを分割します。負の周波数では解析信号にパワーがないことを確認します。

pwelch([x;y].',256,0,[],fs,'centered')
legend('Original','hilbert')

関数 designfilt を使用して 60 次のヒルベルト変換器 FIR フィルターを設計します。400 Hz の遷移幅を指定します。フィルターの周波数応答を可視化します。正弦波シーケンスをフィルター処理し、解析信号の虚数部を近似します。

fo = 60;

d = designfilt('hilbertfir','FilterOrder',fo, ...
       'TransitionWidth',400,'SampleRate',fs);

freqz(d,1024,fs)

hb = filter(d,x);

フィルターの群遅延 grd はフィルターの次数の 1/2 と等価です。この遅延を補正します。虚数部の最初の grd サンプルと実数部の最後の grd サンプルおよび時間ベクトルを削除します。0.01 秒と 0.03 秒の間の結果をプロットします。

grd = fo/2;

y2 = x(1:end-grd) + 1j*hb(grd+1:end);
t2 = t(1:end-grd);

plot(t2,real(y2),t2,imag(y2))
xlim([0.01 0.03])
legend('real','imaginary')
title('FIR Filter')

近似の解析信号の PSD を推定し、hilbert の結果と比較します。

pwelch([y;[y2 zeros(1,grd)]].',256,0,[],fs,'centered')
legend('hilbert','FIR Filter')

詳細

すべて折りたたむ

アルゴリズム

シーケンス x の解析信号は "片側フーリエ変換" です。つまり、変換は負の周波数に対してはゼロになります。解析信号を近似するため、hilbert では、入力シーケンスの FFT が計算され、負の周波数に対応する FFT 係数がゼロに置き換えられて、その結果の逆 FFT が計算されます。

詳細に説明すると、hilbert では以下の 4 ステップのアルゴリズムが使用されます。

  1. 入力シーケンスの FFT が計算され、その結果がベクトル x に格納されます。

  2. 要素 h(i) が以下の値をもつ、ベクトル h が作成されます。

    • i = 1, (n/2)+1 の場合は 1

    • i = 2, 3, ... , (n/2) の場合、2

    • i = (n/2)+2, ... , n の場合 0

  3. xh の要素単位の積が計算されます。

  4. ステップ 3 で得られたシーケンスの逆 FFT が計算され、結果から最初の n 個の要素が返されます。

このアルゴリズムは、まず [2] で紹介されました。この手法は、入力信号 x が有限データ ブロックであることを前提としています。この前提により、x のスペクトル冗長性をこの関数で厳密に排除することが可能になります。FIR フィルター処理に基づいた手法は解析信号の近似しかできませんが、データに連続して作用するという利点があります。この他の FIR フィルターを使用して計算されるヒルベルト変換の例については、「単側波帯振幅変調」を参照してください。

参考文献

[1] Claerbout, Jon F. Fundamentals of Geophysical Data Processing with Applications to Petroleum Prospecting. Oxford, UK: Blackwell, 1985, pp. 59–62.

[2] Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” IEEE Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.

[3] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

参考

| |

R2006a より前に導入

この情報は役に立ちましたか?