ahofft.pas

CT图象重建的DELPHI源代码,转载自一个学长的博客

unit AHoFFT;

{$BOOLEVAL OFF }
{$EXTENDEDSYNTAX ON }
{$IOCHECKS ON }
{$LONGSTRINGS ON }
{$OPENSTRINGS ON }
{$TYPEDADDRESS OFF }
{$VARSTRINGCHECKS ON }
{$WRITEABLECONST OFF }

// {$DEFINE DEMOVERSION}

INTERFACE

{$IFDEF DEMOVERSION}
uses
  Dialogs;
{$ENDIF}

{$IFNDEF DOUBLEFFTFLOAT }
  type
    TFFTFloat = Single;
    TFFTDouble = Double;  // for intermediate calculations
{$ELSE}
  type
    TFFTFloat = Double;
    TFFTDouble = Extended;  // for intermediate calculations
{$ENDIF}

type
  TFFTComplex = record
    real : TFFTFloat;
    imag : TFFTFloat;
  end;

function Cplx(const re,im : TFFTFloat) : TFFTComplex;
function CplxPolar(const r,phi : TFFTFloat) : TFFTComplex;

function CplxArg(const z : TFFTComplex) : TFFTFloat;
function CplxConj(const z : TFFTComplex) : TFFTComplex;

function CplxAdd(const a,b : TFFTComplex) : TFFTComplex;
function CplxSub(const a,b : TFFTComplex) : TFFTComplex;
function CplxMul(const a,b : TFFTComplex) : TFFTComplex;
function CplxQuo(const a,b : TFFTComplex) : TFFTComplex;

function CplxSqrt(const x : TFFTComplex) : TFFTComplex;
function CplxPowerInt(const x : TFFTComplex; exponent : Integer) : TFFTComplex;

procedure FFT(const signal   : Array of TFFTFloat;
              var   spectrum : Array of TFFTComplex);
// Real signal FFT

procedure IFFT(const spectrum : Array of TFFTComplex;
               var   signal   : Array of TFFTFloat);
// Real signal inverse FFT

procedure FFTCplx(const signal   : Array of TFFTComplex;
                  var   spectrum : Array of TFFTComplex);
// Complex signal FFT

procedure IFFTCplx(const spectrum : Array of TFFTComplex;
                   var   signal   : Array of TFFTComplex);

procedure IFFFT(const spectrum : Array of TFFTComplex;
                fStart, fStep : TFFTDouble;
                var signal : Array of TFFTFloat;
                tStart, tStep : TFFTDouble;
                scaleFactor : TFFTFloat);
// Real signal inverse FFFT


procedure FFFTCplx(const signal : Array of TFFTComplex;
                   tStart, tStep : TFFTDouble;
                   var spectrum : Array of TFFTComplex;
                   fStart, fStep : TFFTDouble);
// Complex signal FFFT

procedure IFFFTCplx(const spectrum : Array of TFFTComplex;
                    fStart, fStep : TFFTDouble;
                    var signal : Array of TFFTComplex;
                    tStart, tStep : TFFTDouble;
                    scaleFactor : TFFTFloat);
// Complex signal inverse FFT

procedure ChirpTransform(const signal : Array of TFFTComplex;
                         wr,wi, ar,ai, br,bi, cr,ci : Extended;
                         var spectrum : Array of TFFTComplex);

//--- General weighting windows

type
  TFFTWindowFunc = function(x : TFFTFloat) : TFFTFloat;

function WTriangle(x : TFFTFloat) : TFFTFloat;
// Also called Parzen or Bartlett window.
// Goes linearly from 0 to 1 and back to 0.

function WWelch(x : TFFTFloat) : TFFTFloat;
// This window is formed by a parabola which goes from 0 to 1 and back to 0. 

type
  TFFTWindowInfo = record
                     intw  : TFFTFloat; // integral of window
                     intw2 : TFFTFloat; // integral of window squared
                   end;

function ApplyWindow(var data : array of TFFTFloat;
                     winFunc  : TFFTWindowFunc) : TFFTWindowInfo;

