Main Content

fft

説明

Y = fft(X) は高速フーリエ変換 (FFT) アルゴリズムを使用して、X離散フーリエ変換 (DFT) を計算します。

  • X がベクトルの場合、fft(X) はそのベクトルのフーリエ変換を返します。

  • X が行列の場合、fft(X) は、X の列をベクトルとして扱い、各列のフーリエ変換を返します。

  • X が多次元配列の場合、fft(X) は、サイズが 1 ではない最初の配列次元に沿った値をベクトルとして扱い、各ベクトルのフーリエ変換を返します。

Y = fft(X,n)n 点の DFT を返します。値を指定しない場合、Y のサイズは X と同じです。

  • X がベクトルであり、X の長さが n より短い場合、n の長さになるように X の末尾をゼロで埋めます。

  • X がベクトルであり、X の長さが n を超える場合、X が長さ n で切り捨てられます。

  • X が行列の場合、各列はベクトルの場合と同様に扱われます。

  • X が多次元配列の場合、サイズが 1 でない最初の配列次元がベクトルの場合と同様に扱われます。

Y = fft(X,n,dim) は、次元 dim に沿ったフーリエ変換を返します。たとえば、X が行列の場合、fft(X,n,2) は、各行の n 点のフーリエ変換を返します。

すべて折りたたむ

フーリエ変換を使用して、ノイズに埋もれた信号の周波数成分を求めます。

信号のパラメーターとして、サンプリング周波数 1 kHz、信号の持続期間 1.5 秒を指定します。

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

50 Hz、振幅 0.7 の正弦波と 120 Hz、振幅 1 の正弦波で信号を構成します。

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

平均値 0、分散 4 のホワイト ノイズで信号を乱します。

X = S + 2*randn(size(t));

ノイズを含む信号を時間領域にプロットします。信号 X(t) を見て周波数成分を特定することは困難です。

plot(1000*t(1:50),X(1:50))
title("Signal Corrupted with Zero-Mean Random Noise")
xlabel("t (milliseconds)")
ylabel("X(t)")

Figure contains an axes object. The axes object with title Signal Corrupted with Zero-Mean Random Noise contains an object of type line.

信号のフーリエ変換を計算します。

Y = fft(X);

両側スペクトル P2 を計算します。次に、P2 および偶数の信号長 L に基づいて、片側スペクトル P1 を計算します。

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

周波数領域 f を定義し、片側振幅スペクトル P1 をプロットします。ノイズを追加したため、予期したとおり、振幅は正確に 0.7 と 1 にはなりません。平均的に、信号が長くなるほど周波数がよりよく近似されます。

f = Fs*(0:(L/2))/L;
plot(f,P1) 
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Amplitude Spectrum of X(t) contains an object of type line.

ここで、ノイズを含まない元の信号のフーリエ変換を計算すると、正確な振幅 0.7 と 1.0 が得られます。

Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

