このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

goertzel

2 次 Goertzel アルゴリズムを使用した離散時間フーリエ変換

ページ内をすべて折りたたむ

構文

dft = goertzel(data)
dft = goertzel(data,findx)
dft = goertzel(data,findx,dim)

説明

dft = goertzel(data) は、2 次 Goertzel アルゴリズムを使用して、入力配列 data の離散時間フーリエ変換 (DTFT) を返します。data が複数次元配列である場合、goertzel は、サイズが 1 より大きい最初の配列次元に沿って動作します。詳細については、アルゴリズムを参照してください。

例

dft = goertzel(data,findx) は、findx で指定された周波数インデックスの DTFT を整数または非整数のベクトルとして返します。

例

dft = goertzel(data,findx,dim) では、DTFT が次元 dim に沿って計算されます。次元を入力し、findx の既定値を使用するには、2 番目の引数を空 [] として指定します。

例

例

すべて折りたたむ

ライブ スクリプトを開く

goertzel の出力と Goertzel アルゴリズムの直接型 II の実装による結果を比較します。

50 Hz で 10 秒間サンプリングされ、ホワイト ガウス ノイズに組み込まれたチャープ信号 xn を定義します。測定中、チャープの周波数は 15 Hz から 20 Hz に線形で増加します。

fs = 50;
t = 0:1/fs:10-1/fs;
N = length(t);
xn = chirp(t,15,t(end),20) + randn(1,N)/100;

離散フーリエ変換 (DTFT) の計算で使用する、fs/N の整数倍ではない周波数を定義します。 

f0 = 17.36;
k = N*f0/fs;

2 次 Goertzel アルゴリズムの直接型 II の実装を使用して、xn の DTFT を計算します。

ykn = filter([1 -exp(-2j*pi*k/N)],[1 -2*cos(2*pi*k/N) 1],[xn 0]);
Xk = exp(-2j*pi*k)*ykn(end)
Xk = 
-30.5384 -14.3396i

goertzel 関数を使用して xn の DTFT を計算します。goertzel を呼び出す場合、MATLAB® のベクトルの範囲が 0 から N - 1 ではなく 1 から N であることに注意してください。

dft = goertzel(xn,k+1)
dft = 
-30.5384 -14.3396i

出力を比較します。高精度の結果が得られます。

df = abs(Xk-dft)
df = 
4.3634e-12
ライブ スクリプトを開く

電話のキーパッドの「1」ボタンを押して生成されるトーンの周波数を推定します。

番号「1」を押すと、周波数 697 と 1209 Hz をもつトーンが生成されます。サンプル レート 8 kHz のトーンのサンプルを 205 個生成します。

Fs = 8000;
N = 205;
lo = sin(2*pi*697*(0:N-1)/Fs);
hi = sin(2*pi*1209*(0:N-1)/Fs);
data = lo + hi;

Goertzel アルゴリズムを使用してトーンの離散時間フーリエ変換 (DTFT) を計算します。0 ～ 9 の番号の生成に使用される周波数に対応するインデックスを選択します。

f = [697 770 852 941 1209 1336 1477];
freq_indices = round(f/Fs*N) + 1;   
dtft_data = goertzel(data,freq_indices);

DTFT の振幅をプロットします。

stem(f,abs(dtft_data))

ax = gca;
ax.XTick = f;
xlabel("Frequency (Hz)")
ylabel("DTFT Magnitude")

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel DTFT Magnitude contains an object of type stem.

ライブ スクリプトを開く

goertzel 関数および 0 Hz を中心とする周波数を使用して、複素数値チャープ信号の DTFT を計算します。

1 kHz で 10 秒間サンプリングされた、線形スイープ複素数値チャープと 2 次スイープ実数値チャープで構成される信号を生成します。線形スイープ チャープ信号の瞬時周波数は、最初は -100 Hz で、最後は 200 Hz になります。2 次スイープ チャープ信号の瞬時周波数は、最初は 350 Hz で、最後は 400 Hz になります。初期位相はゼロです。

Fs = 1e3;
t = single(0:1/Fs:10);

y = chirp(t,-100,t(end),200,"linear",0,"complex") + ...
    0.1*chirp(t,350,t(end),400,"quadratic");

Goertzel アルゴリズムを使用して信号の DTFT を計算します。0 Hz を中心とした周波数軸を指定します。DTFT の振幅をデシベル単位でプロットします。

F = linspace(-Fs/2,Fs/2,1001);

N = length(t);
k =  N*F/Fs;
dft = goertzel(y,k+1);

plot(F,db(dft))
xlabel("Frequency (Hz)")
ylabel("Magnitude (dB)")
grid on

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Magnitude (dB) contains an object of type line.

ライブ スクリプトを開く

1.24 kHz、1.26 kHz および 10 kHz の周波数成分をもつ、ノイズを含む余弦を生成します。サンプル レートを 32 kHz に指定します。再現可能な結果が必要な場合は、乱数発生器をリセットします。

rng default

Fs = 32e3;
t = 0:1/Fs:2.96;
x = cos(2*pi*t*10e3) + cos(2*pi*t*1.24e3) + cos(2*pi*t*1.26e3) ...
    + randn(size(t));

