谁有fft代码啊？小弟只有200分，全给了！(200分)

最近数字信号编程遇到个问题，原先编的fft（快速傅立叶变换）只能进行1024*32个数据量的变换。现在想找一个能进行至少1024*512个数据量变换的delphi源代码！请哪位高手给小弟一个！email:chenxiu198228@163.com qq:390617495!希望高手能发小弟一个!

小弟现在课题急需，在线等待啊！ 怎么都没有回答呢？

http://www.simdesign.nl/fft.html,以前有搞过

sonican:你给的这个能支持１０２４*２５６点以上的吗？
你以前用过这个单元吗？

刚刚看了下网页上的说明，好象只支持１０２４点的！

试试下面的代码（将其保存为SkyFFT2.pas）：
{*****************************************************************************
* *
* 下面的FFT（傅立叶变换）函数是从网上找来的！虽然我也学过，但早忘了^o^ *
* *
* Email: iamdream@yeah.net *
* *
*****************************************************************************}

unit SkyFFT2;

interface

uses
Windows, Classes, SysUtils, Math;

type
TComplex = record
Real: Double;
Imag: Double;
end;

PComplexArray = ^TComplexArray;
TComplexArray = array of TComplex;

procedure FFT(Source, Target: TComplexArray
iCount: Integer);
procedure iFFT(Source, Target: TComplexArray
iCount: Integer);

implementation
{
function IsPowerOfTwo ( x: word ): boolean
// 岆 2 腔蹶?
var
i, y: word;
begin
y := 2;
for i := 1 to 31 do begin
if x = y then begin
Result := TRUE;
Exit;
end;
y := y SHL 1;
end;
Result := FALSE;
end;
}

function NumberOfBitsNeeded ( PowerOfTwo: word ): word;
var
i: word;
begin
for i := 0 to 16 do begin
if (PowerOfTwo AND (1 SHL i)) <> 0 then begin
Result := i;
Exit;
end;
end;
Result := 0;
end;


function ReverseBits ( index, NumBits: word ): word;
var
i: word;
begin
Result := 0;
for i := 0 to NumBits-1 do begin
Result := (Result SHL 1) OR (index AND 1);
index := index SHR 1;
end;
end;

procedure DoFFT(Source, Target: TComplexArray
iCount: Integer
AngleNumerator: Double);
var
NumBits, i, j, k, n, BlockSize, BlockEnd: word;
delta_angle, delta_ar: double;
alpha, beta: double;
tr, ti, ar, ai: double;
begin
//if not IsPowerOfTwo(iCount) or (iCount<2) then
// raise Exception.Create('Count is not a positive integer power of 2');

NumBits := NumberOfBitsNeeded (iCount);
for i := 0 to iCount-1 do begin
j := ReverseBits ( i, NumBits );
Target[j] := Source;
end;
BlockEnd := 1;
BlockSize := 2;
while BlockSize <= iCount do begin
delta_angle := AngleNumerator / BlockSize;
alpha := sin ( 0.5 * delta_angle );
alpha := 2.0 * alpha * alpha;
beta := sin ( delta_angle );

i := 0;
while i < iCount do begin
ar := 1.0
(* cos(0) *)
ai := 0.0
(* sin(0) *)

j := i;
for n := 0 to BlockEnd-1 do begin
k := j + BlockEnd;
tr := ar*Target[k].Real - ai*Target[k].Imag;
ti := ar*Target[k].Imag + ai*Target[k].Real;
Target[k].Real := Target[j].Real - tr;
Target[k].Imag := Target[j].Imag - ti;
Target[j].Real := Target[j].Real + tr;
Target[j].Imag := Target[j].Imag + ti;
delta_ar := alpha*ar + beta*ai;
ai := ai - (alpha*ai - beta*ar);
ar := ar - delta_ar;
INC(j);
end;
i := i + BlockSize;
end;
BlockEnd := BlockSize;
BlockSize := BlockSize SHL 1;
end;
end;

procedure FFT(Source, Target: TComplexArray
iCount: Integer);
begin
DoFFT(Source, Target, iCount, 2 * PI);
end;

procedure iFFT(Source, Target: TComplexArray
iCount: Integer);
var
i: Integer;
begin
DoFFT(Source, Target, iCount, -2 * PI);
(* Normalize the resulting time samples... *)
for i := 0 to iCount -1 do begin
Target.Real := Target.Real / iCount;
Target.Imag := Target.Imag / iCount;
end;
end;

end.

dreamisx: 你的代码我调试了下,它只能支持1024*16 点运算.
谢谢哪位可以给我能支持1024*256以上点运算的代码啊,小弟现在毕业环节急死了!跪求各位了!

