ahofft.pas

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

    end;
  end;

  procedure CplxSqrtDbl(zr, zi : TFFTDouble; var rr, ri : TFFTDouble);
    var
      tr : TFFTDouble;
  begin
    if (zr = 0) and (zi = 0) then begin
      rr := 0;
      ri := 0;
    end
    else begin
      tr := sqrt((abs(zr) + sqrt(zr*zr + zi*zi)) / 2);
      if (zr >= 0) then begin
        rr := tr;
        ri := zi / (2*tr);
      end
      else if (zi >= 0) then begin
        rr := zi / (2*tr);
        ri := tr;
      end
      else begin
        rr := -zi / (2*tr);
        ri := -tr;
      end;
    end;
  end;

  const
    nn = 4096; // readjust frequency
  var
    a,b : PArrayOfFloat; // auxiliary arrays
    s,k,t,r,d,p,lk,t2 : Integer;
    tmp : TFFTFloat;
    vr,vi,qr,qi,zr,zi,tr,ti : TFFTDouble;
    rp,ip : ^TFFTFloat;
begin
  if (m <= 0) then Exit;
  if (n <= 0) then begin
    for s := 0 to m-1 do begin
      realOut^[s] := 0;
      imagOut^[s] := 0;
    end;
    Exit;
  end;

  // w := Sqrt(w);
  CplxSqrtDbl(wr,wi, wr,wi);

  // determine size of the auxiliary arrays
  k := 8;
  lk := 3;
  while (k < n) do begin
    k := k*2;
    lk := lk + 1;
  end;
  d := k-(n-1);
  tr := (((m+d-1)/d)*2.0 + 1.0) * k*lk;
  r := 2*k-(n-1);
  zr := (((m+r-1)/r)*2.0 + 1.0) * 2.0*k*(lk+1.0);
  while (zr < tr) do begin
    k := k*2;
    lk := lk + 1;
    d := r;
    tr := zr;
    r := 2*k-(n-1);
    zr := (((m+r-1)/r)*2.0 + 1.0) * 2.0*k*(lk+1.0);
  end;

  // get temp space
  s := 2*k;
  GetMem(a, 2*s*SizeOf(TFFTFloat));
  b := @a^[s];
  try
    // calculate array a
    { The powers of w are calculated iteratvely by the following recurrences:
      w^((s+1)^2) = w^(s^2) * w^(1+2*s);
      w^(1+2*(s+1)) = w^(1+2*s) * w^2
      z[0] := w^(s^2);
      q[0] := w^(1+2*s)
      z[n+1] := z[n] * q[n];
      q[n+1] := q[n] * w^2;  }

    zr := 1;
    zi := 0;
    vr := wr;
    vi := wi;
    // q := w^2;
    CplxPowerIntDbl(wr,wi, 2, qr,qi);
    rp := @realIn[0];
    ip := @imagIn[0];
    for t := 0 to n-1 do begin
      // a[t] := x[t] * (w ^ (t^2));

      // a[t] := in[t] * z
      tr := rp^;
      ti := ip^;
      t2 := t+t;
      a^[t2]   := tr*zr - ti*zi;
      a^[t2+1] := tr*zi + ti*zr;

      // z := z * v;
      tr := zr*vr - zi*vi;
      zi := zr*vi + zi*vr;
      zr := tr;
      // v := v * q;
      tr := vr*qr - vi*qi;
      vi := vr*qi + vi*qr;
      vr := tr;
      if (t mod nn = 0) then begin
        // z := z / Abs(z);
        tr := 1.0 / sqrt(zr*zr + zi*zi);
        zr := zr * tr;
        zi := zi * tr;
        // v := v / Abs(v);
        tr := 1.0 / sqrt(vr*vr + vi*vi);
        vr := vr * tr;
        vi := vi * tr;
      end;
      Inc(rp, realStepIn);
      Inc(ip, imagStepIn);
    end;

    if (ar <> 1.0) or (ai <> 0.0) then begin
      for t := 0 to n-1 do begin
        // a[t] := a[t] * a;
        t2 := t+t;
        tr       := a^[t2]*ar - a^[t2+1]*ai;
        a^[t2+1] := a^[t2]*ai + a^[t2+1]*ar;
        a^[t2]   := tr;
      end;
    end;

    if (br <> 1.0) or (bi <> 0.0) then begin
      zr := 1;
      zi := 0;
      for t := 0 to n-1 do begin
        // a[t] := a[t] * b^k;
        t2 := t+t;
        tr       := a^[t2]*zr - a^[t2+1]*zi;
        a^[t2+1] := a^[t2]*zi + a^[t2+1]*zr;
        a^[t2]   := tr;

        // z := z * b;
        tr := zr*br - zi*bi;
        zi := zr*bi + zi*br;
        zr := tr;

        if (t mod nn = 0) then begin
          // z := z / Abs(z);
          tr := 1.0 / sqrt(zr*zr + zi*zi);
          zr := zr * tr;
          zi := zi * tr;
        end;
      end;
    end;

    for t := n to k-1 do begin
      t2 := t+t;
      a^[t2]   := 0;
      a^[t2+1] := 0;
    end;
    SandeTukeyFFT(a,2, @a[1],2, k,1);

    // loop for as many spectrum values as needed
    p := 0;
    while (p < m) do begin
      if (p + d > m) then d := m-p;

      // calculate array b

      //  q := w^(-2);
      CplxPowerIntDbl(wr,wi, -2, qr,qi);
      s := -(n-1) + p;
      // z := w ^ -(s^2) = (w^(-s))^s
      // v := w ^ -(1 + 2*s)
      CplxPowerIntDbl(wr,wi, -s, zr,zi);
      CplxPowerIntDbl(zr,zi,  s, zr,zi);
      CplxPowerIntDbl(wr,wi, -(1+2*s), vr,vi);
      for t := 0 to n+d-2 do begin
        // b[t] := w ^ -((t-(n-1)+startIndex)^2);

        t2 := t+t;
        b^[t2]   := zr;
        b^[t2+1] := zi;

        // z := z * v;
        tr := zr*vr - zi*vi;
        zi := zr*vi + zi*vr;
        zr := tr;
        // v := v * q;
        tr := vr*qr - vi*qi;
        vi := vr*qi + vi*qr;
        vr := tr;
        if (t mod nn = 0) then begin
          // z := z / Abs(z);
          tr := 1.0 / sqrt(zr*zr + zi*zi);
          zr := zr * tr;
          zi := zi * tr;
          // v := v / Abs(v);
          tr := 1.0 / sqrt(vr*vr + vi*vi);
          vr := vr * tr;
          vi := vi * tr;
        end;
      end;
      for t := n+d-1 to k-1 do begin
        t2 := t+t;
        b^[t2]   := 0;
        b^[t2+1] := 0;
      end;

      // convolution of a and b
      SandeTukeyFFT(b,2, @b[1],2, k,1);
      for t := 0 to k-1 do begin
        t2 := t+t;
        tmp      := a^[t2]  *b^[t2] - a^[t2+1]*b^[t2+1];
        b^[t2+1] := a^[t2+1]*b^[t2] + a^[t2]  *b^[t2+1];
        b^[t2]   := tmp;
      end;
      CooleyTukeyFFT(b,2, @b[1],2, k,-1);

      // rescale and store result

      // q := w^2;
      CplxPowerIntDbl(wr,wi, 2, qr,qi);
      s := p;
      // z := w ^ (s^2) = (w^s)^s
      // v := w ^ (1 + 2*s)
      CplxPowerIntDbl(wr,wi, s, zr,zi);
      CplxPowerIntDbl(zr,zi, s, zr,zi);
      CplxPowerIntDbl(wr,wi, 1+2*s, vr,vi);
      rp := @realOut[p*realStepOut];
      ip := @imagOut[p*imagStepOut];
      for t := 0 to d-1 do begin
        // x[t] := b[r] * (w ^ -((t+startIndex)^2));

        r := (n-1)+t;
        t2 := r+r;
        rp^ := zr*b^[t2] - zi*b^[t2+1];
        ip^ := zi*b^[t2] + zr*b^[t2+1];

        // z := z * v;
        tr := zr*vr - zi*vi;
        zi := zr*vi + zi*vr;
        zr := tr;
        // v := v * q;
        tr := vr*qr - vi*qi;
        vi := vr*qi + vi*qr;
        vr := tr;
        if (t mod nn = 0) then begin
          // z := z / Abs(z);
          tr := 1.0 / sqrt(zr*zr + zi*zi);
          zr := zr * tr;
          zi := zi * tr;
          // v := v / Abs(v);
          tr := 1.0 / sqrt(vr*vr + vi*vi);
          vr := vr * tr;
          vi := vi * tr;
        end;
        Inc(rp, realStepOut);
        Inc(ip, imagStepOut);
      end;
      p := p + d;
    end;

    if (cr <> 1.0) or (ci <> 0.0) then begin
      zr := 1.0 / k;
      zi := 0;
      rp := @realOut[0];
      ip := @imagOut[0];
      for s := 0 to m-1 do begin
        // y[s] := y[s] * (1/k) * c^s;
        tr  := rp^ * zr - ip^ * zi;
        ip^ := rp^ * zi + ip^ * zr;
        rp^ := tr;

        // z := z * c;
        tr := zr*cr - zi*ci;
        zi := zr*ci + zi*cr;
        zr := tr;

        if (s mod nn = 0) then begin
          // z := z / Abs(z) * 1/k;
          tr := 1.0 / (k * sqrt(zr*zr + zi*zi));
          zr := zr * tr;
          zi := zi * tr;
        end;
        Inc(rp, realStepOut);
        Inc(ip, imagStepOut);
      end;
    end
    else begin
      tmp := 1.0 / k;
      rp := @realOut[0];
      ip := @imagOut[0];
      for s := 0 to m-1 do begin
        rp^ := rp^ * tmp;
        ip^ := ip^ * tmp;
        Inc(rp, realStepOut);
        Inc(ip, imagStepOut);
      end;
    end;
  finally
    FreeMem(a);
  end;
