注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

鑫淼梦园的博客

圆你的梦想 从这里开始

 
 
 

日志

 
 

Delphi傅立叶变换  

2014-10-18 01:03:49|  分类: delphixe7 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
傅立叶变换
unit cplx;
interface
type
PReal = ^TReal;
TReal = extended;

PComplex = ^TComplex;
TComplex = record
r : TReal;
i : TReal;
end;

function MakeComplex(x, y: TReal): TComplex;
function Sum(x, y: TComplex) : TComplex;
function Difference(x, y: TComplex) : TComplex;
function Product(x, y: TComplex): TComplex;
function TimesReal(x: TComplex; y: TReal): TComplex;
function PlusReal(x: TComplex; y: TReal): TComplex;
function EiT(t: TReal):TComplex;
function ComplexToStr(x: TComplex): string;
function AbsSquared(x: TComplex): TReal;

implementation

uses SysUtils;

function MakeComplex(x, y: TReal): TComplex;
begin
with result do
begin
r:=x;
i:= y;
end;
end;

function Sum(x, y: TComplex) : TComplex;
begin
with result do
begin
r:= x.r + y.r;
i:= x.i + y.i;
end;
end;

function Difference(x, y: TComplex) : TComplex;
begin
with result do
begin
r:= x.r - y.r;
i:= x.i - y.i;
end;
end;

function EiT(t: TReal): TComplex;
begin
with result do
begin
r:= cos(t);
i:= sin(t);
end;
end;

function Product(x, y: TComplex): TComplex;
begin
with result do
begin
r:= x.r * y.r - x.i * y.i;
i:= x.r * y.i + x.i * y.r;
end;
end;

function TimesReal(x: TComplex; y: TReal): TComplex;
begin
with result do
begin
r:= x.r * y;
i:= x.i * y;
end;
end;

function PlusReal(x: TComplex; y: TReal): TComplex;
begin
with result do
begin
r:= x.r + y;
i:= x.i;
end;
end;

function ComplexToStr(x: TComplex): string;
begin
result:= FloatToStr(x.r)
+ ' + '
+ FloatToStr(x.i)
+ 'i';
end;

function AbsSquared(x: TComplex): TReal;
begin
result:= x.r*x.r + x.i*x.i;
end;

end.

--------------------------------------------------------------------------------
--------------------------------------------------------------------------------

unit cplxfft1;

interface

uses Cplx;

type
PScalar = ^TScalar;
TScalar = TComplex; {Making conversion to real version easier}

PScalars = ^TScalars;
TScalars = array[0..High(integer) div SizeOf(TScalar) - 1]
of TScalar;

const

TrigTableDepth: word = 0;
TrigTable : PScalars = nil;

procedure InitTrigTable(Depth: word);

procedure FFT(Depth: word;
Src: PScalars;
Dest: PScalars);

{REQUIRES allocating

(integer(1) shl Depth) * SizeOf(TScalar)

bytes for Src and Dest before call!}

implementation

procedure DoFFT(Depth: word;
Src: PScalars;
SrcSpacing: word;
Dest: PScalars);
{the recursive part called by FFT when ready}
var j, N: integer; Temp: TScalar; Shift: word;
begin
if Depth = 0 then
begin
Dest^[0]:= Src^[0];
exit;
end;

N:= integer(1) shl (Depth - 1);

DoFFT(Depth - 1, Src, SrcSpacing * 2, Dest);
DoFFT(Depth - 1, @Src^[SrcSpacing], SrcSpacing * 2, @Dest^[N] );

Shift:= TrigTableDepth - Depth;

for j:= 0 to N - 1 do
begin
Temp:= Product(TrigTable^[j shl Shift],
Dest^[j + N]);
Dest^[j + N]:= Difference(Dest^[j], Temp);
Dest^[j] := Sum(Dest^[j], Temp);
end;

end;

procedure FFT(Depth: word;
Src: PScalars;
Dest: PScalars);
var j, N: integer; Normalizer: extended;
begin

N:= integer(1) shl depth;

