ahofft.pas

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

procedure FFTCplx(const signal   : Array of TFFTComplex;
                  var   spectrum : Array of TFFTComplex);
  var
    i, n : Integer;
begin
  n := Length(spectrum);
  if n = 0 then Exit;
  Assert(Length(spectrum) >= Length(signal), 'FFTCplx: spectrum array too short');
  for i := 0 to n-1 do spectrum[i] := signal[i];
  InPlaceFFTCplx(Pointer(@spectrum[0].real),2,
                 Pointer(@spectrum[0].imag),2,  n, false);
end;



procedure IFFTCplx(const spectrum : Array of TFFTComplex;
                   var   signal   : Array of TFFTComplex);
  var
    i, n : Integer;
begin
  n := Length(spectrum);
  if n = 0 then Exit;
  Assert(Length(signal) >= Length(spectrum), 'IFFTCplx: signal array too short');
  for i := 0 to n-1 do signal[i] := spectrum[i];
  InPlaceFFTCplx(Pointer(@signal[0].real),2,
                 Pointer(@signal[0].imag),2,  n, true);
end;



procedure FFT(const signal   : Array of TFFTFloat;
              var   spectrum : Array of TFFTComplex);
  var
    n : Integer;
begin
  n := Length(signal);
  if n = 0 then Exit;
  Assert(Length(spectrum) >= n, 'FFT: spectrum array too short');
  FFTReal(Pointer(@signal[0]),1,
          Pointer(@spectrum[0].real),2,
          Pointer(@spectrum[0].imag),2,  n);
end;



procedure IFFT(const spectrum : Array of TFFTComplex;
               var   signal   : Array of TFFTFloat);
  var
    n : Integer;
begin
  n := Length(spectrum);
  if n = 0 then Exit;
  Assert(Length(signal) >= n, 'IFFT: signal array too short');
  IFFTReal(Pointer(@spectrum[0].real),2,
           Pointer(@spectrum[0].imag),2,
           Pointer(@signal[0]),1,  n);
end;



//=== Fractional FFT ========================================================


procedure FFFTCplx(const signal : Array of TFFTComplex;
                   tStart, tStep : TFFTDouble;
                   var spectrum : Array of TFFTComplex;
                   fStart, fStep : TFFTDouble);
  var
    n,m,i : Integer;
    wr,wi, theta : Extended;
    ar,ai, br,bi, cr,ci : Extended;
begin
  n := Length(signal);
  m := Length(spectrum);
  if m = 0 then Exit;
  if n = 0 then begin
    for i := 0 to m-1 do begin
      spectrum[i].real := 0;
      spectrum[i].imag := 0;
    end;
  end
  else begin
    theta := -2*pi*fStep*tStep;
    wr := cos(theta);
    wi := sin(theta);

    theta := -2*pi*fStart*tStart;
    ar := cos(theta);
    ai := sin(theta);

    theta := -2*pi*fStart*tStep;
    br := cos(theta);
    bi := sin(theta);

    theta := -2*pi*fStep*tStart;
    cr := cos(theta);
    ci := sin(theta);

    Chirp1(Pointer(@signal[0].real),2,
           Pointer(@signal[0].imag),2, n,
           wr,wi, ar,ai, br,bi, cr,ci,
           Pointer(@spectrum[0].real),2,
           Pointer(@spectrum[0].imag),2, m);
  end;
end;



procedure IFFFT(const spectrum : Array of TFFTComplex;
                fStart, fStep : TFFTDouble;
                var signal : Array of TFFTFloat;
                tStart, tStep : TFFTDouble;
                scaleFactor : TFFTFloat);
  var
    n,m,i : Integer;
    wr,wi, theta : Extended;
    ar,ai, br,bi, cr,ci : Extended;
    imag : array of TFFTFloat;
begin
  n := Length(signal);
  m := Length(spectrum);
  if n = 0 then Exit;
  if m = 0 then begin
    for i := 0 to m-1 do begin
      signal[i] := 0;
    end;
  end
  else begin
    theta := 2*pi*fStep*tStep;
    wr := cos(theta);
    wi := sin(theta);

    theta := 2*pi*fStart*tStart;
    ar := cos(theta);
    ai := sin(theta);

    theta := 2*pi*fStep*tStart;
    br := cos(theta);
    bi := sin(theta);

    theta := 2*pi*fStart*tStep;
    cr := cos(theta);
    ci := sin(theta);

    SetLength(imag, n);
    Chirp1(Pointer(@spectrum[0].real),2,
           Pointer(@spectrum[0].imag),2, m,
           wr,wi, ar,ai, br,bi, cr,ci,
           Pointer(@signal[0]),1,
           Pointer(@imag[0]),1, n);

    if scaleFactor <> 1 then begin
      for i := 0 to Length(signal)-1 do signal[i] := signal[i] * scaleFactor;
    end;  
  end;