end;



procedure InPlaceFFTCplx(real : PArrayOfFloat; realStep : Integer;
                         imag : PArrayOfFloat; imagStep : Integer;
                         n : Integer; inverse : Boolean);
  var
    i : Integer;
    wr,wi, theta : TFFTDouble;
    r : TFFTFloat;
    rp,ip : ^TFFTFloat;
begin
  if (n <= 1) then Exit;
  if ((n and (n-1)) = 0) then begin
    // n is a power of 2
    Power2LengthFFT(real, realStep, imag, imagStep,  n, inverse);
  end
  else begin
    theta := -2*Pi/n;
    if (inverse) then theta := -theta;
    wr := cos(theta);
    wi := sin(theta);
    Chirp1(real, realStep, imag, imagStep, n,
              wr, wi, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,
              real, realStep, imag, imagStep, n);
    if (inverse) then begin
      r := 1/n;
      rp := @real^[0];
      ip := @imag^[0];
      for i := 0 to n-1 do begin
        rp^ := rp^ * r;
        ip^ := ip^ * r;
        Inc(rp, realStep);
        Inc(ip, imagStep);
      end;
    end;
  end;
end;



procedure FFTReal(signal : PArrayOfFloat; signalStep : Integer;
                  real : PArrayOfFloat; realStep : Integer;
                  imag : PArrayOfFloat; imagStep : Integer;
                  n : Integer);
  var
    i, n2 : Integer;
    rp,ip,sp,xp : ^TFFTFloat;
