ahofft.pas

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

        Inc(ipk,ni);

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

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

        i := i + n1;
      end;

      wt := wr;
      wr := wr*wpr - wi*wpi + wr;
      wi := wi*wpr + wt*wpi + wi;

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




procedure Power2LengthFFT(real : PArrayOfFloat; realStep : Integer;
                          imag : PArrayOfFloat; imagStep : Integer;
                          n : Integer; inverse : Boolean);
  var
    sign : Integer;
    r : TFFTFloat;
    i : Integer;
    rp,ip : ^TFFTFloat;
begin
  if (n < 2) then Exit;
  sign := 1;
  if (inverse) then sign := -1;

  SandeTukeyFFT(real, realStep, imag, imagStep, n, sign);
  BitReversal(real, realStep, imag, imagStep, n);
  if (inverse) then begin
    r := 1.0/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;



procedure RealFFTInterleave(real : PArrayOfFloat; realStep : Integer;
                            imag : PArrayOfFloat; imagStep : Integer;
                            n : Integer; isign : Integer);
  var
    wr,wi,wpr,wpi,wt,theta : TFFTDouble;
    i : Integer;
    c1,c2,h1r,h1i,h2r,h2i,tr,ti : TFFTFloat;
    rpi,ipi, rpj,ipj : ^TFFTFloat;
begin
  theta := Pi/(-n*isign);
  c1 := 0.5;
  c2 := isign * 0.5;

  // wp := exp(j*theta) - 1
  wpr := -2.0*sqr(sin(0.5*theta)); // = cos(theta) - 1
  wpi := sin(theta);
  // w := (wp+1) ^ i
  wr := 1.0 + wpr;
  wi := wpi;
  rpi := @real^[1*realStep];
  ipi := @imag^[1*imagStep];
  rpj := @real^[(n-1)*realStep];
  ipj := @imag^[(n-1)*imagStep];
  for i := 1 to n div 2 do begin
    // j := n-i;
    // h1r := 1/2 * (data[i] + Conj(data[j]))
    h1r :=  c1 * (rpi^ + rpj^);
    h1i :=  c1 * (ipi^ - ipj^);
    // h2r := +/- 1/2 * (-j) * (data[i] - Conj(data[j]))
    h2r :=  c2 * (ipi^ + ipj^);
    h2i := -c2 * (rpi^ - rpj^);
    // t := w * h2r
    tr := wr*h2r - wi*h2i;
    ti := wr*h2i + wi*h2r;
    // data[i] := h1 + t
    rpi^ :=  h1r + tr;
    ipi^ :=  h1i + ti;
    // data[j] := Conj(h1 - t);
    rpj^ :=  h1r - tr;
    ipj^ := -h1i + ti;
    // w := w * (wp+1)
    wt := wr;
    wr := wr*wpr - wi*wpi + wr;
    wi := wi*wpr + wt*wpi + wi;

    Inc(rpi, realStep);
    Inc(ipi, imagStep);
    Dec(rpj, realStep);
    Dec(ipj, imagStep);
  end;

  // special case data[0] and data[n/2] merged into data[0]
  if (isign > 0) then begin
    h1r := real^[0];
    real^[0] := h1r + imag^[0];
    imag^[0] := h1r - imag^[0];
  end
  else begin
    h1r := real^[0];
    real^[0] := c1 * (h1r + imag^[0]);
    imag^[0] := c1 * (h1r - imag^[0]);
  end;
end;



procedure ChirpZ(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);

  procedure CplxPowerIntExt(zr, zi : Extended;
                            exponent : Integer;
                            var rr, ri : Extended);
    var
      e : Integer;
      xr,xi,xt : Extended;
  begin
    e := exponent;
    if (e < 0) then e := -e;
    xr := zr;
    xi := zi;
    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;
  end;

  procedure CplxSqrtExt(zr, zi : Extended; var rr, ri : Extended);
    var
      tr : Extended;
  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 : Extended;
    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);
  CplxSqrtExt(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;  }

    // q := w^2;
    CplxPowerIntExt(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));

      if (t mod nn = 0) then begin
        // z := w ^ (t^2) = (w^t)^t
        // v := w ^ (1 + 2*t)
        CplxPowerIntExt(wr,wi, t, zr,zi);
        CplxPowerIntExt(zr,zi, t, zr,zi);
        CplxPowerIntExt(wr,wi, 1+2*t, vr,vi);
      end;

      // 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;

      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
      for t := 0 to n-1 do begin
        if (t mod nn = 0) then begin
          // z := b^t
          CplxPowerIntExt(br,bi, t, zr,zi);
        end;

        // a[t] := a[t] * b^t;
        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;
      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);
      CplxPowerIntExt(wr,wi, -2, qr,qi);
      for t := 0 to n+d-2 do begin
        // b[t] := w ^ -((t-(n-1)+startIndex)^2);

        if (t mod nn = 0) then begin
          s := t - (n-1) + p;
          // z := w ^ -(s^2) = (w^(-s))^s
          // v := w ^ -(1 + 2*s)
          CplxPowerIntExt(wr,wi, -s, zr,zi);
          CplxPowerIntExt(zr,zi,  s, zr,zi);
          CplxPowerIntExt(wr,wi, -(1+2*s), vr,vi);
        end;

        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;
      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;
      CplxPowerIntExt(wr,wi, 2, qr,qi);
      rp := @realOut[p*realStepOut];
      ip := @imagOut[p*imagStepOut];
      for t := 0 to d-1 do begin
        // x[t] := b[r] * (w ^ -((t+startIndex)^2));

        if (t mod nn = 0) then begin
          s := t + p;
          // z := w ^ (s^2) = (w^s)^s
          // v := w ^ (1 + 2*s)
          CplxPowerIntExt(wr,wi, s, zr,zi);
          CplxPowerIntExt(zr,zi, s, zr,zi);
          CplxPowerIntExt(wr,wi, 1+2*s, vr,vi);
        end;

        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;

        Inc(rp, realStepOut);
        Inc(ip, imagStepOut);
      end;
      p := p + d;
    end;

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

        // 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;

        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 Chirp1(realIn : PArrayOfFloat; realStepIn : Integer;
                 imagIn : PArrayOfFloat; imagStepIn : Integer;
                 n : Integer;
                 wr,wi, ar,ai, br,bi, cr,ci : TFFTDouble;
                 realOut : PArrayOfFloat; realStepOut : Integer;
                 imagOut : PArrayOfFloat; imagStepOut : Integer;
                 m : Integer);

  procedure CplxPowerIntDbl(zr, zi : TFFTDouble;
                            exponent : Integer;
                            var rr, ri : TFFTDouble);
    var
      e : Integer;
      xr,xi,xt : TFFTDouble;
  begin
    e := exponent;
    if (e < 0) then e := -e;
    xr := zr;
    xi := zi;
    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;