if Depth TrigTableDepth then
InitTrigTable(Depth);

DoFFT(Depth, Src, 1, Dest);

Normalizer:= 1 / sqrt(N) ;

for j:=0 to N - 1 do
Dest^[j]:= TimesReal(Dest^[j], Normalizer);

end;

procedure InitTrigTable(Depth: word);
var j, N: integer;
begin

N:= integer(1) shl depth;
ReAllocMem(TrigTable, N * SizeOf(TScalar));
for j:=0 to N - 1 do
TrigTable^[j]:= EiT(-(2*Pi)*j/N);
TrigTableDepth:= Depth;

end;

initialization

;

finalization
ReAllocMem(TrigTable, 0);

end.

--------------------------------------------------------------------------------
--------------------------------------------------------------------------------

unit DemoForm;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls;

type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
Edit1: TEdit;
Label1: TLabel;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form1: TForm1;

implementation

{$R *.DFM}

uses cplx, cplxfft1, MMSystem;

procedure TForm1.Button1Click(Sender: TObject);
var j: integer; s:string;
src, dest: PScalars;
norm: extended;
d,N,count:integer;
st,et: longint;
begin

d:= StrToIntDef(edit1.text, -1) ;
if d <1 then
raise exception.Create('depth must be a positive integer');


N:= integer(1) shl d ;

GetMem(Src, N*Sizeof(TScalar));
GetMem(Dest, N*SizeOf(TScalar));

for j:=0 to N-1 do
begin
src^[j]:= MakeComplex(random, random);
end;

begin

st:= timeGetTime;
FFT(d, Src, dest);
et:= timeGetTime;

end;

