БПФ или FFT

Сергей Викторович · #1 05.03.2012, 23:37

дайте кто может ссылочку на рабочий исходник БПФ (FFT) чтоб массив однострочный реальных значений обрабатывал, желательно чтоб с описанием. Несколько дней ищу, а нормальный не могу найти

angvelem · #2 05.03.2012, 23:42

Что значит:

Цитата:

...БПФ (FFT)...

ChinYan · #3 05.03.2012, 23:51

Быстрое преобразование Фурье?

Сергей Викторович · #4 05.03.2012, 23:52

да быстрое преобразование фурье

Сергей Викторович · #5 05.03.2012, 23:53

Цитата:

Сообщение от ChinYan

Быстрое преобразование Фурье?

да быстрое преобразование фурье

angvelem · #6 06.03.2012, 00:05

Что то ведь нашёл, может доработать "напильником"?

Сергей Викторович · #7 06.03.2012, 00:37

Цитата:

Сообщение от angvelem

Что то ведь нашёл, может доработать "напильником"?

а можно сторонние файлы с расширением .pas как то к проекту добавить?
Можно ответить : можно или нельзя. Если можно сам попробую найти как.
Есть 2 файла с расширением .pas вроде как выполняющих БПФ, а что с ними делать пока не могу разобратся

вот один

Код:

unit fft;
interface
uses Math, Ap, Sysutils;

procedure FastFourierTransform(var a : TReal1DArray;
     nn : Integer;
     InverseFFT : Boolean);

implementation

procedure FastFourierTransform(var a : TReal1DArray;
     nn : Integer;
     InverseFFT : Boolean);
var
    ii : Integer;
    jj : Integer;
    n : Integer;
    mmax : Integer;
    m : Integer;
    j : Integer;
    istep : Integer;
    i : Integer;
    isign : Integer;
    wtemp : Double;
    wr : Double;
    wpr : Double;
    wpi : Double;
    wi : Double;
    theta : Double;
    tempr : Double;
    tempi : Double;
begin
    if InverseFFT then
    begin
        isign := -1;
    end
    else
    begin
        isign := 1;
    end;
    n := 2*nn;
    j := 1;
    ii:=1;
    while ii<=nn do
    begin
        i := 2*ii-1;
        if j>i then
        begin
            tempr := a[j-1];
            tempi := a[j];
            a[j-1] := a[i-1];
            a[j] := a[i];
            a[i-1] := tempr;
            a[i] := tempi;
        end;
        m := n div 2;
        while (m>=2) and (j>m) do
        begin
            j := j-m;
            m := m div 2;
        end;
        j := j+m;
        Inc(ii);
    end;
    mmax := 2;
    while n>mmax do
    begin
        istep := 2*mmax;
        theta := 2*Pi/(isign*mmax);
        wpr := -2.0*sqr(sin(0.5*theta));
        wpi := sin(theta);
        wr := 1.0;
        wi := 0.0;
        ii:=1;
        while ii<=mmax div 2 do
        begin
            m := 2*ii-1;
            jj:=0;
            while jj<=(n-m) div istep do
            begin
                i := m+jj*istep;
                j := i+mmax;
                tempr := wr*a[j-1]-wi*a[j];
                tempi := wr*a[j]+wi*a[j-1];
                a[j-1] := a[i-1]-tempr;
                a[j] := a[i]-tempi;
                a[i-1] := a[i-1]+tempr;
                a[i] := a[i]+tempi;
                Inc(jj);
            end;
            wtemp := wr;
            wr := wr*wpr-wi*wpi+wr;
            wi := wi*wpr+wtemp*wpi+wi;
            Inc(ii);
        end;
        mmax := istep;
    end;
    if InverseFFT then
    begin
        I:=1;
        while I<=2*nn do
        begin
            a[I-1] := a[I-1]/nn;
            Inc(I);
        end;
    end;
end;


end.

angvelem · #8 06.03.2012, 00:49

Всегда было можно.
Данный юнит запрашивает ещё один ap.pas
P.S. Выдели весь введённый код и нажми кнопку [code]. Впредь не забывай это делать всегда.

Сергей Викторович · #9 06.03.2012, 00:50

