complexs.pas

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

{-------------------------------------------------------------------------------
2003.08.12 根据以下的单元，修改，将所有的缩写改成全拼，定义复数的数组
--------------------------------------------------------------------------------
2003.08.21:
准备开发离散Fourier变换，利用Euler公式将复数的乘法变成加法。有望开发出快速的变换
程序。(没有做)
问题的关键在于复数的处理上，考虑实现为一个类，采用两种形式表示复数：a+ib与re^i()
形式。两种形式有一个转换的过程，还要有一个转换的指示器。
-------------------------------------------------------------------------------}
{ Unit Complexs

  This unit implements complex number arithmetic, including the basic
  operations addition, substraction, multiplication, division,
  magnitude and phase.

  Copyright: Nils Haeck M.Sc. (email: [e-mail])
  For more information visit http://ww.simdesign.nl
  Original date of publication: 10 Mar 2003

  This unit requires these other units:
  - Math:  Delphi mathematics unit
  - Types: Additional mathematical variable types

  ****************************************************************

  The contents of this file are subject to the Mozilla Public
  License Version 1.1 (the "License"); you may not use this file
  except in compliance with the License. You may obtain a copy of
  the License at:
  http://www.mozilla.org/MPL/

  Software distributed under the License is distributed on an
  "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
  implied. See the License for the specific language governing
  rights and limitations under the License.
}
unit Complexs;

interface

uses

  MathTypes, Math;

// Set a complex number
function Complex(Re: TFloat; Im: TFloat): TComplex;

// Add complex numbers (Result = C1 + C2)
function ComplexAdd(const C1, C2: TComplex): TComplex;

// Substract complex numbers (Result = C1 - C2)
function ComplexSub(const C1, C2: TComplex): TComplex;

// Multiply complex numbers (Result = C1 * C2)
function ComplexMultiply(const C1, C2: TComplex): TComplex;

// Scale complex numbers (Result = Scale * C)
function ComplexScale(Scale: TFloat; const C: TComplex): TComplex;

// Get the magnitude of the complex number C
function ComplexMagnitude(const C: TComplex): TFloat;

// Get the phase of the complex number C (in radians, between -pi and pi)
function ComplexPhase(const C: TComplex): TFloat;

// Get the power capacity of the complex number C
function ComplexPower(const C: TComplex): TFloat;

implementation

function Complex(Re: TFloat; Im: TFloat): TComplex;
// Set a complex number
begin
  Result.Re := Re;
  Result.Im := Im;
end;

function ComplexAdd(const C1, C2: TComplex): TComplex;
// Add complex numbers (Result = C1 + C2)
begin
  Result.Re := C1.Re + C2.Re;
  Result.Im := C1.Im + C2.Im;
end;

function ComplexSub(const C1, C2: TComplex): TComplex;
// Substract complex numbers (Result = C1 - C2)
begin
  Result.Re := C1.Re - C2.Re;
  Result.Im := C1.Im - C2.Im;
end;

function ComplexMultiply(const C1, C2: TComplex): TComplex;
// Multiply complex numbers (Result = C1 * C2)
begin
  Result.Re := C1.Re * C2.Re - C1.Im * C2.Im;
  Result.Im := C1.Im * C2.Re + C1.Re * C2.Im;
end;

function ComplexScale(Scale: TFloat; const C: TComplex): TComplex;
// Scale complex numbers (Result = Scale * C)
begin
  Result.Re := Scale * C.Re;
  Result.Im := Scale * C.Im;
end;

function ComplexMagnitude(const C: TComplex): TFloat;
// Get the magnitude of the complex number C
begin
  Result := sqrt(sqr(C.Re) + sqr(C.Im));
end;

function ComplexPhase(const C: TComplex): TFloat;
// Get the phase of the complex number C (in radians, between -pi and pi)
const
  c2Pi  = 2 * pi;
  cPid2 = 0.5 * pi;
begin
  // Both zero
  if (C.Re = 0) and (C.Im = 0) then begin
    Result := 0;
    exit;
  end;

  // Non-zero case
  if abs(C.Re) > abs(C.Im) then begin
    Result := ArcTan(C.Im / C.Re); {-45 to 45 deg, 135 to -135 deg}
    if C.Re < 0 then Result := Result + pi;
  end else begin
    Result := cPid2 - ArcTan(C.Re / C.Im); {45 to 135, -45 to -135}
    if C.Im < 0 then Result := Result + pi;
  end;
  if Result > pi then Result := Result - c2pi;
end;

function ComplexPower(const C: TComplex): TFloat;
// Get the power capacity of the complex number C
begin
    Result := sqr(C.Re) + sqr(C.Im);
end;

end.