fft

高速フーリエ変換

構文

Y = fft(X)

Y = fft(X,n)

Y = fft(X,n,dim)

説明

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

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

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

このプロットは、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)|")

3 つの周波数ピークの振幅を求めるには、Y の fft スペクトルを片側振幅スペクトルに変換します。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)|")

ここで、ノイズを含まない元の信号のフーリエ変換を計算すると、正確な振幅 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)|")

ガウス パルス

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

信号のパラメーターとして、サンプリング周波数 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]) 

関数 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)|")

余弦波

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

信号のパラメーターとして、サンプリング周波数 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

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

正弦波の位相

周波数 15 Hz および 40 Hz の 2 つの正弦波で構成される信号を作成します。最初の正弦波は位相が $-\pi /4$ の余弦波で、2 番目の正弦波は位相が $\pi /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

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

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

FFT の内挿

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

信号のパラメーターとして、サンプリング周波数 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)")

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

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

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

新しい入力長として、元の信号長の次の 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)|")

入力引数

X — 入力配列

ベクトル | 行列 | 多次元配列

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

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

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

n — 変換の長さ

[] (既定値) | 非負の整数スカラー

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

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

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

dim — 演算の対象の次元

正の整数スカラー

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

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

  

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

  

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

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

出力引数

Y — 周波数領域の表現

ベクトル | 行列 | 多次元配列

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

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

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) は逆フーリエ変換をそれぞれ実装します。長さ n の X および Y の変換は、次式で定義されます。

$$\begin{array}{l}Y(k)={\displaystyle \sum _{j=1}^{n}X}(j)W_n^{(j-1)\text{}(k-1)}\\ X(j)=\frac{1}{n}{\displaystyle \sum _{k=1}^{n}Y}(k)W_n{}^{-(j-1)\text{}(k-1)},\end{array}$$

ここで、

$$W_n=e^{(-2\pi i)/n}$$

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

ヒント

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

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

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

アルゴリズム

FFT 関数 (fft、fft2、fftn、ifft、ifft2、ifftn) は、FFTW というライブラリに基づいています[1] [2]。

参照

[1] FFTW (https://www.fftw.org)

[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.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• 入力引数 X では、次のようになります。

  • 次元を指定しない場合、コード ジェネレーターは、可変サイズであるかサイズが 1 と等しくない入力配列の最初の次元に沿って演算を行います。この次元がコード生成時に可変サイズで実行時に 1 になる場合、ランタイム エラーが発生することがあります。このエラーを回避するには、次元を指定します。

• MEX 出力の場合、MATLAB^® Coder™ は MATLAB が FFT アルゴリズムに使用するライブラリを使用します。スタンドアロン C/C++ コードの場合、コード ジェネレーターは既定で、FFT ライブラリの呼び出しを生成する代わりに FFT アルゴリズムのコードを生成します。特定のインストールされた FFTW ライブラリの呼び出しを生成するには、FFT ライブラリ コールバック クラスを指定します。FFT ライブラリ コールバック クラスの詳細については、coder.fftw.StandaloneFFTW3Interface (MATLAB Coder) を参照してください。

• MATLAB Function ブロックのシミュレーションの場合、シミュレーション ソフトウェアは MATLAB が FFT アルゴリズムに使用するライブラリを使用します。C/C++ コード生成の場合、コード ジェネレーターは既定で、FFT ライブラリの呼び出しを生成する代わりに FFT アルゴリズム用のコードを生成します。特定のインストールされた FFTW ライブラリの呼び出しを生成するには、FFT ライブラリ コールバック クラスを指定します。FFT ライブラリ コールバック クラスの詳細については、coder.fftw.StandaloneFFTW3Interface (MATLAB Coder) を参照してください。

• コード置換ライブラリ (CRL) を使用して、Neon 拡張を含む ARM^® Cortex^®-A Processors で実行される最適化されたコードを生成できます。最適化されたコードを生成するには、Embedded Coder^® Support Package for ARM Cortex-A Processors (Embedded Coder) をインストールしなければなりません。ARM Cortex-A で生成されたコードは Ne10 ライブラリを使用します。詳細については、Ne10 Conditions for MATLAB Functions to Support ARM Cortex-A Processors (Embedded Coder) を参照してください。

• コード置換ライブラリ (CRL) を使用して、ARM Cortex-M Processors で実行される最適化されたコードを生成できます。最適化されたコードを生成するには、Embedded Coder Support Package for ARM Cortex-M Processors (Embedded Coder) をインストールしなければなりません。ARM Cortex-M で生成されたコードでは、CMSIS ライブラリを使用します。詳細については、CMSIS Conditions for MATLAB Functions to Support ARM Cortex-M Processors (Embedded Coder) を参照してください。

GPU コード生成

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

使用上の注意および制限については、「C/C++ コード生成」セクションを参照してください。GPU コード生成にも同様の、使用上の注意および制限が適用されます。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

fft 関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

fft 関数は GPU 配列入力をサポートしますが、次の使用上の注意および制限があります。

• 虚数部がすべてゼロであっても、出力 Y は常に複素数です。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

分散配列

Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。

fft 関数は分散配列をサポートしますが、次の使用上の注意および制限があります。

• 分散配列では、並列 FFT アルゴリズムを使用する代わりに fft は単一のワーカーでベクトルを収集し、素数長の FFT を実行します。大きい素数長ベクトルの FFT では、メモリ不足のエラーが発生する場合があります。

詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

関数

• fft2 | fftn | fftw | fftshift | ifft | interpft

トピック

• フーリエ変換

授業用リソース

• Fourier Analysis (MathWorks Teaching Resources)
• Convolution in Digital Signal Processing (MathWorks Teaching Resources)