вот второй

Код:

unit Ap;

interface

uses Math;


// constants

const
    MachineEpsilon = 5E-16;
    MaxRealNumber = 1E300;
    MinRealNumber = 1E-300;

// arrays

    Complex = record
        X, Y: Double;
    end;

    TInteger1DArray     = array of LongInt;
    TReal1DArray        = array of Double;
    TComplex1DArray     = array of Complex;
    TBoolean1DArray     = array of Boolean;

    TInteger2DArray     = array of array of LongInt;
    TReal2DArray        = array of array of Double;
    TComplex2DArray     = array of array of Complex;
    TBoolean2DArray     = array of array of Boolean;


/////////////////////////////////////////////////////////////////////////
// Functions/procedures
/////////////////////////////////////////////////////////////////////////
function AbsReal(X : Extended):Extended;
function AbsInt (I : Integer):Integer;
function RandomReal():Extended;
function RandomInteger(I : Integer):Integer;
function Sign(X:Extended):Integer;

function DynamicArrayCopy(const A: TInteger1DArray):TInteger1DArray;overload;
function DynamicArrayCopy(const A: TReal1DArray):TReal1DArray;overload;
function DynamicArrayCopy(const A: TComplex1DArray):TComplex1DArray;overload;
function DynamicArrayCopy(const A: TBoolean1DArray):TBoolean1DArray;overload;

function DynamicArrayCopy(const A: TInteger2DArray):TInteger2DArray;overload;
function DynamicArrayCopy(const A: TReal2DArray):TReal2DArray;overload;
function DynamicArrayCopy(const A: TComplex2DArray):TComplex2DArray;overload;
function DynamicArrayCopy(const A: TBoolean2DArray):TBoolean2DArray;overload;

function AbsComplex(const Z : Complex):Double;
function Conj(const Z : Complex):Complex;
function CSqr(const Z : Complex):Complex;

function C_Complex(const X : Double):Complex;
function C_Opposite(const Z : Complex):Complex;
function C_Add(const Z1 : Complex; const Z2 : Complex):Complex;
function C_Mul(const Z1 : Complex; const Z2 : Complex):Complex;
function C_AddR(const Z1 : Complex; const R : Double):Complex;
function C_MulR(const Z1 : Complex; const R : Double):Complex;
function C_Sub(const Z1 : Complex; const Z2 : Complex):Complex;
function C_SubR(const Z1 : Complex; const R : Double):Complex;
function C_RSub(const R : Double; const Z1 : Complex):Complex;
function C_Div(const Z1 : Complex; const Z2 : Complex):Complex;
function C_DivR(const Z1 : Complex; const R : Double):Complex;
function C_RDiv(const R : Double; const Z2 : Complex):Complex;
function C_Equal(const Z1 : Complex; const Z2 : Complex):Boolean;
function C_NotEqual(const Z1 : Complex; const Z2 : Complex):Boolean;
function C_EqualR(const Z1 : Complex; const R : Double):Boolean;
function C_NotEqualR(const Z1 : Complex; const R : Double):Boolean;

implementation

/////////////////////////////////////////////////////////////////////////
// Functions/procedures
/////////////////////////////////////////////////////////////////////////
function AbsReal(X : Extended):Extended;
begin
    Result := Abs(X);
end;

function AbsInt (I : Integer):Integer;
begin
    Result := Abs(I);
end;

function RandomReal():Extended;
begin
    Result := Random;
end;

function RandomInteger(I : Integer):Integer;
begin
    Result := Random(I);
end;

function Sign(X:Extended):Integer;
begin
    if X>0 then
        Result := 1
    else if X<0 then
        Result := -1
    else
        Result := 0;
end;

/////////////////////////////////////////////////////////////////////////
// dynamical arrays copying
/////////////////////////////////////////////////////////////////////////
function DynamicArrayCopy(const A: TInteger1DArray):TInteger1DArray;overload;
var
    I:  Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
        Result[i]:=A[i];
end;

function DynamicArrayCopy(const A: TReal1DArray):TReal1DArray;overload;
var
    I:  Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
        Result[i]:=A[i];
end;

