Delphiでお手軽プログラミング

Delphiでお手軽プログラミングメニュー

Delphiで高速フーリエ変換(FFT)を行う


高速フーリエ変換のユニット


//高速フーリエ変換
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 WinFlapTop(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));
  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 WinFlapTop(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);//ブラックマン・ナットール窓関数
  //WinFlapTop(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)法を使うのがいいそうです。

Copyright 2019 Mam