end;



procedure IFFFTCplx(const spectrum : Array of TFFTComplex;
                    fStart, fStep : TFFTDouble;
                    var signal : Array of TFFTComplex;
                    tStart, tStep : TFFTDouble;
                    scaleFactor : TFFTFloat);
  var
    n,m,i : Integer;
    wr,wi, theta : Extended;
    ar,ai, br,bi, cr,ci : Extended;
begin
  n := Length(signal);
  m := Length(spectrum);
  if n = 0 then Exit;
  if m = 0 then begin
    for i := 0 to m-1 do begin
      signal[i].real := 0;
      signal[i].imag := 0;
    end;
  end
  else begin
    theta := 2*pi*fStep*tStep;
    wr := cos(theta);
    wi := sin(theta);

    theta := 2*pi*fStart*tStart;
    ar := cos(theta);
    ai := sin(theta);

    theta := 2*pi*fStep*tStart;
    br := cos(theta);
    bi := sin(theta);

    theta := 2*pi*fStart*tStep;
    cr := cos(theta);
    ci := sin(theta);

    Chirp1(Pointer(@spectrum[0].real),2,
           Pointer(@spectrum[0].imag),2, m,
           wr,wi, ar,ai, br,bi, cr,ci,
           Pointer(@signal[0].real),2,
           Pointer(@signal[0].imag),2, n);

    for i := 0 to Length(signal)-1 do begin
      signal[i].real := signal[i].real * scaleFactor;
      signal[i].imag := signal[i].imag * scaleFactor;
    end;
  end;
end;



procedure ChirpTransform(const signal : Array of TFFTComplex;
                         wr,wi, ar,ai, br,bi, cr,ci : Extended;
                         var spectrum : Array of TFFTComplex);
  var
    n,m : Integer;
begin
  n := Length(signal);
  m := Length(spectrum);
  Chirp(Pointer(@signal[0].real),2,
        Pointer(@signal[0].imag),2, n,
        wr,wi, ar,ai, br,bi, cr,ci,
        Pointer(@spectrum[0].real),2,
        Pointer(@spectrum[0].imag),2, m);
end;


//=== Data Windows ==========================================================


function ApplyWindow(var data : array of TFFTFloat;
                     winFunc  : TFFTWindowFunc) : TFFTWindowInfo;
  var
    n,i : Integer;
    r : TFFTFloat;
    w,ws,wss : TFFTDouble;
begin
  ws := 0;
  wss := 0;
  n := Length(data);
  if n > 0 then begin
    r := 1/n;
    for i := 0 to n-1 do begin
      w := winFunc((i+0.5)*r);
      data[i] := data[i] * w;
      ws := ws + w;
      wss := wss + w*w;
    end;
  end;
  if n = 0 then begin
    Result.intw := 1;
    Result.intw2 := 1;
  end
  else begin
    Result.intw := ws / n;
    Result.intw2 := wss / n;
  end;
end;



function ApplySumCosineWindow(var data : array of TFFTFloat;
                              const coeff : array of Double) : TFFTWindowInfo;
  var
    w : TFFTDouble;
    d,t : TFFTDouble;
    i,j,k,n : Integer;
    s1,c1,s2,c2,s3,c3 : TFFTDouble; // c1:cos(i*d), c2:cos(d)-1, c3:cos(k*i*d);
    coeffSum : TFFTDouble;
begin
  n := Length(data);

  if n > 0 then begin
    d := pi*2/n;
    s2 := sin(d);
    c2 := -2*Sqr(sin(0.5*d)); // = cos(d) - 1

    i := n div 2 - 1;
    if Odd(n) then begin
      coeffSum := coeff[0];
      for k := 1 to High(coeff) do coeffSum := coeffSum + coeff[k];
      data[i+1] := coeffSum * data[i+1];
      j := i+2;
      s1 := s2;
      c1 := c2 + 1;
    end
    else begin
      j := i+1;
      s1 := sin(d/2);
      c1 := cos(d/2);
    end;

    while i >= 0 do begin
      s3 := s1;
      c3 := c1;
      w := coeff[0] + coeff[1]*c3;
      for k := 2 to High(coeff) do begin
        t  := c3*c1 - s3*s1;
        s3 := c3*s1 + s3*c1;
        c3 := t;
        w := w + coeff[k]*c3;
      end;

      data[i] := w * data[i];
      data[j] := w * data[j];

      t  := c1*c2 - s1*s2 + c1;
      s1 := c1*s2 + s1*c2 + s1;
      c1 := t;

      Dec(i);
      Inc(j);
    end;
  end;

  Result.intw := coeff[0];
  Result.intw2 := 0;
  for k := 1 to High(coeff) do Result.intw2 := Result.intw2 + sqr(coeff[k]);
  Result.intw2 := sqr(coeff[0]) + Result.intw2/2;
