[SOLVED] Is there FFT and hanning window code for Lazarus?

Forum > Packages and Libraries

(1/1)

mas steindorff:
I'm looking for an FFT and "Hanning windowing" signal processing code to add to one of my projects.

I remember (and found) my old Numerical Recipes for pascal tool box I got back in the D4 days but the code needs work to run with the today's FPC. [and I would need to find an old floppy disk drive to access the code :( ]

It seems to me this code has been around long enough to be built in to some math library some where but I can't find it.

I searched and found out there uses to be an FFT example included with Lazarus back with fpc 2.2.2 but it was a dll interface. That example is not included with 2.4.2 and the website it calls out longer supports free pascal .

I though I'd ask here before I reinvent the wheel.

MAS

Laksen:
I made some functions for it once making use of high level functionality. Surprisingly it was still really fast

There's a lot of extra functionality, but if you want you can just pry out the stuff you need
http://j-software.dk/dsp.zip

mas steindorff:
Thanks Laksen, it looks like the code also has some other vector functions I was figuring I would have to write myself like mean() and median() as well ::).

depending on how fast my end code is, this may be the 2nd to final nail for some python code >:D I' hoping to do away with. Now all I need is a HDF5 file writer as the final nail!

Ocye:
"Numerical recipes" published the whole code not only in C/C++ but as well in standard Pascal. Have a look for the keywords with your preferred search engine.
BTW: I know Hann and Hamming window functions (and some other) but not Hanning. Some old stuff that needs to be checked carefully, perhaps it helps a little bit:

--- Code: ---
procedure Windowing(var data:array of double;Taps:integer;WF:TWindowFunktion);
var i:integer;
begin
for i:=0 to high(data) do
case WF of
wf_None : ;
wf_Hamming : data[i]:=data[i]*(0.54-0.46*cos((i*2*PI)/Taps));
wf_Blackman : data[i]:=data[i]*(0.42-0.5*cos((i*2*PI)/Taps)+0.08*cos((i*4*PI)/Taps));
wf_KaiserBessel : data[i]:=data[i]*(0.4021-0.4986*cos((i*2*PI)/Taps)+0.0981*cos((i*4*PI)/Taps)-0.0012*cos((i*6*PI)/Taps));
wf_FlatTop : data[i]:=data[i]*(0.2155-0.4159*cos((i*2*PI)/Taps)+0.2780*cos((i*4*PI)/Taps)-0.0836*cos((i*6*PI)/Taps)+0.0070*cos((i*8*PI)/Taps));
end;//case
end;

procedure DFT(data: array of RealType;isign:integer);
var i,i1,i2,i3,i4,np3, Count : integer;
c1,c2,h1r,h2r,h1i,h2i : RealType;
wr,wi,wpr,wpi,wtemp,theta : RealType;
begin
Count:=length(data);
if (Count mod 4<>0) then raise EMathError.Create('Fehler in DFT (N kein Vielfaches von 4)');
c1:=0.5;
theta:=PI/(Count/2);
if (isign=1) then
begin
c2:=-0.5;
four(data,1);
end else
begin
c2:=0.5;
theta:=-theta;
end;
wtemp:=sin(0.5*theta);
wpr:=-2.0*wtemp*wtemp;
wpi:=sin(theta);
wr:=1.0+wpr;
wi:=wpi;
np3:=Count+3;
for i:=2 to Count div 4 do
begin
// i4:=1+(i3=np3-(i2=1+(i1=i+i-1)));
h1r:=c1*(data[i1]+data[i3]);
h1i:=c1*(data[i2]-data[i4]);
h2r:=-c2*(data[i2]+data[i4]);
h2i:=c2*(data[i1]-data[i3]);
data[i1]:=h1r+wr*h2r-wi*h2i;
data[i2]:= h1i+wr*h2i+wi*h2r;
data[i3]:= h1r-wr*h2r+wi*h2i;
data[i4]:=-h1i+wr*h2i+wi*h2r;
// wr:=(wtemp=wr)*wpr-wi*wpi+wr;
wi:=wi*wpr+wtemp*wpi+wi;
end;
if (isign=1) then
begin
// data[1]:=(h1r=data[1])+data[2];
data[2]:=h1r-data[2];
end else
begin
// data[1]:=c1*((h1r=data[1])+data[2]);
data[2]:=c1*(h1r-data[2]);
four(data,-1);
end;
end;

