fft
高速フーリエ変換
説明
例
ノイズを含む信号
フーリエ変換を使用して、ノイズに埋もれた信号の周波数成分を求めて、ピーク周波数の振幅を求めます。
信号のパラメーターとして、サンプリング周波数 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 つの正弦波で構成される信号を作成します。最初の正弦波は位相が の余弦波で、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
の変換は、次式で定義されます。
ここで、
は 1 の n 乗根の 1 つです。
ヒント
関数
fft
の実行時間は、変換の長さに依存します。変換の長さが小さな素因数 (7 以下) のみからなる場合は、素数である場合や大きな素因数からなる場合よりも実行時間がかなり速くなります。ほとんどの
n
の値について、実数入力 DFT の計算時間は複素数入力 DFT の約半分になります。ただし、n
が大きな素因数をもつ場合、速度の差はほとんどありません。ユーティリティ関数
fftw
を使用して、fft
の処理速度を向上できます。この関数は、特定のサイズと次元をもつ FFT の計算に使用されるアルゴリズムの最適化を制御します。
参照
[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++ コードを生成します。
使用上の注意事項および制限事項:
可変サイズ データに関連した制限については、ツールボックス関数のコード生成に対する可変サイズの制限 (MATLAB Coder)を参照してください。
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® コードを生成します。
スレッドベースの環境
MATLAB® の backgroundPool
を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool
を使用してコードを高速化します。
この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。
GPU 配列
Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。
関数 fft
は GPU 配列を部分的にサポートしています。この関数の一部の構文は、入力データを gpuArray
(Parallel Computing Toolbox) として指定すると、GPU で実行されます。使用上の注意事項および制限事項:
虚数部がすべてゼロであっても、出力
Y
は常に複素数です。
詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。
分散配列
Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。
使用上の注意事項および制限事項:
分散配列では、並列 FFT アルゴリズムを使用する代わりに
fft
は単一のワーカーでベクトルを収集し、素数長の FFT を実行します。大きい素数長ベクトルの FFT では、メモリ不足のエラーが発生する場合があります。
詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。
バージョン履歴
R2006a より前に導入
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)