高速フーリエ変換(FFT) ～Delphiソースコード集
高速フーリエ変換のユニット
//高速フーリエ変換
unit UFFT;
interface
uses System.SysUtils,System.Math;
type
TFFTData=array of Extended;
//高速フーリエ変換(入力信号,出力実部,出力虚部)
procedure fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
procedure _fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
//*******************窓関数*****************************
//ハン窓関数
procedure WinHanning(var data:TFFTData);
//ハミング窓関数
procedure WinHamming(var data:TFFTData);
//ガウス窓関数
procedure WinGauss(var data:TFFTData;m:integer=1);
//ブラックマンハリス窓関数
procedure WinBlackmanHarris(var data:TFFTData);
//ブラックマンナットール窓関数
procedure WinBlackmanNuttall(var data:TFFTData);
//フラットトップ窓関数
procedure WinFlatTop(var data:TFFTData);
//半正弦波窓関数
procedure WinHalfSin(var data:TFFTData);
implementation
procedure fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
var i :Integer;
InN:Integer;//入力データ数
n :Integer;//補正後データ数
begin
InN:=Length(InRe);
//データ数が2の乗数に満たない場合は0のデータを追加する
i:=1;
while InN > Power(2,i) do inc(i);
n:=trunc(IntPower(2,i));
if InN < n then
begin
setlength(InRe,n);
for i := InN to n-1 do InRe[i]:=0;
end;
//高速フーリエ変換
_fft(InRe,OutRe,OutIm);
end;
procedure _fft(InRe:TFFTData;var OutRe,OutIm:TFFTData);
var n:Integer;
i:integer;
InIm:TFFTData; //複素数の虚部
ct1,ct2,ct3:integer;
TmpRe,TmpIm:extended;
nfft:array[0..3] of integer;
fcos,fsin:TFFTData;
tmp:extended;
noblk:integer;
cntb:array[0..1] of integer;
begin
n:=Length(InRe);
setlength(InIm,n);
for i := 0 to n-1 do InIm[i]:=0;
ct2:=1;
for ct1 := 1 to length(InIm)-2 do
begin
TmpRe:=0;
TmpIm:=0;
if ct1<ct2 then
begin
TmpRe:=InRe[ct1-1];
InRe[ct1-1]:=InRe[ct2-1];
InRe[ct2-1]:=TmpRe;
TmpIm:=InIm[ct1-1];
InIm[ct1-1]:=InIm[ct2-1];
InIm[ct2-1]:=TmpIm;
end;
ct3:=length(InIm) div 2;
while ct3<ct2 do
begin
ct2:=ct2-ct3;
ct3:=ct3 div 2;
end;
ct2:=ct2+ct3;
end;
//誤差調整
nfft[0]:=floor(Log2(length(InIm))/Log2(2)+0.0000001);
SetLength(fcos,n);
SetLength(fsin,n);
fcos[0]:=1;
fsin[0]:=0;
for ct1 := 1 to nfft[0] do
begin
nfft[2]:=floor(System.math.Power(2,ct1));
nfft[1]:=n div nfft[2];
nfft[3]:=nfft[2] div 2;
for ct2 := 1 to nfft[3] do
begin
tmp:=-Pi/nfft[3]*ct2;
fcos[ct2]:=cos(tmp);
fsin[ct2]:=sin(tmp);
end;
for ct2 := 1 to nfft[1] do
begin
noblk:=nfft[2]*(ct2-1);
for ct3 := 0 to nfft[3]-1 do
begin
cntb[0]:=noblk+ct3;
cntb[1]:=cntb[0]+nfft[3];
TmpRe:=InRe[cntb[1]]*fcos[ct3]-InIm[cntb[1]]*fsin[ct3];
TmpIm:=InIm[cntb[1]]*fcos[ct3]+InRe[cntb[1]]*fsin[ct3];
InRe[cntb[1]]:=InRe[cntb[0]]-TmpRe;
InIm[cntb[1]]:=InIm[cntb[0]]-TmpIm;
InRe[cntb[0]]:=InRe[cntb[0]]+TmpRe;
InIm[cntb[0]]:=InIm[cntb[0]]+TmpIm;
end;
end;
end;
setlength(OutRe,length(InRe) div 2);
setlength(OutIm,length(InIm) div 2);
for i := 0 to (n div 2)-1 do
begin
OutRe[i]:=InRe[i]; //実部
OutIm[i]:=InIm[i]; //虚部
end;
end;
//ハン窓関数
procedure WinHanning(var data:TFFTData);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=(
0.5 - 0.5*Cos(2*Pi*i/(n-1))
)*data[i];
end;
end;
//ハミング窓関数
procedure WinHamming(var data:TFFTData);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=(
0.54 - 0.46 * Cos(2*Pi*i/(n-1))
)*data[i];
end;
end;
//ガウス窓関数
procedure WinGauss(var data:TFFTData;m:integer=1);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=Exp(
-2 * power(m,2) / power(n-1,2) * power(i-(n-1)/2,2)
)*data[i];
end;
end;
//ブラックマンハリス窓関数
procedure WinBlackmanHarris(var data:TFFTData);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=(0.35875-0.48829*cos(2*Pi*i/(n-1))
+0.14128*cos(4*Pi*i/(n-1))
-0.01168*cos(6*Pi*i/(n-1))
)*data[i];
end;
end;
//ブラックマンナットール窓関数
procedure WinBlackmanNuttall(var data:TFFTData);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=(0.3635819-0.4891775*cos(2*Pi*i/(n-1))
+0.1365995*cos(4*Pi*i/(n-1))
-0.0106411*cos(6*Pi*i/(n-1))
)*data[i];
end;
end;
//フラットトップ窓関数
procedure WinFlatTop(var data:TFFTData);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=(1-1.930*Cos(2*Pi*i/(n-1))
+1.290*Cos(4*Pi*i/(n-1))
-0.388*Cos(6*Pi*i/(n-1))
+0.032*Cos(8*Pi*i/(n-1))
)*data[i];
end;
end;
//半正弦波窓関数
procedure WinHalfSin(var data:TFFTData);
var i,n:integer;
begin
n:=length(data);
for i := 0 to n-1 do
begin
data[i]:=Sin(pi*i/(n-1))*data[i];
end;
end;
end.
参考URL https://gijyutsu-keisan.com/excel/numcal/fourier/fourier_2_2.php高速フーリエ変換は、入力データ数が2のn乗個であることが前提の為、関数fft内で入力データ数が2のn乗個になるように調整してから関数_fftを呼び出しています。 また、高速フーリエ変換は、入力データの開始と終了が連続していることが前提となっており、そうでない場合は窓関数を使用することになります。
高速フーリエ変換を使ってみる
unit Unit1;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.ExtCtrls;
type
TForm1 = class(TForm)
Button1: TButton;
Button2: TButton;
procedure Button2Click(Sender: TObject);
procedure Button1Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses UFFT,system.math;
procedure TForm1.Button1Click(Sender: TObject);
var InRe:TFFTData;
//複素数の実部と虚部 i = √-1
// z = Re + Im * i
OutRe,OutIm:TFFTData;
n,i:integer;
bs:extended;// サンプリング間隔(秒)
hz:extended;// サンプリング周波数（Hz）= 1/bs
bhz:extended;
OutHz:extended; //周波数
OutAmp:extended; //振幅
begin
n:=200; //サンプル数
bs:=0.01;//サンプリング間隔（10ミリ秒毎のデータ）
hz:=1/bs;//サンプリング周波数（100Hz）
//テストデータの作成
setlength(InRe,n);
for i := 0 to n-1 do
begin
InRe[i]:=4*sin(i*bs*2*pi*10)+ //10Hz
5*sin(i*bs*2*pi*20)+ //20Hz
6*sin(i*bs*2*pi*40); //40Hz
end;
//窓関数処理が必要な（データの開始と終了が連続でない）場合は
//いずれかの窓関数処理を行う
WinHanning(InRe); //ハン窓関数
//WinHamming(InRe); //ハミング窓関数
//WinGauss(InRe,1); //ガウス窓関数 1～3σ
//WinBlackmanHarris(InRe); //ブラックマン・ハリス窓関数
//WinBlackmanNuttall(InRe);//ブラックマン・ナットール窓関数
//WinFlatTop(InRe); //フラットトップ窓関数
//WinHalfSin(InRe); //半正弦波窓関数
//高速フーリエ変換を行う
fft(InRe,OutRe,OutIm);
bhz:=hz/(length(OutRe)*2);
for i := 0 to length(OutRe)-1 do
begin
//周波数
OutHz :=bhz*i;
//振幅 = √(実部^2 + 虚部^2)
OutAmp:=Sqrt(power(OutRe[i],2)+power(OutIm[i],2));
end;
end;
end.
サンプリング周波数Fsの場合、ナイキスト周波数 = Fs/2
FFTはナイキスト周波数を境にして左右対称の鏡像データが得られるので、FFTで得られる最大周波数はナイキスト周波数となる。
振幅はFFTで得られた実数部の2乗と虚数部の2乗の和を平方根した値。
ちなみに位相はFFTで得られた虚数部÷実数部 をアークタンジェントした値。
上記例ではサンプリング間隔10ms(サンプリング周波数100Hz、ナイキスト周波数50Hz)で 振幅4の10Hz、振幅5の20Hz、振幅6の40Hzの合成データを与えています。
また、FFTでは入力データの最初と最後が連続していることが前提なので、 この例ではハン窓関数（ハニング窓関数）を用いてデータを補正してからFFT関数を呼び出しています。
データが連続（一定間隔でデータを取得）していれば高速フーリエ変換で周波数解析ができるのですが、 一様にサンプリングされず欠損データがある場合は Lomb-Scargle ピリオドグラム(Lomb-Scargle Periodogram)法 を使うのがいいそうです。