procedure SWAP(var a, b:double);
var temp:double;
begin
temp:=(a);
a:=b;
b:=temp;
end;

procedure FFT(var data: array of double; nn:integer; invers: boolean);
(* Programs using routine FOUR1 must define type
TYPE
gldarray = ARRAY [1..nn2] OF real;
in the calling routine, where nn2=nn+nn. *)
var ii,jj,n,mmax,m,j,istep,i: integer;
wtemp,wr,wpr,wpi,wi,theta: double;
tempr,tempi: real;
begin
n := 2*nn;
j := 1;
for ii := 1 to nn do
begin
i := 2*ii-1;
if (j > i) then
begin
tempr := data[j];
tempi := data[j+1];
data[j] := data[i];
data[j+1] := data[i+1];
data[i] := tempr;
data[i+1] := tempi
end;
m := n div 2;
while ((m >= 2) and (j > m)) do
begin
j := j-m;
m := m div 2
end;
j := j+m
end;
mmax := 2;
while (n > mmax) do
begin
istep := 2*mmax;
if invers then
theta := 6.28318530717959/mmax else
theta := 6.28318530717959/-mmax;
wpr := -2.0*sqr(sin(0.5*theta));
wpi := sin(theta);
wr := 1.0;
wi := 0.0;
for ii := 1 to (mmax div 2) do
begin
m := 2*ii-1;
for jj := 0 to ((n-m) div istep) do
begin
i := m + jj*istep;
j := i+mmax;
tempr := real(wr)*data[j]-real(wi)*data[j+1];
tempi := real(wr)*data[j+1]+real(wi)*data[j];
data[j] := data[i]-tempr;
data[j+1] := data[i+1]-tempi;
data[i] := data[i]+tempr;
data[i+1] := data[i+1]+tempi
end;
wtemp := wr;
wr := wr*wpr-wi*wpi+wr;
wi := wi*wpr+wtemp*wpi+wi
end;
mmax := istep
end;
end;

procedure four(data : array of RealType; isign: integer);
var tmpdata : array of RealType;
alpha, beta, theta,
temp,sw,cw : double;
i,j,Count : integer;
begin
Count:=length(data);
setlength(tmpdata,2*Count);
ZeroMemory(tmpdata,sizeof(tmpdata));
for i:=0 to Count-1 do
begin
theta:=isign*2.0*PI*(i/Count);
temp:=sin(0.5*theta);
alpha:=2.0*sqr(temp);
beta:=sin(theta);
cw:=1.0;
sw:=0.0;
for j:=0 to Count-1 do
begin
tmpdata[2*i]:=tmpdata[2*i]+(cw*data[2*j]-sw*data[2*j+1]);
tmpdata[2*i+1]:=tmpdata[2*i+1]+(cw*data[2*j+1]+sw*data[2*j]);
// cw:=(temp=cw)-(alpha*cw+beta*sw);
sw:=sw-(alpha*sw-beta*temp);
end;
end;
if (isign=-1) then for i:=0 to 2*Count-1 do tmpdata[i]:=tmpdata[i]/Count;
for i:=0 to 2*Count-1 do data[i]:=tmpdata[i];
setlength(tmpdata,0);
end;

--- End code ---

Blaazen:
Of some reason, "Hanning" is name of function in Matlab. His name was Julius von Hann.

Thanks for code Ocye, it is useful.

I used windowing for design of FIR filters (linear frequency filters). If I will have time I will make the unit usable (user freindly) and I will publish it too.

Navigation

[0] Message Index