begin
  if (n <= 0) then Exit;
  if Odd(n) then begin
    // Do a complex FFT with imaginary parts set to 0
    rp := @real[0];
    ip := @imag[0];
    sp := @signal[0];
    for i := 0 to n-1 do begin
      rp^ := sp^;
      ip^ := 0;
      Inc(rp, realStep);
      Inc(ip, imagStep);
      Inc(sp, signalStep);
    end;
    InPlaceFFTCplx(real, realStep, imag, imagStep, n, false);
  end
  else begin
    // Interleave odd and even samples of the signal as real
    // and imaginary parts of a complex signal with half of the length
    n2 := n div 2;
    rp := @real[0];
    ip := @imag[0];
    sp := @signal[0];
    for i := 0 to n2-1 do begin
      rp^ := sp^;
      Inc(rp, realStep);
      Inc(sp, signalStep);
      ip^ := sp^;
      Inc(ip, imagStep);
      Inc(sp, signalStep);
    end;
    InPlaceFFTCplx(real, realStep, imag, imagStep, n2, false);
    // reconstruct the spectrum of the original real signal
    RealFFTInterleave(real, realStep, imag, imagStep, n2, 1);
    // mirror the symmetric spectrum
    rp := @real[1*realStep];
    sp := @real[(n-1)*realStep];
    ip := @imag[1*imagStep];
    xp := @imag[(n-1)*imagStep];
    for i := 1 to n2-1 do begin
      sp^ :=  rp^;
      Inc(rp, realStep);
      Dec(sp, realStep);
      xp^ :=  -ip^;
      Inc(ip, imagStep);
      Dec(xp, imagStep);
    end;

    // separate the values combined by RealFFTInterleave
    real^[n2*realStep] := imag^[0];
    imag^[n2*imagStep] := 0;
    imag^[0] := 0;
  end;