end;



function WTriangle(x : TFFTFloat) : TFFTFloat;
begin
  Result := 1 - System.Abs(2*x-1);
end;



function WWelch(x : TFFTFloat) : TFFTFloat;
  var
    t : TFFTDouble;
begin
  t := 2*x-1;
  Result := 1 - t*t;
end;



function WSumCosine(const coeff : array of TFFTDouble; x : TFFTFloat) : TFFTFloat;
  var
    j : Integer;
    w,t : TFFTDouble;
    s1,c1,s2,c2 : TFFTDouble;
begin
  t := pi*(2*x-1);
  // the code below does the following:
  //   w := coeff[0];
  //   for j := 1 to High(coeff) do 
  //     w := w + coeff[j]*cos(j*t);

  s1 := sin(t);
  c1 := cos(t);
  s2 := s1;
  c2 := c1;
  w := coeff[0] + coeff[1]*c2;
  for j := 2 to High(coeff) do begin
    t  := c2*c1 - s2*s1;
    s2 := c2*s1 + s2*c1;
    c2 := t;
    w := w + coeff[j]*c2;
  end;
  Result := w;
end;



procedure TaylorWinCoeff(sideLobedB : Single; var coeff : array of Double);
  var
    em,b,c,p,q,s : Double;
    i,j,k : Integer;
begin
  sideLobedB := System.abs(sideLobedB);
  em := exp(-0.11512925465 * sideLobedB);
  k := High(coeff);

  b := ln(1/em + sqrt(sqr(1/em)-1))/Pi;
  c := sqrt(sqr(k+1)/(sqr(b)+sqr(k+0.5)));
  coeff[0] := 1;
  s := 1;
  for j := 1 to k do begin
    p := 1;
    for i := 1 to k do begin
      p := p * (1 - sqr(j/c)/(sqr(b)+sqr(i-0.5)));
    end;
    q := 1;
    for i := 1 to k do begin
      if i <> j then q := q * (1-sqr(j/i));
    end;
    if Odd(j) then
      coeff[j] := p/q
    else begin
      coeff[j] := -p/q;
    end;
    s := s + coeff[j];
  end;
  s := 1/s;
  for j := 0 to k do coeff[j] := coeff[j] * s;
end;


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


function NextPowerOf2(n : Cardinal) : Cardinal;
  const
     limit = Cardinal(High(Cardinal) div 2);
begin
  Assert(n <= limit, 'NextPowerOf2: number too big');
  Result := 1;
  while Result < n do Result := Result * 2;
end;



procedure RealToCplx(const realSignal    : Array of TFFTFloat;
                     var   complexSignal : Array of TFFTComplex);
  var
    n,i : Integer;
begin
  n := Length(realSignal);
  if Length(complexSignal) < n then n := Length(complexSignal);
  for i := 0 to n-1 do begin
    complexSignal[i].real := realSignal[i];
    complexSignal[i].imag := 0;
  end;
  for i := n to Length(complexSignal)-1 do begin
    complexSignal[i].real := 0;
    complexSignal[i].imag := 0;
  end;
end;



procedure RealToCplx2(const realSignal    : Array of TFFTFloat;
                      const imagSignal    : Array of TFFTFloat;
                      var   complexSignal : Array of TFFTComplex);
  var
    n,i : Integer;
begin
  n := Length(realSignal);
  if Length(complexSignal) < n then n := Length(complexSignal);
  for i := 0 to n-1 do complexSignal[i].real := realSignal[i];
  for i := n to Length(complexSignal)-1 do complexSignal[i].real := 0;

  n := Length(imagSignal);
  if Length(complexSignal) < n then n := Length(complexSignal);
  for i := 0 to n-1 do complexSignal[i].imag := imagSignal[i];
  for i := n to Length(complexSignal)-1 do complexSignal[i].imag := 0;
end;



procedure RealOfCplx(const complexSignal : Array of TFFTComplex;
                     var   realSignal    : Array of TFFTFloat);
  var
    n,i : Integer;
begin
  n := Length(complexSignal);
  if Length(realSignal) < n then n := Length(realSignal);
  for i := 0 to n-1 do realSignal[i] := complexSignal[i].real;
  for i := n to Length(realSignal)-1 do realSignal[i] := 0;
end;



procedure ImagOfCplx(const complexSignal : Array of TFFTComplex;
                     var   imagSignal    : Array of TFFTFloat);
  var
    n,i : Integer;
begin
  n := Length(complexSignal);
  if Length(imagSignal) < n then n := Length(imagSignal);
  for i := 0 to n-1 do imagSignal[i] := complexSignal[i].imag;
  for i := n to Length(imagSignal)-1 do imagSignal[i] := 0;
end;

{$IFDEF DEMOVERSION}

initialization
{$ENDIF}

END.