fourier.pas

以面向对象方法实现的数值算法类库

{-------------------------------------------------------------------------------
2004.05.08 用于实现傅里叶变换
-------------------------------------------------------------------------------}
unit Fourier;

interface

uses
  SysUtils, Classes, Math, Mathtypes, Vector, Complexs;

type

  TOnGetDataEvent = procedure(index : integer; var Value : TComplex) of Object;

  EFastFourierError = class(Exception);

  TFourier = class(TComponent)
  private
    { Private declarations }
    // 数据缓冲区
    FBufferReal     : TVector;
    FBufferImage    : TVector;
    FOnGetData      : TOnGetDataEvent;

    function  IsPowerOfTwo ( x: word ): boolean;
    procedure FourierAnalysis_2(FFT: Boolean = True);

  protected
    { Protected declarations }
    function  GetTransformedData(idx : integer) : TComplex;     // 复数结果
    function  GetTransformedMagnitude(idx : integer) : TFloat;  // 复数的模
    function  GetTransformedPhase(idx : integer) : TFloat;      // 复数的相位
    function  GetTransformedComplexArray: TComplexArray;        // 返回结果复数数组
    function  GetTransformedPowerArray: TVector;                // 返回功率数组
    function  GetTransformedFrequenceArray: TVector;            // 返回对应频率数组

  public
    { Public declarations }
    // 原始数据
    procedure SetSourceData(aX: Array of Double; aY: Array of Double); overload;
    procedure SetSourceData(aX: TFloatArray);     overload;
    procedure SetSourceData(aX: TVector);         overload;
    procedure SetSourceData(aX, aY: TVector);     overload;
    procedure SetSourceData(aX: TComplexArray);   overload;

    // 执行积分变换
    procedure FFT;
    procedure iFFT;
    // 变换结果的输出
    property  Data[idx : integer]       : TComplex  read GetTransformedData;
    property  Magnitude[idx : integer]  : TFloat    read GetTransformedMagnitude;
    // 构造函数与析构函数
    constructor
              create(AOwner : TComponent); override;
    destructor
              destroy;                     override;

  published
    { Published declarations }
    property  OnGetData   : TOnGetDataEvent read FOnGetData write FOnGetData;
    {
    property  ComplexArray: TComplexArray read GetTransformedComplexArray;
    }
    property  RealArray: TVector read FBufferReal;
    property  ImageArray: TVector read FBufferImage;
    property  PowerArray: TVector read GetTransformedPowerArray;
    property  FrequenceArray: TVector read GetTransformedFrequenceArray;
  end;

procedure Register;

implementation

procedure Register;
begin
  RegisterComponents('Lxp', [TFourier]);
end;

{ TFourierTransform }

constructor TFourier.create(AOwner: TComponent);
begin
    inherited create(AOwner);
    FBufferReal := TVector.Create;
    FBufferImage := TVector.Create;
end;

destructor TFourier.destroy;
begin
    FBufferReal := Nil;
    FBufferImage := Nil;
    inherited Destroy;
end;

procedure TFourier.FFT;
begin
    FourierAnalysis_2;//(FBufferReal.Values, FBufferImage.Values);
end;

procedure TFourier.FourierAnalysis_2(FFT: Boolean);
var
  i, i0, i1, j, l1, ns, k, n, m: integer;
  s, c, s1, c1, sc, x1, y1, t: real;
begin
  // 数据点数
  n := Length(FBufferReal.Values);
  m := Length(FBufferImage.Values);
  if Not IsPowerOfTwo(m) then
    raise EFastFourierError.Create('数据序列长度不对！');

  if m<>n then begin
    raise EFastFourierError.Create('数据序列长度不同！');
    Exit;
  end;

  sc := Pi;

  j := 0;
  for i:=0 to n-2 do begin
    if i<j then begin
      t := FBufferReal.Values[i];
      FBufferReal.Values[i] := FBufferReal.Values[j];
      FBufferReal.Values[j] := t;

      t := FBufferImage.Values[i];
      FBufferImage.Values[i] := FBufferImage.Values[j];
      FBufferImage.Values[j] := t;
    end;
    k := n div 2;
    while k<=j do begin
      j := j - k;
      k := k div 2;
    end;
    j := j + k;
  end;

  ns := 1;
  while ns<= n/2 do begin      // 主循环
    c1 := cos(sc);
    if FFT then s1 := sin(-1.0 * sc) else s1 := sin(1.0 * sc);
    c := 1.0;
    s := 0.0;
    for l1:=0 to ns-1 do begin
      i0 := l1;
      while i0 <n do begin
        i1 := i0 + ns;
        x1 := FBufferReal.Values[i1] * c - FBufferImage.Values[i1] * s;
        y1 := FBufferImage.Values[i1] * c + FBufferReal.Values[i1] * s;
        FBufferReal.Values[i1] := FBufferReal.Values[i0] - x1;
        FBufferImage.Values[i1] := FBufferImage.Values[i0] - y1;
        FBufferReal.Values[i0] := FBufferReal.Values[i0] + x1;
        FBufferImage.Values[i0] := FBufferImage.Values[i0] + y1;
        i0 := i0 + 2 * ns;
      end;
      t := c1 * c - s1 * s;
      s := s1 * c + c1 * s;
      c := t;
    end;
    ns := ns * 2;
    sc := sc / 2;
  end;

  if FFT then begin
    for i := 0 to n-1 do begin
      FBufferReal.Values[i] := FBufferReal.Values[i] / n;
      FBufferImage.Values[i] := FBufferImage.Values[i] / n;
    end;
  end;