function DynamicArrayCopy(const A: TComplex1DArray):TComplex1DArray;overload;
var
    I:  Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
        Result[i]:=A[i];
end;

function DynamicArrayCopy(const A: TBoolean1DArray):TBoolean1DArray;overload;
var
    I:  Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
        Result[i]:=A[i];
end;

function DynamicArrayCopy(const A: TInteger2DArray):TInteger2DArray;overload;
var
    I,J:    Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
    begin
        SetLength(Result[i], High(A[i])+1);
        for J:=Low(A[i]) to High(A[i]) do
            Result[I,J]:=A[I,J];
    end;
end;

function DynamicArrayCopy(const A: TReal2DArray):TReal2DArray;overload;
var
    I,J:    Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
    begin
        SetLength(Result[i], High(A[i])+1);
        for J:=Low(A[i]) to High(A[i]) do
            Result[I,J]:=A[I,J];
    end;
end;

function DynamicArrayCopy(const A: TComplex2DArray):TComplex2DArray;overload;
var
    I,J:    Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
    begin
        SetLength(Result[i], High(A[i])+1);
        for J:=Low(A[i]) to High(A[i]) do
            Result[I,J]:=A[I,J];
    end;
end;

function DynamicArrayCopy(const A: TBoolean2DArray):TBoolean2DArray;overload;
var
    I,J:    Integer;
begin
    SetLength(Result, High(A)+1);
    for I:=Low(A) to High(A) do
    begin
        SetLength(Result[i], High(A[i])+1);
        for J:=Low(A[i]) to High(A[i]) do
            Result[I,J]:=A[I,J];
    end;
end;

/////////////////////////////////////////////////////////////////////////
// complex numbers
/////////////////////////////////////////////////////////////////////////
function AbsComplex(const Z : Complex):Double;
var
    W : Double;
    XABS : Double;
    YABS : Double;
    V : Double;
begin
    XABS := AbsReal(Z.X);
    YABS := AbsReal(Z.Y);
    W := Max(XABS, YABS);
    V := Min(XABS, YABS);
    if V=0 then
    begin
        Result := W;
    end
    else
    begin
        Result := W*SQRT(1+Sqr(V/W));
    end;
end;


function Conj(const Z : Complex):Complex;
begin
    Result.X := Z.X;
    Result.Y := -Z.Y;
end;


function CSqr(const Z : Complex):Complex;
begin
    Result.X := Sqr(Z.X)-Sqr(Z.Y);
    Result.Y := 2*Z.X*Z.Y;
end;


function C_Complex(const X : Double):Complex;
begin
    Result.X := X;
    Result.Y := 0;
end;


function C_Opposite(const Z : Complex):Complex;
begin
    Result.X := -Z.X;
    Result.Y := -Z.Y;
end;


function C_Add(const Z1 : Complex; const Z2 : Complex):Complex;
begin
    Result.X := Z1.X+Z2.X;
    Result.Y := Z1.Y+Z2.Y;
end;


function C_Mul(const Z1 : Complex; const Z2 : Complex):Complex;
begin
    Result.X := Z1.X*Z2.X-Z1.Y*Z2.Y;
    Result.Y := Z1.X*Z2.Y+Z1.Y*Z2.X;
end;


function C_AddR(const Z1 : Complex; const R : Double):Complex;
begin
    Result.X := Z1.X+R;
    Result.Y := Z1.Y;
end;


function C_MulR(const Z1 : Complex; const R : Double):Complex;
begin
    Result.X := Z1.X*R;
    Result.Y := Z1.Y*R;
end;


function C_Sub(const Z1 : Complex; const Z2 : Complex):Complex;
begin
    Result.X := Z1.X-Z2.X;
    Result.Y := Z1.Y-Z2.Y;
end;


function C_SubR(const Z1 : Complex; const R : Double):Complex;
begin
    Result.X := Z1.X-R;
    Result.Y := Z1.Y;
end;


function C_RSub(const R : Double; const Z1 : Complex):Complex;
begin
    Result.X := R-Z1.X;
    Result.Y := -Z1.Y;
end;


function C_Div(const Z1 : Complex; const Z2 : Complex):Complex;
var
    A : Double;
    B : Double;
    C : Double;
    D : Double;
    E : Double;
    F : Double;