Memo1.Lines.Add('N = ' + IntToStr(N));
Memo1.Lines.Add('expected norm: ' +#9+ FloatToStr(N*2/3));

norm:=0;
for j:=0 to N-1 do norm:= norm + AbsSquared(src^[j]);
Memo1.Lines.Add('Data norm: '+#9+FloatToStr(norm));
norm:=0;
for j:=0 to N-1 do norm:= norm + AbsSquared(dest^[j]);
Memo1.Lines.Add('FT norm: '+#9#9+FloatToStr(norm));


Memo1.Lines.Add('Time in FFT routine: '+#9
+ inttostr(et - st)
+ ' ms.');
Memo1.Lines.Add(' ');

FreeMem(Src);
FreeMem(DEst);
end;

end.

--------------------------------------------------------------------------------
**** The real version:

****
--------------------------------------------------------------------------------

unit cplxfft2;

interface

type
PScalar = ^TScalar;
TScalar = extended;

PScalars = ^TScalars;
TScalars = array[0..High(integer) div SizeOf(TScalar) - 1]
of TScalar;

const

TrigTableDepth: word = 0;
CosTable : PScalars = nil;
SinTable : PScalars = nil;

procedure InitTrigTables(Depth: word);

procedure FFT(Depth: word;
SrcR, SrcI: PScalars;
DestR, DestI: PScalars);

{REQUIRES allocating

(integer(1) shl Depth) * SizeOf(TScalar)

bytes for SrcR, SrcI, DestR and DestI before call!}

implementation

procedure DoFFT(Depth: word;
SrcR, SrcI: PScalars;
SrcSpacing: word;
DestR, DestI: PScalars);
{the recursive part called by FFT when ready}
var j, N: integer; 
TempR, TempI: TScalar;
Shift: word;
c, s: extended;
begin
if Depth = 0 then
begin
DestR^[0]:= SrcR^[0];
DestI^[0]:= SrcI^[0];
exit;
end;

N:= integer(1) shl (Depth - 1);

DoFFT(Depth - 1, SrcR, SrcI, SrcSpacing * 2, DestR, DestI);
DoFFT(Depth - 1,
@SrcR^[srcSpacing],
@SrcI^[SrcSpacing],
SrcSpacing * 2,
@DestR^[N],
@DestI^[N]);

Shift:= TrigTableDepth - Depth;

for j:= 0 to N - 1 do
begin
c:= CosTable^[j shl Shift];
s:= SinTable^[j shl Shift];

TempR:= c * DestR^[j + N] - s * DestI^[j + N];
TempI:= c * DestI^[j + N] + s * DestR^[j + N];

DestR^[j + N]:= DestR^[j] - TempR;
DestI^[j + N]:= DestI^[j] - TempI;

DestR^[j]:= DestR^[j] + TempR;
DestI^[j]:= DestI^[j] + TempI;
end;

end;

procedure FFT(Depth: word;
SrcR, SrcI: PScalars;
DestR, DestI: PScalars);
var j, N: integer; Normalizer: extended;
begin

N:= integer(1) shl depth;

if Depth TrigTableDepth then
InitTrigTables(Depth);

DoFFT(Depth, SrcR, SrcI, 1, DestR, DestI);

Normalizer:= 1 / sqrt(N) ;

for j:=0 to N - 1 do
begin
DestR^[j]:= DestR^[j] * Normalizer;
DestI^[j]:= DestI^[j] * Normalizer;
end;

end;

procedure InitTrigTables(Depth: word);
var j, N: integer;
begin

N:= integer(1) shl depth;
ReAllocMem(CosTable, N * SizeOf(TScalar));
ReAllocMem(SinTable, N * SizeOf(TScalar));
for j:=0 to N - 1 do
begin
CosTable^[j]:= cos(-(2*Pi)*j/N);
SinTable^[j]:= sin(-(2*Pi)*j/N);
end;
TrigTableDepth:= Depth;

end;

initialization

;

finalization
ReAllocMem(CosTable, 0);
ReAllocMem(SinTable, 0);

end.

--------------------------------------------------------------------------------
--------------------------------------------------------------------------------

unit demofrm;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics,
Controls, Forms, Dialogs, cplxfft2, StdCtrls;

type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
Edit1: TEdit;
Label1: TLabel;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form1: TForm1;

implementation

{$R *.DFM}

uses MMSystem;

procedure TForm1.Button1Click(Sender: TObject);
var SR, SI, DR, DI: PScalars;
j,d,N:integer;
st, et: longint;
norm: extended;
begin

d:= StrToIntDef(edit1.text, -1) ;
if d <1 then
raise exception.Create('depth must be a positive integer');

N:= integer(1) shl d;

GetMem(SR, N * SizeOf(TScalar));
GetMem(SI, N * SizeOf(TScalar));
GetMem(DR, N * SizeOf(TScalar));
GetMem(DI, N * SizeOf(TScalar));

for j:=0 to N - 1 do
begin
SR^[j]:=random;
SI^[j]:=random;
end;

st:= timeGetTime;
FFT(d, SR, SI, DR, DI);
et:= timeGetTime;

memo1.Lines.Add('N = '+inttostr(N));
memo1.Lines.Add('expected norm: '+#9+FloatToStr(N*2/3));

norm:=0;
for j:=0 to N - 1 do
norm:= norm + SR^[j]*SR^[j] + SI^[j]*SI^[j];
memo1.Lines.Add('Data norm: '+#9+FloatToStr(norm));

norm:=0;
for j:=0 to N - 1 do
norm:= norm + DR^[j]*DR^[j] + DI^[j]*DI^[j];
memo1.Lines.Add('FT norm: '+#9#9+FloatToStr(norm));

memo1.Lines.Add('Time in FFT routine: '+#9+inttostr(et-st));
memo1.Lines.add('');
(*for j:=0 to N - 1 do
Memo1.Lines.Add(FloatToStr(SR^[j])
+ ' + '
+ FloatToStr(SI^[j])
+ 'i');

for j:=0 to N - 1 do
Memo1.Lines.Add(FloatToStr(DR^[j])
+ ' + '
+ FloatToStr(DI^[j])
+ 'i');*)

FreeMem(SR, N * SizeOf(TScalar));
FreeMem(SI, N * SizeOf(TScalar));
FreeMem(DR, N * SizeOf(TScalar));
FreeMem(DI, N * SizeOf(TScalar));
end;

end.
  评论这张
 
阅读(372)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017