function ApplySumCosineWindow(var data : array of TFFTFloat;
                              const coeff : array of Double) : TFFTWindowInfo;
// Apply the sum cosine window function defined by 'coeff' to data.
// Replaces the values in data.


function WSumCosine(const coeff : array of TFFTDouble; x : TFFTFloat) : TFFTFloat;

const  // info: first sidelobe level [dB]; sidelobe asymptotic behaviour [db/Oct]
  wCoeffHann     : array[0..1] of Double = (0.5, 0.5);  // 31.5dB; 18dB/Oct
  wCoeffHamming  : array[0..1] of Double = (0.54, 0.46);  // 43dB; 6dB/Oct
  wCoeffBlackman : array[0..2] of Double = (0.42, 0.5, 0.08);  // 58.11dB; 18dB/Oct
  wCoeffNuttall  : array[0..3] of Double = (0.355768, 0.487396,
                                      0.144232, 0.012604); // 93.3dB; 18dB/Oct

  wCoeffFlatTop2 : array[0..2] of Double = (0.28235,
                                         0.52105, 0.19659); // 43.2dB; 6dB/Oct
  wCoeffFlatTop3 : array[0..3] of Double = (0.241096, 0.460841,
                                        0.255381, 0.041872); // 66.8dB; 6dB/Oct
  wCoeffFlatTop4 : array[0..4] of Double = (0.209671, 0.407331,
                             0.281225, 0.092669, 0.009104); // 90.5dB; 18dB/Oct

  wCoeffTaylor60 : array[0..11] of Double = (0.46979160648, 0.4897086008,
                       0.040543436437, 0.00046498453306, -0.00089292984937,
                       0.00063096192603, -0.00039727029523, 0.00023660292287,
                       -0.00013197019847, 6.6200331699e-05, -2.7248084479e-05,
                       7.0249981132e-06); // 60dB; 6dB/Oct

  wCoeffTaylor80 : array[0..15] of Double = (0.40987520229, 0.49726668004,
                       0.091629538934, 0.0011747429958, 2.8285040313e-05,
                       6.6896112991e-05, -8.0865998033e-05, 7.0758365573e-05,
                       -5.6300382963e-05, 4.2878214708e-05, -3.1677796476e-05,
                       2.2686692688e-05, -1.5592228864e-05, 1.005086136e-05,
                       -5.7617932293e-06, 2.4786542281e-06); // 80dB; 6dB/Oct


procedure TaylorWinCoeff(sideLobedB : Single; var coeff : array of Double);

//=== Auxiliary procedures ==================================================

const // multiply with these to convert from rad to deg and back
  radToDeg = 180/pi;
  degToRad = pi/180;


function NextPowerOf2(n : Cardinal) : Cardinal;
// if 'n' is a power of 2 return n, otherwise return
// the next greater power of 2. 


procedure RealToCplx(const realSignal    : Array of TFFTFloat;
                     var   complexSignal : Array of TFFTComplex);
// Convert a real signal into a complex one.


procedure RealToCplx2(const realSignal    : Array of TFFTFloat;
                      const imagSignal    : Array of TFFTFloat;
                      var   complexSignal : Array of TFFTComplex);
// Convert two real signals into one complex one.


procedure RealOfCplx(const complexSignal : Array of TFFTComplex;
                     var   realSignal    : Array of TFFTFloat);
// Extract the real part of a complex signal.

procedure ImagOfCplx(const complexSignal : Array of TFFTComplex;
                     var   imagSignal    : Array of TFFTFloat);

IMPLEMENTATION //************************************************************

//=== Complex numbers =======================================================

function Cplx(const re,im : TFFTFloat) : TFFTComplex;
begin
  Result.real := re;
  Result.imag := im;
end;



function CplxPolar(const r,phi : TFFTFloat) : TFFTComplex;
begin
  Result.real := r * Cos(phi);
  Result.imag := r * Sin(phi);
end;




