Форум программистов [Powered by Invision Power Board]

Код

unit Unit1;

interface

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

type
  TForm1 = class(TForm)
    btn1: TButton;
    Memo1: TMemo;
    imgpb1: TImage;
    procedure FormCreate(Sender: TObject);
    procedure btn1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses math;

const
   MaxSize = 512;
   Ampl = 511;

var
 int1024sin2pidiv1024,
 int1024cos2pidiv1024: array[0..MaxSize - 1] of integer;
 re,
 im                  : array[0..MaxSize - 1] of integer;

procedure FFT1024();
var i,j,k,L,N, ind, md2,mdd,mdd1, i1,i2:integer;
    a,b,c,d,batr,bati:integer;
    nm : array[0..MaxSize - 1] of integer;
begin
  N:=MaxSize;
  L:=round(ln(N)/ln(2));
  ind:=0;        //перестановка отсчетов в битреверсном порядке как для re[] так и для im[]
  nm[0]:=0;
  for i:=0 to L-1 do begin
    for j:=0 to (1 shl i)-1 do begin
      if ind>nm[ind] then begin
        a:=re[ind];                                //
        re[ind]:=re[nm[ind]];               // перестановка для re[]
        re[nm[ind]]:=a;
        a:=im[ind];
        im[ind]:=im[nm[ind]];             // перестановка для im[]
        im[nm[ind]]:=a;
      end;
      ind:=ind+1;
      nm[ind]:=nm[j]+1 shl (L-1-i);
    end;
  end;     //конец перестановки отсчетов в битреверсном порядке
  md2:=N shr 1;
  for k:=0 to L-1 do begin     //цикл по стадиям  (10 стадий)
    mdd:=md2 shr k;
    mdd1:=md2 shr(L-k-1);
    for i:=0 to mdd-1 do begin  //цикл по группам в каждой стадии  (512x1, 256x2, 128x4, 64x8, 16x32,…, 1x512)
      i1:=2*mdd1*i;
      i2:=2*mdd1*i+mdd1;
      for j:=0 to mdd1-1 do begin   //цикл по парам внутри группы (выполнение бабочки для каждой пары в группе)
        a:=re[j+i1];
        b:=im[j+i1];
        c:=re[j+i2];
        d:=im[j+i2];
        ind:=j*mdd;
        batr:=(c*int1024cos2pidiv1024[ind] +
               d*int1024sin2pidiv1024[ind])div Ampl;
        bati:=(d*int1024cos2pidiv1024[ind] -
               c*int1024sin2pidiv1024[ind])div Ampl;
        re[j+i1]:=a+batr;     //бабочка
        im[j+i1]:=b+bati;    //бабочка
        re[j+i2]:=a-batr;      //бабочка
        im[j+i2]:=b-bati;     //бабочка
      end;
    end;
  end;

end;


procedure IFFT1024();
var i,j,k,L,N, ind, md2,mdd,mdd1,
    i1,i2:integer;
    a,b,c,d,batr,bati:integer;
    nm:array[0..MaxSize - 1] of integer;

begin
  N:= MaxSize;
  L:=round(ln(N)/ln(2));
  ind:=0;            //перестановка отсчетов в битреверсном порядке как для re[] так и для im[]
  nm[0]:=0;
  for i:=0 to L-1 do begin
    for j:=0 to (1 shl i)-1 do begin
      if ind>nm[ind] then begin
        a:=re[ind];
        re[ind]:=re[nm[ind]];
        re[nm[ind]]:=a;
        a:=im[ind];
        im[ind]:=im[nm[ind]];
        im[nm[ind]]:=a;
      end;
      ind:=ind+1;
      nm[ind]:=nm[j]+1 shl (L-1-i);
    end;
  end;                 //конец перестановки отсчетов в битреверсном порядке

  md2:=N shr 1;
  for k:=0 to L-1 do begin      //цикл по стадиям  (10 стадий)
    mdd:=md2 shr k;
    mdd1:=md2 shr(L-k-1);
    for i:=0 to mdd-1 do begin //цикл по группам в каждой стадии  (512x1, 256x2, 128x4, 64x8, 16x32,…, 1x512)
      i1:=2*mdd1*i;
      i2:=2*mdd1*i+mdd1;
      for j:=0 to mdd1-1 do begin //цикл по парам внутри группы (выполнение бабочки для каждой пары в группе)
        a:=re[j+i1];
        b:=im[j+i1];
        c:=re[j+i2];
        d:=im[j+i2];

        ind:=j*mdd;
        batr:=(c*int1024cos2pidiv1024[ind] -
               d*int1024sin2pidiv1024[ind]) div Ampl;
        bati:=(d*int1024cos2pidiv1024[ind] +
               c*int1024sin2pidiv1024[ind]) div Ampl;

        re[j+i1] :=a+batr;
        im[j+i1] :=b+bati;
        re[j+i2] :=a-batr;
        im[j+i2] :=b-bati;
      end;
    end;
  end;

  for k:=0 to N-1 do begin  //деление на 1024 (арифметический сдвиг вправо на 10 дв. разрядов)
    re[k]:=re[k] div N;
    im[k]:=im[k] div N;
  end;

   for i := 0 to MaxSize - 1
    do
      form1.memo1.Lines.Add(IntTostr(re[i]));
    form1.memo1.Lines.SaveToFile('c:\delpi.txt');

end;


procedure TForm1.FormCreate(Sender: TObject);
var
  i : integer;
begin
  for i:=0 to MaxSize - 1
  do
    begin  //массивы заполняются один раз при инициализации.
      int1024cos2pidiv1024[i]:=round(Ampl*cos(2*pi*i/MaxSize)); //один период косинуса:1024 отсчетов
      int1024sin2pidiv1024[i]:=round(Ampl*sin(2*pi*i/MaxSize)); //один период синуса:1024 отсчетов
      re[i] := round(100 * exp(-power(i - 200, 2) / 80) *
                           sin(40 * 2 * pi * i / MaxSize) -
                      80 * exp(-power(i - 210, 2) / 20) *
                           sin(45 * 2 * pi * i / MaxSize) -
                      30 * exp(-power(i - 205, 2) / 10) *
                           sin(41.5 * 2 * pi * i / MaxSize));
      im[i] := 0;
    end;
end;



procedure TForm1.btn1Click(Sender: TObject);
var
  i : integer;
begin
  imgpb1.Canvas.Pen.Color := $0000FF;
  imgpb1.Canvas.MoveTo(0, 120);
  for i := 0 to MaxSize - 1
  do
    imgpb1.Canvas.LineTo(i, 120 - re[i]);

  FFT1024();

{  for i := 0 to MaxSize - 1
  do
    begin
      re[i] := 0;
      im[i] := 0;
    end;
                  }
  //re[4] := 10000;

  IFFT1024();

  imgpb1.Canvas.Pen.Color := $009900;
  imgpb1.Canvas.MoveTo(0, 350);
  for i := 0 to MaxSize - 1
  do
    imgpb1.Canvas.LineTo(i, 400 {350} - re[i]);
                  // round(sqrt(re[i] * re[i] + im[i]*im[i])*2) );
  {
    for i := 0 to MaxSize - 1
    do
      memo1.Lines.Add(IntTostr(re[i]) + '  ' +
                      intToStr(im[i]));
    memo1.Lines.SaveToFile('c:\delpi.txt');
    }
end;

end.