unit1.pas

进行检测的频扑图的例子.设计到算法和显示fft变化.

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls,Math;

type
  compx=record
    real:Double;
    Imag:Double;
  end;
  TForm1 = class(TForm)
    Button1: TButton;
    procedure Button1Click(Sender: TObject);
  private
    function EE(b1,b2:compx):compx;
    procedure FFT(xin:array of compx;N:integer);
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

{ TForm1 }

function TForm1.EE(b1, b2: compx): compx;
var
  b3:compx;
begin
  b3.real:=b1.real*b2.real-b1.Imag*b2.Imag;
  b3.Imag:=b1.real*b2.Imag+b1.Imag*b2.real;

  result:=b3;
end;

procedure TForm1.FFT(xin: array of compx; N: integer);
var
  f,m,LH,nm,i,k,j,L:integer;
  p , ps:double;
  le,B,ip:integer;
  pi:real;
  v,w,t:compx;
begin
  LH:=N div 2;
  f:=N;

  for m:=1 to high(integer) do  //(m=1;(f=f/2)!=1;m++){;}
  begin
    f:=f div 2;
    if f=1 then break;
  end;

  nm:=N-2;
  j:=N div 2;

  //*变址运算*/
  for i:=1 to nm do
  begin
    if i<j then
    begin
      t:=xin[j];
      xin[j]:=xin[i];
      xin[i]:=t;
    end;

    k:=LH;
    while j>=k do
    begin
      j:=j-k;
      k:=k div 2;
    end;
    j:=j+k;
  end;

  for L:=1 to m do  //(L=1;L<=m;L++)
  begin
    le:=Trunc(power(2,L));
    B:=le div 2;
    pi:=3.14159;
    for j:=0 to B-1 do //(j=0;j<=B-1;j++)
    begin
      p:=power(2,m-L)*j;
      ps:=2*pi/N*p;
      w.real:=cos(ps);
      w.imag:=-sin(ps);
      for i:=j to N-1 do  //(i=j;i<=N-1;i=i+le)
      begin
        ip:=i+B;
        t:=EE(xin[ip],w);
        xin[ip].real:=xin[i].real-t.real;
        xin[ip].imag:=xin[i].imag-t.imag;
        xin[i].real:=xin[i].real+t.real;
        xin[i].imag:=xin[i].imag+t.imag;
      end;
    end;
  end;
end;

procedure TForm1.Button1Click(Sender: TObject);
const
  pp=3.14159;
var
  result:array[0..257] of real;
  s:array[0..257] of compx;
  Num:integer;
  i:integer;
begin
  Num:=16;
  for i:=0 to 16 do  //(i=0;i<16;i++)
  begin
    s[i].real:=sin(pp*i/32);
    s[i].imag:=0;
  end;
  FFT(s,Num);

  for i:=0 to 16 do  //(i=0;i<16;i++)
  begin
    showmessage(Format('%.4f',[s[i].real]));
    showmessage(format('+%.4fj\n',[s[i].imag]));
    result[i]:=sqrt(power(s[i].real,2)+power(s[i].imag,2));
  end;
end;

end.