function CplxArg(const z : TFFTComplex) : TFFTFloat;
begin
  // returns -pi .. +pi
  if z.real = 0 then begin
    if z.imag = 0 then
      Result := 0
    else if z.imag < 0 then
      Result := -pi/2
    else begin
      Result := pi/2;
    end;
  end
  else if z.real < 0 then begin
    if z.imag > 0 then
      Result := ArcTan(z.imag/z.real) + pi
    else begin
      Result := ArcTan(z.imag/z.real) - pi
    end;
  end
  else begin
    Result := ArcTan(z.imag/z.real);
  end;
end;



function CplxConj(const z : TFFTComplex) : TFFTComplex;
begin
  Result.real :=  z.real;
  Result.imag := -z.imag;
end;



function CplxAdd(const a,b : TFFTComplex) : TFFTComplex;
begin
  Result.real := a.real + b.real;
  Result.imag := a.imag + b.imag;
end;



function CplxSub(const a,b : TFFTComplex) : TFFTComplex;
begin
  Result.real := a.real - b.real;
  Result.imag := a.imag - b.imag;
end;



function CplxMul(const a,b : TFFTComplex) : TFFTComplex;
begin
  Result.real := a.real*b.real - a.imag*b.imag;
  Result.imag := a.real*b.imag + a.imag*b.real;
end;



function CplxQuo(const a,b : TFFTComplex) : TFFTComplex;
  var
    d : Single;
begin
  d := 1 / (b.real*b.real + b.imag*b.imag);
  Result.real := (a.real*b.real + a.imag*b.imag) * d;
  Result.imag := (a.imag*b.real - a.real*b.imag) * d;
end;



function CplxPowerInt(const x : TFFTComplex; exponent : Integer) : TFFTComplex;
  var
    e : Integer;
    xr,xi,xt,rr,ri : TFFTFloat;
begin
  e := exponent;
  if (e < 0) then e := -e;
  xr := x.real;
  xi := x.imag;
  rr := 1;
  ri := 0;
  while (e > 0) do begin
    while ((e and 1) = 0) do begin
      e := e div 2;
      // x := x*x
      xt := (xr+xi)*(xr-xi);
      xi := 2 * xi * xr;
      xr := xt;
    end;
    e := e - 1;
    // r := r*x
    xt := rr*xr - ri*xi;
    ri := ri*xr + rr*xi;
    rr := xt;
  end;
  if (exponent < 0) then begin
    // r := 1/r
    xt := 1 / (Sqr(rr) + Sqr(ri));
    rr :=  rr * xt;
    ri := -ri * xt;
  end;
  Result.real := rr;
  Result.imag := ri;
end;



function CplxSqrt(const x : TFFTComplex) : TFFTComplex;
  var
    tr : TFFTFloat;
begin
  if (x.real = 0) and (x.imag = 0) then begin
    Result.real := 0;
    Result.imag := 0;
  end
  else begin
    tr := sqrt((abs(x.real) + sqrt(Sqr(x.real)+Sqr(x.imag))) / 2);
    if (x.real >= 0) then begin
      Result.real := tr;
      Result.imag := x.imag / (2*tr);
    end
    else if (x.imag >= 0) then begin
      Result.real := x.imag / (2*tr);
      Result.imag := tr;
    end
    else begin
      Result.real := -x.imag / (2*tr);
      Result.imag := -tr;
    end;
  end;
end;


//=== Low-level routines ====================================================

type
  TArray_OfFloat = Array[0..High(Longint) div SizeOf(TFFTFloat)-1] of TFFTFloat;
  PArrayOfFloat = ^TArray_OfFloat;


procedure BitReversal(real : PArrayOfFloat; realStep : Integer;
                      imag : PArrayOfFloat; imagStep : Integer;
                      n : Integer);
  var
    i,j,m : Integer;
    pr,pi,pj : ^TFFTFloat;
    swap : TFFTFloat;
