研究傅立叶变换，找到的代码贴出来共享。。

[ksaiy]

unit FFT1;
{ Unit for Real Fast Fourier Transformation
Julian Ziersch, CIS = 100744,2101
------------------------------------------------------------------------
real_fft(tau,
tau_2: Integer;
VAR y:array of single;
direkt:Boolean);

******* Usage to perform FFT : real_fft()
direkt=FALSE
size = tau = 2^tau_2
y : array[1..size] of single = Input

Output: y[0]=a0
y[k]=a(k+1)/2
with k=1,3,...,tau_2-1

if you need the real values you have to use the procedure
fft_makereal(var y:array of single;size:Integer);
Input: y : array[1..size] of single = Output of previous FFT
Output: y[0] = a0
y[1..size/2] = abs(FFT)

******* Usage to perform reverse FFT : real_fft()
direkt=TRUE
size = tau = 2^tau_2
y : array[1..size] of single = Output of previous FFT

Output: y[1..size] = Real Values

You can delete some of the y[...] values before you preform the reverse FFT.
If you do so and use WAV-Values for Input, you will get a digital filter.
You can use the procedure
delete_fft(from,to,count:Integer; VAR y:array of single);
to delete y[1..from-1] and y[to+1..size-1]
---------------------------------------------------------------------------}

{$U-} (* Pentium Sicherheits-Abfrage abschalten *)

{Jede Megen Optimierung moeglich, precis Abfrage, fft,
fft braucht fuer 512 Werte derzeit 28ms (Pentium 75) }

interface
uses SysUtils;

(* Real-FFT : *)
procedure real_fft(tau,tau_2:Integer;VAR y:array of single;direkt:Boolean);
procedure fft_makereal(var y:array of single;Size:Integer);
procedure delete_fft(von,bis,anz:Integer;VAR y:array of single);
procedure FFT_World(i:integer; A:array of single;
var XMax:Integer;
var YMax,YMin:single;
precis:single);
function DoubleToLong(d:single):Longint;
implementation

{Reele FFT mit Umkehrfunktion ----------------------------------------
Eingabeparameter:
direct:Boolean FALSE: Berechnung der diskkreten Fourierkoeffizienten
TRUE: Berechnung der Funktionswerte
tau:Integer 2 hoch tau ist Anzahl der Funktionswerte
tau_2:Integer Anzahl der Funktionswerte
VAR y: enthaelt Werte und zwar fuer direct=
array[tau_2] FALSE : Funktionswerte
og single TRUE: Diskrete Fourierkoeffizienten
y[0]=a0
y[k]=a(k+1)/2
mit k=1,3,...,tau_2-1
also in der Reihenfolge a0 a1 b1 a2 b2
Ausgabeparameter:
y:array[tau_2] fuer direct=
FALSE: diskrete Fourierkoeffizienten (s.o)
TRUE: Funktionswerte
Benoetigt: sin cos und PI:=3.141...
tau_2 MUSS 2 hoch n sein !!!
-----------------------------------------------------------------------}

procedure real_fft(tau,tau_2:Integer;VAR y:array of single;direkt:Boolean);
var
md_2,md_4, sigma,sigma_2, j_2, min_n,n_min_0,n_min_1 : Integer;
ind_1, ind_2, k, j, n, l, k_2, tau_2_2k : Integer;
faktor, arg, arg_m, arg_md_2, y_hilf : single;
ew_r, ew_i, eps_r, eps_i, ur, ui, wr, wi : single;
rett, hilf_1, hilf_2, hilf_3, hilf_4 : single;
begin
md_2 := tau_2 div 2;
md_4 := md_2 div 2;
faktor := 1.0 / md_2;
arg_md_2 := 2.0 * PI * faktor;
arg_m := 0.5 * arg_md_2;
if direkt=TRUE then faktor := 1.0;

(* Zusammenfassen der reelen Daten zur FFT halber Laenge : *
* ........................................................*)
if direkt=TRUE then
begin
y_hilf := y[1];
y[1] := y[0] - y[tau_2-1];
y[0] := y[0] + y[tau_2-1];