end;

function TFourier.GetTransformedData(idx: integer): TComplex;
begin
    Result.Re := FBufferReal.Values[idx];
    Result.Im := FBufferImage.Values[idx];
end;

function TFourier.GetTransformedMagnitude(idx: integer): TFloat;
Var
    aC: TComplex;
begin
    aC.Re := FBufferReal.Values[idx];
    aC.Im := FBufferImage.Values[idx];
    Result := ComplexMagnitude(aC);
end;

procedure TFourier.iFFT;
begin
    FourierAnalysis_2(False);
end;

function TFourier.IsPowerOfTwo(x: word): boolean;
var
    i, y:  word;
begin
    y := 2;
    for i := 1 to 31 do begin
        if x = y then begin
            IsPowerOfTwo := TRUE;
            exit;
        end;
        y := y SHL 1;
    end;

    IsPowerOfTwo := FALSE;
end;

procedure TFourier.SetSourceData(aX: TComplexArray);
Var
    i,N: Integer;
begin
    N := Length(aX);
    FBufferReal.Resize(N);
    FBufferImage.Resize(N);
    for i:=0 to N-1 do begin
        FBufferReal.Values[i] := aX[i].Re;
        FBufferImage.Values[i] := aX[i].Im;
    end;
end;

function TFourier.GetTransformedPhase(idx: integer): TFloat;
Var
    aC: TComplex;
begin
    aC.Re := FBufferReal.Values[idx];
    aC.Im := FBufferImage.Values[idx];
    Result := ComplexPhase(aC);
end;

function TFourier.GetTransformedComplexArray: TComplexArray;
Var
    i,N: Integer;
begin
    N := FBufferReal.Count;
    Result := nil;
    SetLength(Result, N);
    for i:=0 to N-1 do begin
      Result[i].Re := FBufferReal.Values[i];
      Result[i].Im := FBufferImage.Values[i];
    end;
end;

function TFourier.GetTransformedPowerArray: TVector;
Var
    i,N: Integer;
begin
    N := FBufferReal.Count div 2;
    Result := TVector.Create;
    Result.Resize(N);
    for i:=0 to N-1 do begin
      Result.Values[i] := FBufferReal.Values[i] * FBufferReal.Values[i] +
                   FBufferImage.Values[i] * FBufferImage.Values[i];
    end;
end;

function TFourier.GetTransformedFrequenceArray: TVector;
Var
    i,N: Integer;
    Nyquist: double;
begin
    Nyquist := 0.5;
    N := FBufferReal.Count div 2;
    Result := TVector.Create;
    Result.Resize(N);
    for i:=0 to N-1 do begin
      Result.Values[i] := i / N * Nyquist;
    end;
end;

procedure TFourier.SetSourceData(aX: TVector);
Var
    i,N: Integer;
begin
    N := aX.Count;
    FBufferReal.Resize(N);
    FBufferImage.Resize(N);
    for i:=0 to N-1 do begin
      FBufferReal.Values[i] := aX.Values[i];
      FBufferImage.Values[i] := 0;
    end;
end;

procedure TFourier.SetSourceData(aX, aY: TVector);
Var
    i,N: Integer;
begin
    N := aX.Count;
    FBufferReal.Resize(N);
    FBufferImage.Resize(N);
    for i:=0 to N-1 do begin
      FBufferReal.Values[i] := aX.Values[i];
      FBufferImage.Values[i] := aY.Values[i];
    end;
end;

procedure TFourier.SetSourceData(aX: Array of Double; aY: Array of Double);
Var
    i,N,M: Integer;
begin
    N := Length(aX);
    M := Length(aY);
    FBufferReal.Resize(N);
    FBufferImage.Resize(N);
    for i:=0 to N-1 do begin
      FBufferReal.Values[i] := aX[i];
      if M>0 then      FBufferImage.Values[i] := aY[i]
      else             FBufferImage.Values[i] := 0;
    end;
end;

procedure TFourier.SetSourceData(aX: TFloatArray);
Var
    i,N: Integer;
begin
    N := Length(aX);
    FBufferReal.Resize(N);
    FBufferImage.Resize(N);
    for i:=0 to N-1 do begin
      FBufferReal.Values[i] := aX[i];
      FBufferImage.Values[i] := 0;
    end;
end;

end.