plot(f,P1) 
title("Single-Sided Amplitude Spectrum of S(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Amplitude Spectrum of S(t) contains an object of type line.

ガウス パルスを時間領域から周波数領域に変換します。

信号のパラメーターとして、サンプリング周波数 44.1 kHz、信号の持続期間 1 ミリ秒を指定します。標準偏差が 0.1 ミリ秒のガウス パルスを作成します。

Fs = 44100;         % Sampling frequency
T = 1/Fs;           % Sampling period
t = -0.5:T:0.5;     % Time vector
L = length(t);      % Signal length

X = 1/(0.4*sqrt(2*pi))*(exp(-t.^2/(2*(0.1*1e-3)^2)));

パルスを時間領域にプロットします。

plot(t,X)
title("Gaussian Pulse in Time Domain")
xlabel("Time (t)")
ylabel("X(t)")
axis([-1e-3 1e-3 0 1.1]) 

Figure contains an axes object. The axes object with title Gaussian Pulse in Time Domain contains an object of type line.

関数 fft の実行時間は、変換の長さに依存します。変換の長さが小さな素因数のみからなる場合、実行時間は大きな素因数からなる場合よりもかなり速くなります。

この例では、信号長 L は 44,101 で、非常に大きな素数です。関数 fft のパフォーマンスを向上させるには、入力長として、元の信号長の次の 2 のべき乗を指定します。この入力長を使用して fft を呼び出すと、変換の長さが指定した値になるように、パルス X の末尾にゼロがパディングされます。

n = 2^nextpow2(L);

ガウス パルスを周波数領域に変換します。

Y = fft(X,n);

周波数領域を定義し、固有周波数をプロットします。

f = Fs*(0:(n/2))/n;
P = abs(Y/n).^2;

plot(f,P(1:n/2+1)) 
title("Gaussian Pulse in Frequency Domain")
xlabel("f (Hz)")
ylabel("|P(f)|^2")

Figure contains an axes object. The axes object with title Gaussian Pulse in Frequency Domain contains an object of type line.

時間領域および周波数領域で複数の余弦波を比較します。

信号のパラメーターとして、サンプリング周波数 1 kHz、信号の持続期間 1 秒を指定します。

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

周波数がスケーリングされた余弦波を各行が表す行列を作成します。その結果の X は 3 行 1000 列の行列です。1 行目の波の周波数は 50、2 行目の波の周波数は 150、3 行目の波の周波数は 300 です。

x1 = cos(2*pi*50*t);          % First row wave
x2 = cos(2*pi*150*t);         % Second row wave
x3 = cos(2*pi*300*t);         % Third row wave

X = [x1; x2; x3];

X の各行の最初から 100 個の要素を順番に 1 つの Figure にプロットし、それらの周波数を比較します。

for i = 1:3
    subplot(3,1,i)
    plot(t(1:100),X(i,1:100))
    title("Row " + num2str(i) + " in the Time Domain")
end

Figure contains 3 axes objects. Axes object 1 with title Row 1 in the Time Domain contains an object of type line. Axes object 2 with title Row 2 in the Time Domain contains an object of type line. Axes object 3 with title Row 3 in the Time Domain contains an object of type line.

X の行に沿って、つまり、各信号に対して fft を使用するように引数 dim を指定します。

dim = 2;

信号のフーリエ変換を計算します。

Y = fft(X,L,dim);

各信号の両側スペクトルと片側スペクトルを計算します。

P2 = abs(Y/L);
P1 = P2(:,1:L/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);

1 つの Figure 内で、周波数領域に各行の片側振幅スペクトルをプロットします。

for i=1:3
    subplot(3,1,i)
    plot(0:(Fs/L):(Fs/2-Fs/L),P1(i,1:L/2))
    title("Row " + num2str(i) + " in the Frequency Domain")
end

Figure contains 3 axes objects. Axes object 1 with title Row 1 in the Frequency Domain contains an object of type line. Axes object 2 with title Row 2 in the Frequency Domain contains an object of type line. Axes object 3 with title Row 3 in the Frequency Domain contains an object of type line.

周波数 15 Hz および 40 Hz の 2 つの正弦波で構成される信号を作成します。最初の正弦波は位相が -π/4 の余弦波で、2 番目の正弦波は位相が π/2 の余弦波です。この信号を 100 Hz で 1 秒間サンプリングします。

Fs = 100;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*15*t - pi/4) + cos(2*pi*40*t + pi/2);

信号のフーリエ変換を計算します。変換の振幅を周波数の関数としてプロットします。

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*Fs;
stem(f,abs(z))
title("Double-Sided Amplitude Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

Figure contains an axes object. The axes object with title Double-Sided Amplitude Spectrum of x(t) contains an object of type stem.

振幅の小さい変換の値を除去して、変換の位相を計算します。この位相を周波数の関数としてプロットします。

tol = 1e-6;
z(abs(z) < tol) = 0;

theta = angle(z);

stem(f,theta/pi)
title("Phase Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("Phase/\pi")
grid

Figure contains an axes object. The axes object with title Phase Spectrum of x(t) contains an object of type stem.

ゼロでパディングして、信号のフーリエ変換を内挿します。

信号のパラメーターとして、サンプリング周波数 80 Hz、信号の持続期間 0.8 秒を指定します。

Fs = 80;
T = 1/Fs;
L = 65;
t = (0:L-1)*T;

2 Hz の正弦波信号とさらに高い高調波の重ね合わせを作成します。この信号には 2 Hz の余弦波、4 Hz の余弦波、6 Hz の正弦波が含まれます。

X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t);

この信号を時間領域にプロットします。

plot(t,X)
title("Signal superposition in time domain")
xlabel("t (ms)")
ylabel("X(t)")

Figure contains an axes object. The axes object with title Signal superposition in time domain contains an object of type line.

信号のフーリエ変換を計算します。

Y = fft(X);

信号の片側振幅スペクトルを計算します。

f = Fs*(0:(L-1)/2)/L;
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);

周波数領域にこの片側振幅スペクトルをプロットします。信号の時間サンプリングが非常に短いため、フーリエ変換の周波数分解能には 4 Hz 付近でのピーク周波数を表示するための十分な精度がありません。

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Spectrum of Original Signal contains an object of type line.

ピーク周波数をより正確に評価するには、元の信号をゼロでパディングし、解析ウィンドウの長さを大きくします。この方法により、信号のフーリエ変換が自動的に内挿され、周波数分解能も高くなります。