不知道这个是不是一样只支持1024*16的，基-2算法，我看不懂...还是发你看看吧

(* 快速付里叶变换
pr-双精度浮点数一维数组，长度是n。
当l=0时，存放n个采样输入的实部，返回时存放傅理叶变换的模；
当l=1时，存放傅理叶变换的n个实部，返回存放逆付里叶变换的模
pi-双精度浮点数一维数组，长度是n。
当l=0时，存放n个采样输入的虚部，返回时存放傅理叶变换的幅角；
当l=1时，存放傅理叶变换的n个虚部，返回存放逆付里叶变换的幅角。幅角单位为实数，不是角度
n-整数变量。输入的点数。
k-整数变量。满足n=2^k。
fr-双精度浮点数一维数组，长度是n。
当l=0时，返回时存放傅理叶变换的实部；
当l=1时，返回逆付里叶变换的实部
fi-双精度浮点数一维数组，长度是n。
当l=0时，返回时存放傅理叶变换的虚部；
当l=1时，返回逆付里叶变换的虚部
l-整数变量。
当l=0时，表示本函数计算傅里叶变换；
当l=1时，表示本函数计算逆付里叶变换
il-整数变量。
当il=0时，表示本函数不要求计算傅里叶或逆付里叶变换的模和幅角；
当il=1时，表示本函数要求计算傅里叶变换或逆付里叶变换的模和幅角 *)
procedure kkfft(var pr: array of double
var pi: array of double
n: integer
k: integer;
var fr: array of double
var fi: array of double
l: integer;il: integer);
var
it,m,is0,i,j,nv,l0: integer;
p,q,s,vr,vi,poddr,poddi: double;
begin
for it := 0 to n - 1 do
begin
m := it;
is0 := 0;
for i := 0 to k-1 do
begin
j := m div 2;
is0 := 2 * is0 + (m - 2 * j);
m := j;
end;
fr[it] := pr[is0];
fi[it] := pi[is0];
end;
pr[0] := 1.0;
pi[0] := 0.0;
p := 6.283185306 / (1.0 * n);
pr[1] := cos(p);
pi[1] := -sin(p);
if (l <> 0) then pi[1] := -pi[1];
for i := 2 to n-1 do
begin
p := pr[i-1] * pr[1];
q := pi[i-1] * pi[1];
s := (pr[i-1] + pi[i-1]) * (pr[1] + pi[1]);
pr := p - q;
pi := s - p - q;
end;
it := 0;
while it <= n-2 do
begin
vr := fr[it];
vi := fi[it];
fr[it] := vr + fr[it+1];
fi[it] := vi + fi[it+1];
fr[it+1] := vr - fr[it+1];
fi[it+1] := vi - fi[it+1];
it := it + 2;
end;
m := n div 2;
nv := 2;
for l0 := k-2 downto 0 do
begin
m := m div 2;
nv := 2 * nv;
it := 0;
while it <= (m-1)*nv do
begin
for j := 0 to (nv div 2)-1 do
begin
p := pr[m*j] * fr[it+j+nv div 2];
q := pi[m*j] * fi[it+j+nv div 2];
s := pr[m*j] + pi[m*j];
s := s * (fr[it+j+nv div 2] + fi[it+j+nv div 2]);
poddr := p - q;
poddi := s - p - q;
fr[it+j+nv div 2] := fr[it+j] - poddr;
fi[it+j+nv div 2] := fi[it+j] - poddi;
fr[it+j] := fr[it+j] + poddr;
fi[it+j] := fi[it+j] + poddi;
end;
it := it + nv;
end;
end;
{if l <> 0 then
begin
for i := 0 to n-1 do
begin
fr := fr / (1.0 * n);
fi := fi / (1.0 * n);
end;
end;
if il <> 0 then
begin
for i := 0 to n-1 do
begin
pr := sqrt(fr*fr+fi*fi);
if abs(fr) < 0.000001*abs(fi) then
begin
if fi*fr > 0 then
//pi=90.0;
pi := 6.283185306 / 4
else
//pi=-90.0;
pi := -6.283185306 / 4;
end
else
//pi = atan(fi/fr) * 360.0 / 6.283185306;
pi := ArcTan(fi/fr);
end
end
}
end;

木桩:你的代码我调试了下，可以支持1024*256点运算，但是这个代码的精度不高，点数越多，后面的误差就越大。我试着改掉里面的一些代码，但是精度的问题还是没解决！
因为我的算法要求精度是至少小数后面6位！

继续在线等待更好的算法！

谢谢大家对本问题的回答！感谢大家！分已送出！

结题！