begin
  j := 0;
  pr := @real[0];
  pi := @imag[0];
  for i := 0 to n-1 do begin
    // loop invariant: j is i in reversed bit order

    // swap places i and j if not already done so before
    if (j > i) then begin
      // swap real data
      pj := @real^[realStep*j];
      swap := pr^;
      pr^ := pj^;
      pj^ := swap;
      // swap imag data
      pj := @imag^[imagStep*j];
      swap := pi^;
      pi^ := pj^;
      pj^ := swap;
    end;

    // increment j reversely
    m := n div 2;
    while (m >= 2) and (j >= m) do begin
      j := j - m;
      m := m div 2;
    end;
    j := j + m;

    // increment pointers
    Inc(pr, realStep);
    Inc(pi, imagStep);
  end;
end;



procedure CooleyTukeyFFT(real : PArrayOfFloat; realStep : Integer;
                         imag : PArrayOfFloat; imagStep : Integer;
                         n : Integer; isign : Integer);
  var
    n1,n2,i,j : Integer;
    tempr, tempi : TFFTFloat;
    nr,ni : Integer;
    rpi,ipi,rpk,ipk,rp,ip : ^TFFTFloat;
    // use double precision for the trigonometric recurrences
    wpr, wpi, wr, wi, wt, theta : TFFTDouble;
begin
  nr := realStep;
  ni := imagStep;
  n1 := 1;
  while (n > n1) do begin
    n2 := n1;
    n1 := n1 * 2;

    // initialize for trigonometric recurrence
    theta := Pi/(-isign*n2);
    wpr := -2.0 * Sqr(Sin(0.5*theta)); // = cos(theta)-1
    wpi := sin(theta);
    wr := 1.0;
    wi := 0.0;

    rp := @real^[0];
    ip := @imag^[0];
    for j := 0 to n2-1 do begin
      i := j;
      rpi := rp;
      ipi := ip;
      while i < n do begin
        // the code below does the following:  ( with complex numbers )
        //   k := i + n2;
        //   temp := w * data[k];
        //   data[k] := data[i] - temp;
        //   data[i] := data[i] + temp;

        rpk := rpi;
        Inc(rpk,nr);
        ipk := ipi;
        Inc(ipk,ni);

        tempr := wr*rpk^ - wi*ipk^;
        tempi := wr*ipk^ + wi*rpk^;
        rpk^  := rpi^ - tempr;
        ipk^  := ipi^ - tempi;
        rpi^  := rpi^ + tempr;
        ipi^  := ipi^ + tempi;

        rpi := rpk;
        Inc(rpi,nr);
        ipi := ipk;
        Inc(ipi,ni);

        i := i + n1;
      end;

      // trigonometric recurrence
      wt := wr;
      wr := wr*wpr - wi*wpi + wr;
      wi := wi*wpr + wt*wpi + wi;

      Inc(rp, realStep);
      Inc(ip, imagStep);
    end;

    nr := nr * 2;
    ni := ni * 2;
  end;
end;



procedure SandeTukeyFFT(real : PArrayOfFloat; realStep : Integer;
                        imag : PArrayOfFloat; imagStep : Integer;
                        n : Integer; isign : Integer);
  var
    n1,n2,i,j : Integer;
    tempr, tempi : TFFTFloat;
    nr,ni : Integer;
    rpi,ipi,rpk,ipk,rp,ip : ^TFFTFloat;
    // use double precision for the trigonometric recurrences
    wpr, wpi, wr, wi, wt, theta : TFFTDouble;
begin
  nr := n*realStep;
  ni := n*imagStep;
  n2 := n;
  while (n2 > 1) do begin
    n1 := n2;
    n2 := n2 div 2;
    nr := nr div 2;
    ni := ni div 2;

    // initialize for trigonometric recurrence
    theta := Pi/(-isign*n2);
    wpr := -2.0 * Sqr(Sin(0.5*theta)); // = cos(theta)-1
    wpi := sin(theta);
    wr := 1.0;
    wi := 0.0;

    rp := @real^[0];
    ip := @imag^[0];
    for j := 0 to n2-1 do begin
      i := j;
      rpi := rp;
      ipi := ip;
      while i < n do begin
        rpk := rpi;
        Inc(rpk,nr);
        ipk := ipi;