begin
    A := Z1.X;
    B := Z1.Y;
    C := Z2.X;
    D := Z2.Y;
    if AbsReal(D)<AbsReal(C) then
    begin
        E := D/C;
        F := C+D*E;
        Result.X := (A+B*E)/F;
        Result.Y := (B-A*E)/F;
    end
    else
    begin
        E := C/D;
        F := D+C*E;
        Result.X := (B+A*E)/F;
        Result.Y := (-A+B*E)/F;
    end;
end;


function C_DivR(const Z1 : Complex; const R : Double):Complex;
begin
    Result.X := Z1.X/R;
    Result.Y := Z1.Y/R;
end;


function C_RDiv(const R : Double; const Z2 : Complex):Complex;
var
    A : Double;
    C : Double;
    D : Double;
    E : Double;
    F : Double;
begin
    A := R;
    C := Z2.X;
    D := Z2.Y;
    if AbsReal(D)<AbsReal(C) then
    begin
        E := D/C;
        F := C+D*E;
        Result.X := A/F;
        Result.Y := -A*E/F;
    end
    else
    begin
        E := C/D;
        F := D+C*E;
        Result.X := A*E/F;
        Result.Y := -A/F;
    end;
end;


function C_Equal(const Z1 : Complex; const Z2 : Complex):Boolean;
begin
    Result := (Z1.X=Z2.X) and (Z1.Y=Z2.Y);
end;


function C_NotEqual(const Z1 : Complex; const Z2 : Complex):Boolean;
begin
    Result := (Z1.X<>Z2.X) or (Z1.Y<>Z2.Y);
end;

function C_EqualR(const Z1 : Complex; const R : Double):Boolean;
begin
    Result := (Z1.X=R) and (Z1.Y=0);
end;

function C_NotEqualR(const Z1 : Complex; const R : Double):Boolean;
begin
    Result := (Z1.X<>R) or (Z1.Y<>0);
end;


end.

angvelem · #10 06.03.2012, 01:05

Выдели весь введённый код и нажми кнопку [code] - (#). Впредь не забывай это делать ВСЕГДА.

angvelem · #11 06.03.2012, 02:13

Возиожно этот вариант проще:

Код:

unit cplxfft2;

interface

type
  PScalar = ^TScalar;
  TScalar = Extended;

  PScalars = ^TScalars;
  TScalars = array[0..High(Integer) div SizeOf(TScalar) - 1] of TScalar;

const
  TrigTableDepth : word = 0;
  CosTable : PScalars = NIL;
  SinTable : PScalars = NIL;

procedure InitTrigTables(Depth : Word);

procedure FFT(Depth : Word; SrcR, SrcI : PScalars; DestR, DestI : PScalars);

{Перед вызовом Src и Dest ТРЕБУЕТСЯ распределение (integer(1) shl Depth) * SizeOf(TScalar) байт памяти!}

implementation

procedure DoFFT(Depth : Word; SrcR, SrcI : PScalars; SrcSpacing : Word;
 DestR, DestI : PScalars);
{рекурсивная часть, вызываемая при готовности FFT}
var
  J, N : Integer;
  TempR, TempI : TScalar;
  Sh : Word;
  c, s : Extended;
begin
  if Depth = 0 then
  begin
    DestR^[0] := SrcR^[0];
    DestI^[0] := SrcI^[0];
    Exit;
  end;

  N := Integer(1) shl (Depth - 1);

  DoFFT(Depth - 1, SrcR, SrcI, SrcSpacing * 2, DestR, DestI);
  DoFFT(Depth - 1, @SrcR^[srcSpacing], @SrcI^[SrcSpacing], SrcSpacing * 2, @DestR^[N], @DestI^[N]);

  Sh := TrigTableDepth - Depth;

  for J := 0 to N - 1 do
  begin
    c := CosTable^[J shl Sh];
    s := SinTable^[J shl Sh];

    TempR := c * DestR^[J + N] - s * DestI^[J + N];
    TempI := c * DestI^[J + N] + s * DestR^[J + N];

    DestR^[J + N] := DestR^[J] - TempR;
    DestI^[J + N] := DestI^[J] - TempI;

    DestR^[J] := DestR^[J] + TempR;
    DestI^[J] := DestI^[J] + TempI;
  end;
end;

procedure FFT(Depth : Word; SrcR, SrcI : PScalars; DestR, DestI : PScalars);
var
  J, N : Integer;
  Normalizer : Extended;
begin
  N := Integer(1) shl Depth;

  if Depth TrigTableDepth then
    InitTrigTables(Depth);

  DoFFT(Depth, SrcR, SrcI, 1, DestR, DestI);

  Normalizer := 1 / sqrt(N);

  for J := 0 to N - 1 do
  begin
    DestR^[J] := DestR^[J] * Normalizer;
    DestI^[J] := DestI^[J] * Normalizer;
  end;
end;

procedure InitTrigTables(Depth : Word);
var
  J, N : Integer;
begin
  N := Integer(1) shl Depth;
  ReAllocMem(CosTable, N * SizeOf(TScalar));
  ReAllocMem(SinTable, N * SizeOf(TScalar));

  for J := 0 to N - 1 do
  begin
    CosTable^[J] := cos(-(2 * Pi) * J / N);
    SinTable^[J] := sin(-(2 * Pi) * J / N);
  end;
  TrigTableDepth := Depth;
end;

finalization
  ReAllocMem(CosTable, 0);
  ReAllocMem(SinTable, 0);

end.

unit demofrm;

interface

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

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    Edit1: TEdit;
    Label1: TLabel;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.DFM}

