FFT算法实现 Delphi / Windows SDK/APIhttp://www.delphi2007.net/DelphiBase/html/delphi_20061210124307209.html
type TComplexData = record //自定义复数类型
vReal,vImaginary: Double;
end;
TComplexDataArray = array of TComplexData; //实数数组
TExtendedArray = array of Extended; //复数数组
function ComplexAdd(c1,c2:TComplexData):TComplexData; //复数相加
function ComplexMinus(c1,c2:TComplexData):TComplexData; //复数相减
function ComplexMultiply(c1,c2:TComplexData):TComplexData; //复数相乘
procedure ComplexArrayWrap(var a1,a2:TComplexDataArray); //复数数组互换
procedure DoFFT(TD:TExtendedArray; var FD:TExtendedArray; r:Integer);
procedure Tmainform.DoFFT(TD:TExtendedArray; var FD:TExtendedArray; r:Integer);
var
i,j,k : Integer; //循环变量
bfsize,p : Integer; //中间变量
angle : Double; //角度
begin
count:= 1 shl r; //傅立叶变换点数
// 分配运算所需存储器
SetLength(W, count div 2);
SetLength(X1,count);
SetLength(X2,count);
// 计算加权系数
for i:=0 to count div 2 -1 do
begin
angle := -i * PI * 2 / count;
W[i].vReal := cos(angle);
W[i].vImaginary := sin(angle);
end;
// 将时域点写入X1
for i:=0 to count-1 do
begin
X1[i].vReal:=TD[i];
X1[i].vImaginary:=0;
end;
// 采用蝶形算法进行快速傅立叶变换
for k:=0 to r-1 do
begin
for j:=0 to (1 shl k)-1 do
begin
bfsize:= 1 shl (r-k);
for i:=0 to bfsize div 2 -1 do
begin
p := j * bfsize;
X2[i+p]:=ComplexAdd(X1[i+p],X1[i+p+bfsize div 2]);
X2[i+p+bfsize div 2]:=ComplexMultiply(
ComplexMinus(X1[i+p],X1[i+p+ bfsize div 2]),
W[i*(1 shl k)]);
end;
end;
ComplexArrayWrap(X1,X2);
end;
// 重新排序
Setlength(Fd,count div 2);
for j:=0 to count div 2 - 1 do
begin
p := 0;
for i:=0 to r-1 do
begin
if (j and (1 shl i))>0 then
p:=p+(1 shl (r-i-1));
end;
FD[j]:=X1[p].vReal;
end;
// 释放内存
Setlength(W,0);
Setlength(X1,0);
Setlength(X2,0);
end;
上面是fft的算法实现代码,但是我在调用的时候,发现TD,FD输入参数有些问题,
我原始数据是一个桥梁检测的信号,一个分量是频率,一个分量是幅值,我要如何通过FFT算法,把这个信号图,分解成为频率图和幅度图?
请教各位高手,谢谢!
网上有很多快速傅立叶变换的例子. 有DELPHI实现的算法库
有源码库的
噢 来凑个热闹
现在的问题是不知道怎么初始化参数,也就是TD和FD,有没有这方面的实例?谢谢!
我的Email是:chenamo5776@126.com