新しい入力長として、元の信号長の次の 2 のべき乗を指定します。信号 X の末尾をゼロでパディングし、長さを大きくします。ゼロでパディングされた信号のフーリエ変換を計算します。

n = 2^nextpow2(L);
Y = fft(X,n);

パディングされた信号の片側振幅スペクトルを計算します。信号長 n が 65 から 128 に拡大したため、周波数分解能は Fs/n (0.625 Hz) となります。

f = Fs*(0:(n/2))/n;
P2 = abs(Y/L);
P1 = P2(1:n/2+1);
P1(2:end-1) = 2*P1(2:end-1);

パディングされた信号の片側スペクトルをプロットします。この新しいスペクトルでは、0.625 Hz の周波数分解能の範囲内で 2 Hz、4 Hz、6 Hz 付近のピーク周波数が表示されます。

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Padded Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Spectrum of Padded Signal contains an object of type line.

入力引数

すべて折りたたむ

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

X が 0 行 0 列の空の行列である場合、fft(X) は 0 行 0 列の空の行列を返します。

データ型: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical
複素数のサポート: あり

変換の長さ。[] または非負の整数スカラーとして指定します。変換の長さとして正の整数スカラーを指定すると、fft のパフォーマンスが向上することがあります。長さは通常 2 のべき乗、あるいは小さい素数 (7 以下の素因数) の積に因数分解可能な値です。n が信号の長さ未満である場合、fftn 番目の要素から後の残りの信号値を無視し、切り捨て後の結果を返します。n0 の場合、fft は空の行列を返します。

例: n = 2^nextpow2(size(X,1))

データ型: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

演算の対象の次元。正の整数のスカラーとして指定します。次元を指定しない場合、既定値はサイズが 1 より大きい最初の配列次元です。

  • fft(X,[],1)X の列に沿って演算し、各列のフーリエ変換を返します。

    fft(X,[],1) column-wise operation

  • fft(X,[],2)X の行に沿って演算し、各行のフーリエ変換を返します。

    fft(X,[],2) row-wise operation

dimndims(X) よりも大きい場合、fft(X,[],dim)X を返します。n が指定された場合、fft(X,n,dim) は次元 dim に沿って埋め込みまたは切り捨てを行うことにより、X を長さ n にします。

データ型: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

出力引数

すべて折りたたむ

周波数領域の表現。ベクトル、行列または多次元配列として返されます。

X の型が single である場合、fft はネイティブ レベルの単精度で計算し、Y の型も single になります。それ以外の場合、Ydouble 型として返されます。

Y のサイズは次のとおりです。

  • Y = fft(X) または Y = fft(X,[],dim) の場合、Y のサイズは X のサイズに等しくなります。

  • Y = fft(X,n,dim) の場合、size(Y,dim) の値は n に等しくなりますが、その他すべての次元のサイズは X と同じままです。

X が実数の場合、Y は共役対称になり、Y の特異点の数は ceil((n+1)/2) になります。

データ型: double | single

詳細

すべて折りたたむ

ベクトルの離散フーリエ変換

Y = fft(X) はフーリエ変換、X = ifft(Y) は逆フーリエ変換をそれぞれ実装します。長さ nX および Y の変換は、次式で定義されます。

Y(k)=j=1nX(j)Wn(j1)(k1)X(j)=1nk=1nY(k)Wn(j1)(k1),

ここで、

Wn=e(2πi)/n

は 1 の n 乗根の 1 つです。

ヒント

  • 関数 fft の実行時間は、変換の長さに依存します。変換の長さが小さな素因数 (7 以下) のみからなる場合は、素数である場合や大きな素因数からなる場合よりも実行時間がかなり速くなります。

  • ほとんどの n の値について、実数入力 DFT の計算時間は複素数入力 DFT の約半分になります。ただし、n が大きな素因数をもつ場合、速度の差はほとんどありません。

  • ユーティリティ関数 fftw を使用して、fft の処理速度を向上できます。この関数は、特定のサイズと次元をもつ FFT の計算に使用されるアルゴリズムの最適化を制御します。

アルゴリズム

FFT 関数 (fftfft2fftnifftifft2ifftn) は、FFTW [1][2]と呼ばれるライブラリに基づいています。

参照

[2] Frigo, M., and S. G. Johnson. “FFTW: An Adaptive Software Architecture for the FFT.” Proceedings of the International Conference on Acoustics, Speech, and Signal Processing. Vol. 3, 1998, pp. 1381-1384.

拡張機能

GPU コード生成
GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

バージョン履歴

R2006a より前に導入