uses
  MMSystem;

procedure TForm1.Button1Click(Sender: TObject);
var
  SR, SI, DR, DI : PScalars;
  J, D, N : Integer;
  st, et : Longint;
  norm : Extended;
begin
  D := StrToIntDef(edit1.text, -1);
  if d < 1 then
    raise
      exception.Create('глубина рекурсии должны быть положительным целым числом');

  N := Integer(1) shl D;

  GetMem(SR, N * SizeOf(TScalar));
  GetMem(SI, N * SizeOf(TScalar));
  GetMem(DR, N * SizeOf(TScalar));
  GetMem(DI, N * SizeOf(TScalar));

  for J := 0 to N - 1 do
  begin
    SR^[J] := Random;
    SI^[J] := Random;
  end;

  st := TimeGetTime;
  FFT(D, SR, SI, DR, DI);

  et := TimeGetTime;

  memo1.Lines.Add('N = ' + inttostr(N));
  memo1.Lines.Add('норма ожидания: ' + #9 + FloatToStr(N * 2 / 3));

  norm := 0;
  for J := 0 to N - 1 do
    norm := norm + SR^[J] * SR^[J] + SI^[J] * SI^[J];
  memo1.Lines.Add('норма данных: ' + #9 + FloatToStr(norm));

  norm := 0;
  for J := 0 to N - 1 do
    norm := norm + DR^[J] * DR^[J] + DI^[J] * DI^[J];
  memo1.Lines.Add('норма FT: ' + #9#9 + FloatToStr(norm));

  memo1.Lines.Add('Время расчета FFT: ' + #9 + inttostr(et - st));
  memo1.Lines.add('');

  FreeMem(SR, N * SizeOf(TScalar));
  FreeMem(SI, N * SizeOf(TScalar));
  FreeMem(DR, N * SizeOf(TScalar));
  FreeMem(DI, N * SizeOf(TScalar));
end;

end.

Опции темы	Поиск в этой теме
Версия для печати Отправить на Email	Поиск в этой теме: Расширенный поиск
Опции просмотра
Линейный вид Комбинированный вид Древовидный вид

	Сайт	Форум	FAQ	Соглашения	Прочее
	Новости Исходники Компоненты Статьи Добавить исходник Поиск	Правила Поиск Сообщения за сегодня Пользователи Кто на форуме?	Delphi FAQ - 5000 статей DRKB - 4500 статей (ZIP-архив) Справочник функций и процедур Delphi Delphi Programming Guide Indy in Depth. Глубины Indy	Конфиденциальность Соглашение пользователя Правила обработки cookie	О сайте Donate Контакты