Recent

Author Topic: [SOLVED] Is there FFT and hanning window code for Lazarus?  (Read 12037 times)

mas steindorff

  • Sr. Member
  • ****
  • Posts: 468
[SOLVED] Is there FFT and hanning window code for Lazarus?
« on: September 01, 2011, 08:42:49 pm »
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
« Last Edit: September 12, 2011, 07:07:29 pm by mas steindorff »
windows 7/10 - laz 2.0 / 1.2.6 general releases

Laksen

  • Hero Member
  • *****
  • Posts: 686
    • J-Software
Re: Is there FFT and hanning window code for Lazarus?
« Reply #1 on: September 02, 2011, 12:09:04 am »
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

  • Sr. Member
  • ****
  • Posts: 468
Re: Is there FFT and hanning window code for Lazarus?
« Reply #2 on: September 02, 2011, 01:15:35 am »
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!
windows 7/10 - laz 2.0 / 1.2.6 general releases

Ocye

  • Hero Member
  • *****
  • Posts: 518
    • Scrabble3D
Re: Is there FFT and hanning window code for Lazarus?
« Reply #3 on: September 12, 2011, 02:45:59 pm »
"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: [Select]

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;
Lazarus 1.7 (SVN) FPC 3.0.0

Blaazen

  • Hero Member
  • *****
  • Posts: 3034
  • POKE 54296,15
    • Eye-Candy Controls
Re: Is there FFT and hanning window code for Lazarus?
« Reply #4 on: September 12, 2011, 04:44:27 pm »
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.
Lazarus 2.3.0 rmain-2_3-280-g5db1922e37 FPC 3.3.1 r40507 x86_64-linux-qt Chakra, Qt 4.8.7/5.13.2, Plasma 5.17.3
Lazarus 1.8.2 r57369 FPC 3.0.4 i386-win32-win32/win64 Wine 3.21

Try Eye-Candy Controls: https://sourceforge.net/projects/eccontrols/files/

 

TinyPortal © 2005-2018