Main Content

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

goertzel

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

説明

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

dft = goertzel(data,findx) では、findx で指定された周波数インデックスの DFT が返されます。

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

すべて折りたたむ

電話のキーパッドの「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 アルゴリズムを使用してトーンの離散フーリエ変換 (DFT) を計算します。0 ~ 9 の番号の生成に使用される周波数に対応するインデックスを選択します。

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

DFT 振幅をプロットします。

stem(f,abs(dft_data))

ax = gca;
ax.XTick = f;
xlabel('Frequency (Hz)')
ylabel('DFT Magnitude')

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

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

dB 単位表記で離散フーリエ変換の振幅をプロットします。

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

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel DFT 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 次元配列として返されます。

アルゴリズム

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

hk(n)=ej2πkej2πkn/Nu(n)ej2πkWNknu(n),

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

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

ここで、

yk(n)=m=0Nx(m)hk(nm)

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

Hk(z)=(1WNkz1)ej2πk12cos(2πkN)z1+z2,

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

goertzel の出力と Goertzel アルゴリズムの直接実装による結果を比較します。入力信号には、50 Hz で 10 秒間サンプリングされ、ホワイト ガウス ノイズに組み込まれたチャープを使用します。測定中、チャープの周波数は 15 Hz から 20 Hz に線形で増加します。fs/N の整数倍ではない周波数で離散フーリエ変換を計算します。goertzel を呼び出す場合、MATLAB® のベクトルは 0 から N – 1 ではなく 1 から N で実行されることに注意してください。高精度の結果が得られます。

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

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

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

dft = goertzel(xn,k+1);

df = abs(Xk-dft)
df =
   4.3634e-12

代替方法

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

  • fft: いくつかの周波数のみで DFT が必要な場合は、Goertzel アルゴリズムよりも非効率的です。周波数 log2N (N は入力信号の長さ) よりも大きい周波数で変換を評価する必要がある場合、fftgoertzel よりも効率的です。

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

参照

[1] Burrus, C. Sidney, and Thomas W. Parks. DFT/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 より前に導入

参考

|