ew_r := cos(arg_m); (* ew_r/ew_i = Real/Imag.teil der *)
ew_i := sin(arg_m); (* tau-ten Einheitswurzel *)

eps_r := 1.0; (* Real bzw. Imaginaerteil der *)
eps_i := 0.0; (* tau_2-ten Einheitswurzel *)

k := 1;
while k<md_4 do (* for(k=1;k<md_4;k++) *)
begin
k_2 := 2 * k;
tau_2_2k := tau_2 - k_2;
rett := eps_r;
eps_r := rett * ew_r - eps_i * ew_i;
eps_i := rett * ew_i + eps_i * ew_r;

hilf_1 := 0.5 * (eps_r * (y_hilf - y[tau_2_2k-1])
+ eps_i * (y[k_2] + y[tau_2_2k]));
hilf_2 := 0.5 * (eps_i * (y_hilf - y[tau_2_2k-1])
- eps_r * (y[k_2] + y[tau_2_2k]));

hilf_3 := 0.5 * (y_hilf + y[tau_2_2k-1]);
hilf_4 := 0.5 * (y[k_2] - y[tau_2_2k]);
y_hilf := y[k_2+1];
y[k_2] := hilf_3 - hilf_2;
y[k_2+1] := hilf_1 - hilf_4;
y[tau_2_2k] := hilf_2 + hilf_3;
y[tau_2_2k+1] := hilf_1 + hilf_4;
inc(k);
end; (*while*)
y[md_2+1] := y[md_2];
y[md_2] := y_hilf;
end; (* if direkt *)


(* Umspeicherung mit der Bit-Umkehrfunktion, gleichzeitig *
* Normierung fals direkt=FALSE ..........................*)
k := 0; j := 0;
while j<md_2 do (* for(j=k=0;j<md_2;j++) *)
begin
k := j;
sigma := 0;
n :=1;
while n<tau do (* for(n=1,sigma=0.0;n<tau;n++) *)
begin
sigma := 2 * sigma + (k and 1);
k := k div 2;
n := n+1;
end; (*while n *)

if j<=sigma then
begin
j_2 := 2 * j;
sigma_2 := 2 * sigma;
ur := y[j_2];
ui := y[j_2+1];
y[j_2] := y[sigma_2] * faktor;
y[j_2+1]:= y[sigma_2+1] * faktor;
y[sigma_2] := ur * faktor;
y[sigma_2+1] := ui * faktor;
end; (*if*)
inc(j);
end; (*while j*)

