Main Content

fft

説明

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

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

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

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

Y = fft(X,n)n 点の DFT を返します。

  • 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

振幅 0.8 の DC オフセット、振幅 0.7 の 50 Hz 正弦波、振幅 1 の 120 Hz 正弦波で信号を構成します。

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

平均値 0、分散 4 のランダム ノイズで信号を乱します。

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

ノイズを含む信号を時間領域にプロットします。周波数成分は、プロットでは視覚的に明らかではありません。

plot(1000*t,X)
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, xlabel t (milliseconds), ylabel X(t) contains an object of type line.

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

Y = fft(X);

フーリエ変換は複素数を含むため、fft スペクトルの複素数の大きさをプロットします。

plot(Fs/L*(0:L-1),abs(Y),"LineWidth",3)
title("Complex Magnitude of fft Spectrum")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Figure contains an axes object. The axes object with title Complex Magnitude of fft Spectrum, xlabel f (Hz), ylabel |fft(X)| contains an object of type line.

このプロットは、DC オフセットの 0 Hz のピークを含む 5 つの周波数ピークを示しています。この例では、信号は 0 Hz、50 Hz、120 Hz の 3 つの周波数ピークをもつと予想されます。このプロットの後半は、0 Hz のピークを含まない形で、前半と左右対称になっています。この理由は、時間領域信号の離散フーリエ変換は周期的な性質をもっており、スペクトルの前半が正の周波数、後半が負の周波数となり、最初の要素はゼロ周波数用に予約されているからです。

実数の信号の場合、fft スペクトルは両側スペクトルであり、正の周波数のスペクトルは負の周波数のスペクトルの複素共役です。正と負の周波数における fft スペクトルを示すには、fftshift を使用できます。L の偶数長の場合、周波数領域は負のナイキスト周波数 -Fs/2 から Fs/2-Fs/L までで、間隔または周波数分解能は Fs/L です。

plot(Fs/L*(-L/2:L/2-1),abs(fftshift(Y)),"LineWidth",3)
title("fft Spectrum in the Positive and Negative Frequencies")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Figure contains an axes object. The axes object with title fft Spectrum in the Positive and Negative Frequencies, xlabel f (Hz), ylabel |fft(X)| contains an object of type line.

3 つの周波数ピークの振幅を求めるには、Yfft スペクトルを片側振幅スペクトルに変換します。fft 関数には、元の信号と変換後の信号の間のスケーリング係数 L が含まれているため、L で除算して Y を再スケーリングします。fft スペクトルの複素数の大きさを取ります。両側振幅スペクトル P2 は、正の周波数のスペクトルが負の周波数のスペクトルの複素共役であり、時間領域信号のピーク振幅が半分になります。片側スペクトルに変換するには、両側スペクトル P2 の前半を取ります。正の周波数のスペクトルに 2 を乗算します。P1(1)P1(end) に 2 を乗算する必要はありません。これらの振幅はそれぞれゼロ周波数とナイキスト周波数に対応し、負の周波数では複素共役のペアをもたないからです。

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

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

f = Fs/L*(0:(L/2));
plot(f,P1,"LineWidth",3) 
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), xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

ここで、ノイズを含まない元の信号のフーリエ変換を計算すると、正確な振幅 0.8、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,"LineWidth",3) 
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), xlabel f (Hz), ylabel |P1(f)| 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, xlabel Time (t), ylabel X(t) 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/sqrt(n)).^2;

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

Figure contains an axes object. The axes object with title Gaussian Pulse in Frequency Domain, xlabel f (Hz), ylabel |P(f)| 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), xlabel Frequency (Hz), ylabel |y| 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), xlabel Frequency (Hz), ylabel Phase/ pi 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, xlabel t (ms), ylabel X(t) 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, xlabel f (Hz), ylabel |P1(f)| 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, xlabel f (Hz), ylabel |P1(f)| 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 より前に導入