周波数ベクトルを生成します。Goertzel アルゴリズムを使用して DFT を計算します。周波数の範囲を 1.2 ～ 1.3 kHz に制限します。

N = (length(x)+1)/2;
f = (Fs/2)/N*(0:N-1);
indxs = find(f>1.2e3 & f<1.3e3);
X = goertzel(x,indxs);

平均二乗スペクトルをデシベル単位でプロットします。

plot(f(indxs)/1e3,mag2db(abs(X)/length(X)))

title('Mean Squared Spectrum')
xlabel('Frequency (kHz)')
ylabel('Power (dB)')
grid

Figure contains an axes object. The axes object with title Mean Squared Spectrum, xlabel Frequency (kHz), ylabel Power (dB) contains an object of type line.

ライブ スクリプトを開く

3.2 kHz で 10 秒間サンプリングされ、ホワイト ガウス ノイズに組み込まれた 2 チャネルの信号を生成します。信号の最初のチャネルは 124 Hz の正弦波です。2 番目のチャネルは周波数が 126 Hz の複素指数です。時間軸が 3 番目の次元に沿うように信号を 3 次元配列に形状変更します。

fs = 3.2e3;
t = 0:1/fs:10-1/fs;

x = [cos(2*pi*t*124);exp(2j*pi*t*126)] + randn(2,length(t))/100;
x = permute(x,[3 1 2]);

size(x)
ans = 1×3

           1           2       32000

Goertzel アルゴリズムを使用して、信号の離散時間フーリエ変換を計算します。周波数範囲を 120 ～ 130 Hz に制限します。

N = (length(x)+1)/2;
f = (fs/2)/N*(0:N-1);
indxs = find(f>=120 & f<=130);

X = goertzel(x,indxs,3);

離散時間フーリエ変換の振幅をデシベル単位でプロットします。

plot(f(indxs),mag2db(abs(squeeze(X))))
xlabel('Frequency (Hz)')
ylabel('DTFT Magnitude (dB)')
grid

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel DTFT Magnitude (dB) contains 2 objects of type line.

入力引数

すべて折りたたむ

入力配列。ベクトル、行列、または N 次元配列として指定します。

例: sin(2*pi*(0:255)/4) は、正弦波を行ベクトルとして指定します。

例: sin(2*pi*[0.1;0.3]*(0:39))' は、2 チャネルの正弦波を指定します。

データ型: single | double
複素数のサポート: あり

周波数インデックス。ベクトルとして指定します。このインデックスは、fs/N の整数倍または非整数倍に対応できます。ここで、fs はサンプル レート、N は信号長です。

データ型: single | double

動作する対象の次元。正の整数スカラーとして指定します。

データ型: single | double

出力引数

すべて折りたたむ

離散フーリエ変換。ベクトル、行列、または N 次元配列として返されます。

ヒント

  • いくつかの周波数でのみ DTFT が必要な場合は、Goertzel アルゴリズムを使用するのが最も効率の良い方法です。

  • 次の方法で DTFT を計算することもできます。

    • fft — log2N 個よりも多い周波数における変換結果を求める必要がある場合、fft は goertzel よりも効率的です。ここで、N は入力信号の長さです。

    • czt — czt は、円または螺旋の等高線上の入力信号のチャープ Z 変換を計算し、特殊ケースとして DTFT を含みます。

アルゴリズム

Goertzel アルゴリズムは、次のインパルス応答をもつ N 点の入力 x(n), n = 0, 1, …, N – 1 の畳み込みとして離散時間フーリエ変換 X(k) を実装します。

hk(n)=e−j2πk ej2πkn/N u(n)≡e−j2πk WN−kn u(n),

ここで、単位ステップのシーケンス u(n) は、n ≥ 0 の場合は 1 であり、それ以外の場合は 0 となります。k は整数である必要はありません。周波数 f = kfs/N のとき (fs はサンプル レート)、変換は次の値をもちます。

X(k)=yk(n)|n=N,

ここで、

yk(n)=∑m=0Nx(m) hk(n−m)

かつ x(N) = 0 です。インパルス応答の Z 変換は次のようになります。

Hk(z)=(1−WNkz−1)e−j2πk1−2cos(2πkN) z−1+z−2

これを使用して直接型 II を次のように実装します。

参照

[1] Burrus, C. Sidney, and Thomas W. Parks. DTFT/FFT and Convolution Algorithms: Theory and Implementation. New York: John Wiley & Sons, 1985.

[2] Proakis, John G., and Dimitris G. Manolakis. Digital Signal Processing: Principles, Algorithms, and Applications. 3rd Edition. Upper Saddle River, NJ: Prentice Hall, 1996.

[3] Sysel, Petr, and Pavel Rajmic. “Goertzel Algorithm Generalized to Non-Integer Multiples of Fundamental Frequency.” EURASIP Journal on Advances in Signal Processing. Vol. 2012, Number 1, December 2012, pp. 56-1–56-8. https://doi.org/10.1186/1687-6180-2012-56.

拡張機能

すべて展開する

バージョン履歴

R2006a より前に導入

参考

|