(* Durchfuehrung der FFT halber Laenge .....................*
* min_n : 2^(tau-1-n)
* n_min_1 : 2^(n-1)
* n_min_0 : 2^n
* wr / wi : Real/Imaginaerteil ((tau_2)/2 ^ 2^min_n
* eps_r/i : Real/imaginaerteil ((tau_2)/2)^(l*2^min_n)
*.........................................................*)

[fansnaf]

学习

[ksaiy]

min_n := md_2;
n_min_1 := 1;
n := 1;
while n<tau do (* for(n=1;n<tau;n++) *)
begin
min_n := min_n div 2;
n_min_0 := 2 * n_min_1;
arg := arg_md_2 * min_n;
wr := cos(arg);
if direkt=TRUE then wi := sin(arg)
else wi := -sin(arg);
eps_r := 1.0;
eps_i := 0.0;

l := 0;
while l<n_min_1 do (* for(l=0;l<n_min_1;l++) *)
begin
j := 0;
while j<=md_2-n_min_0 do (* for(j=0;j<=md2-n_min0;j+=n_min_0) *)
begin
ind_1 := (j + l) * 2;
ind_2 := ind_1 + n_min_0;
ur := y[ind_2] * eps_r - y[ind_2+1] * eps_i;
ui := y[ind_2] * eps_i + y[ind_2+1] * eps_r;
y[ind_2] := y[ind_1] - ur;
y[ind_2+1] := y[ind_1+1] - ui;
y[ind_1] := y[ind_1] + ur;
y[ind_1+1] := y[ind_1+1] + ui;
j := j+n_min_0;
end; (* while j*)
rett := eps_r;
eps_r := rett * wr - eps_i * wi;
eps_i := rett * wi + eps_i * wr;
inc(l);
end; (*while l*)
n_min_1 := n_min_0;
inc(n);
end; (*while n*)

(* Trennung zusammengefasst transformierter Daten falls *
* direkt=FALSE.........................................*)

if direkt=FALSE then
begin
y_hilf := y[tau_2-1];
y[tau_2-1] := 0.5 * (y[0] - y[1]);
y[0] := 0.5 * (y[0] + y[1]);

(* ew_r/i = Real/Imag von tau_2-ten Einheitswurzel *)
ew_r := cos(arg_m);
ew_i := -sin(arg_m);
(* eps_r/i = R/I von (tau_2-ten Einheitswurzel)^k *)
eps_r := 1.0;
eps_i := 0.0;

k := 1;
while k<md_4 do (* for(k=1;k<md_4;k++) *)
begin
rett := eps_r;
eps_r := rett * ew_r - eps_i * ew_i;
eps_i := rett * ew_i + eps_i * ew_r;
ind_1 := k * 2;
ind_2 := tau_2 - ind_1;
hilf_1 := 0.5 * (eps_i * (y[ind_1] - y[ind_2])
+ eps_r *(y[ind_1+1] + y_hilf));
hilf_2 := 0.5 * (eps_r * (y[ind_1] - y[ind_2])
- eps_i *(y[ind_1+1] + y_hilf));
hilf_3 := 0.5 * (y[ind_1] + y[ind_2]);
hilf_4 := 0.5 * (y[ind_1+1]-y_hilf);
y_hilf := y[ind_2-1];
y[ind_1-1] := hilf_1 + hilf_3;
y[ind_1] := hilf_2 - hilf_4;
y[ind_2-1] := hilf_3 - hilf_1;
y[ind_2] := hilf_2 + hilf_4;
inc(k);
end; (* while *)
y[md_2-1] := y[md_2];
y[md_2] := y_hilf;
end; (*if direkt*)
end; (* real_fft *)

(* Fuert Bandpass aus *)
procedure delete_fft(von,bis,anz:Integer;VAR y:array of single);
var
i:Integer;
begin
for i:=1 to von-1 do y[i]:=0.0;
for i:=bis+1 to anz-1 do y[i]:=0.0;
end; (* delete *)

(* Ermittelt X und Y Bereiche fuer Bildschirmausgabe *)
procedure FFT_World(i:integer; A:array of single;
var XMax:Integer;
var YMax,YMin:single;
precis:single);
var
j:integer;
begin
XMax:=1;
YMax:=A[2];
YMin:=YMax;
for j:=1 to i do
begin
if abs(A[j])>precis then XMax:=j;
if A[j]>YMax then YMax:=A[j]
else if A[j]<YMin then YMin:=A[j];
end;
end;



(* Wandelt reele Zahl in Integer um *)
function DoubleToLong(d:single):Longint;
begin
DoubleToLong:=Round(d);
end; (*DoubleToLong*)

{ Wandelt FFT Ergebnis r[i],i[i] Array in Realbetraege
Groesse ist dann Size/2!
y[0] bleibt
y[1,3,5 ...] :=betrag }
procedure fft_makereal(var y:array of single;Size:Integer);
var
i,j:Integer;
begin
i:=1;
j:=0;
while i<Size do
begin
y[j]:=sqrt(sqr(y[i])+sqr(y[i+1]));
inc(i);
inc(i);
inc(j);
end;
end; (* fft_makereal *)
end.

[huanyi]

不错，学习，呵呵

[beyondtkl]

XUEXI

[XuDunYu]

UP

[psp2003]

关注