end;



procedure IFFTReal(real : PArrayOfFloat; realStep : Integer;
                   imag : PArrayOfFloat; imagStep : Integer;
                   signal : PArrayOfFloat; signalStep : Integer;
                   n : Integer);
  var
    i,n2 : Integer;
    tmp : PArrayOfFloat;
    rp,ip,sp : ^TFFTFloat;
begin
  if (n <= 0) then Exit;
  if Odd(n) then begin
    // Do a complex FFT and take out the real parts
    GetMem(tmp, n*sizeof(TFFTFloat));
    try
      rp := @real[0];
      ip := @imag[0];
      sp := @signal[0];
      for i := 0 to n-1 do begin
        sp^ := rp^;
        tmp^[i] := ip^;
        Inc(rp, realStep);
        Inc(ip, imagStep);
        Inc(sp, signalStep);
      end;
      InPlaceFFTCplx(signal,signalStep, tmp,1, n, true);
    finally
      FreeMem(tmp);
    end;
  end
  else begin
    n2 := n div 2;
    rp := @real[0];
    ip := @imag[0];
    sp := @signal[0];
    for i := 0 to n2-1 do begin
      sp^ := rp^;
      Inc(rp, realStep);
      Inc(sp, signalStep);
      sp^ := ip^;
      Inc(ip, imagStep);
      Inc(sp, signalStep);
    end;
    signal^[signalStep] := real[n2*realStep];
    RealFFTInterleave(signal,              2*signalStep,
                      @signal[signalStep], 2*signalStep, n2, -1);
    InPlaceFFTCplx   (signal,              2*signalStep,
                      @signal[signalStep], 2*signalStep, n2, true);
  end;
end;



procedure Chirp(realIn : PArrayOfFloat; realStepIn : Integer;
                imagIn : PArrayOfFloat; imagStepIn : Integer;
                n : Integer;
                wr,wi, ar,ai, br,bi, cr,ci : Extended;
                realOut : PArrayOfFloat; realStepOut : Integer;
                imagOut : PArrayOfFloat; imagStepOut : Integer;
                m : Integer);
  const
    threshold = 1e-15;
begin
  if (Abs(wr*wr + wi*wi - 1) < threshold)
       and (Abs(ar*ar + ai*ai - 1) < threshold)
       and (Abs(br*br + bi*bi - 1) < threshold)
       and (Abs(cr*cr + ci*ci - 1) < threshold)  then begin
    // Consider the length of the complex numbers as 1 and
    // use Chirp1 (faster and better numeric stability)
    Chirp1(realIn, realStepIn, imagIn, imagStepIn, n,
           wr,wi, ar,ai, br,bi, cr,ci,
           realOut, realStepOut, imagOut, imagStepOut, m);
  end
  else begin
    // Use ChirpZ for all the other complex numbers
    ChirpZ(realIn, realStepIn, imagIn, imagStepIn, n,
           wr,wi, ar,ai, br,bi, cr,ci,
           realOut, realStepOut, imagOut, imagStepOut, m);
  end;
end;


//=== Basic FFT algorithms ==================================================