Fast Fourier Transform

  1. Hi all ,
    This FFT program may help some of you in your Engineering studies , it was a lot of fun to write and use ..
    I wrote it in pascal .
    All the instructions are included at the beginning..
    I think my approach to reordering the twiddle factors is innovative as it is written in assembly. and will work on any Intel processor..
    Code (Text):
    {This program does a Fast Fourier Transform , Decimation in Frequency ..
    it works with a data set of powers of two..It can work in real time ,
     but the data window has to be in a power of 2..It has been a while
    since i looked at it and ten years since i wrote it..
    Changing N , changes the number of samples , puting in 1024 for
    example will give you 512 unique output points ,in procedure scope you
    could uncomment the mirror image for a textbook output..this program just
    works with data ..ie no graphics..
    procedure sinewave creates a file of four different frequencies 63 , 31 ,
    15 & 7  for testing purposes..
    procedure squarewave doesnt use a disk file..
    Bit Reverse or Reorder is my own idea it is an assembly program to
    reorder the twiddle factors and you must assemble it into an object file before compiling in pascal ..and have it in the same directory as FFT..
    if you have any questions i'll try to answer them willibur @comcast .net  ..
    Reference : Digital Signal Processing Design Chapter 5 Spectral Analysis and the FFT
    sorry i didnt write down the author..}

    Program FFT(Input,Output,Indata,Outdata);
    Const
       N = 256;
       Max= 2* N-1;
    Type
       signal= array[0..Max] Of Real;
    Var
      X:Signal;
      Outdata:Text;

    {----------------------------  Square Wave --------------------------------}

    Procedure SquareWave;  { Input Array }
    Var
      I:Integer;
    Begin
         For I:=0 To N Div 2-1 Do
          Begin
            X[2*I]    := 1.0;
            X[2*I+1]  := 0.0;
          End;

         For I:=N Div 2 To N-1 Do
          Begin
            X[2*I]:= 0;
            X[2*I+1]:= 0;
          End;
    End;

    {-----------------------------  Sine Wave  -------------------------------}
    Procedure SineWave;
    Var
    Y,A:Real;
    I:Integer;
    Begin
        Assign(Outdata,'Scope1.dat');
        Rewrite(outdata);
        For I:= 0 to N-1 Do
          Begin
           A:= 2*Pi*I/N;

           X[2*I]:= sin(63 *2*Pi*I/N)+
                    sin(31 *2*Pi*I/N)+
                    sin(15 *2*Pi*I/N)+
                    sin(7  *2*Pi*I/N);

            X[2*I+1]:=0;
            Y:= 20*X[2*I];
            Writeln(Outdata,Round(Y));
          End;
          close(outdata);
    End;




    {-----------------------  Butterfly  -------------------------------------}
    Procedure Butterfly2f(Var x:Signal; p,q:Integer; wr,wi:Real);
    { Radix 2 DIF butterfly }
    Var
    tr,ti:Real;
    Begin
    tr:=       x[2*p]  - x[2*q];
    ti:=       x[2*p+1]- x[2*q+1];
    x[2*p]:=   x[2*p]  + x[2*q];
    X[2*p+1]:= x[2*p+1]+ x[2*q+1];
    x[2*q]:=   tr*wr - ti*wi;
    x[2*q+1]:= tr*wi + ti*wr;
    End;

    {--------------------------  DIF FFT  ------------------------------------}
    Procedure FFTdif(Var X:Signal; N:Integer);
    {  Radix 2 DIF FFT }
    Var
    P,q,n1,n2,j:Integer;
    wr,wi:Real;
    Begin
      n2:= N;
      Repeat
          n1:=n2;
          n2:= n2 div 2;
          j:=0;
          repeat
              wr:=  Cos(2*Pi*j/n1);
              wi:= -Sin(2*Pi*j/n1);
              p:=j;
              Repeat
                  q:= p + n2;
                  butterfly2f(x,p,q,wr,wi);
                  p:=p + n1;
              until p >= n;
              j := j + 1;
          until j = n2;
      until n2=1;
    End;

    {-------------------------------  Bit Reverse  ------------------------------}

    {$L Reorder}
    Function Reorder(I,P:word): word; external;
    { external assembly function reorder.obj }

    {-------------------------------  Screen  ------------------------------------}

    Procedure Screen;
    Var
    I,J:Integer;
     Begin
       Writeln;

       Writeln('   N      Xout(r)       Xout(i)       Magnitude ');
       Writeln;
       For J:=0 To N-1 Do
         Begin
           I:=Reorder(J,N);
           Write(J:4{,Xin[2*I]:14:6,Xin[2*I+1]:14:6});
           Write(X[2*I]:14:6,X[2*I+1]:14:6);
           Writeln(Sqrt(Sqr(X[2*I]) + Sqr(X[2*I+1]) ):14:6);
         End;
     End;

    {-------------------------------  Scope disk file  -------------------------------}
    Procedure Scope;
    Var
    I,J:Integer;
    Magnitude:Real;
    Begin
       Assign(Outdata,'Scope0.dat');
       Rewrite(Outdata);

    {   For J:= N Div 2 To N-1 Do
         Begin
           I:=Reorder(J,N);
           Magnitude:= Sqrt(Sqr(X[2*I]) + Sqr(X[2*I+1]));
           Writeln(Outdata,Round(Magnitude));
         End;
    }
       For J:= 0 To N Div 2 -1  Do
         Begin
           I:=Reorder(J,N);
           Magnitude:= Sqrt(Sqr(X[2*I]) + Sqr(X[2*I+1]));
           Writeln(Outdata,Round(Magnitude));
         End;
       Close(Outdata);
    End;


    {-------------------------------  FFT MAIN  ---------------------------------}
    Begin
         { SquareWave; }
         SineWave;
         FFTdif(X,N);
        { Screen; }
         Scope;
    End.
     
    Code (Text):

     

    Code:
    CODE    Segment Word Public
            Assume CS:Code
            Public Reorder
    Reorder Proc near
            Push BP
            Mov bp,sp
            Mov ax,ss:[BP+6]  ; Index
            mov cx,ss:[BP+4]  ; Points in fft
            shr cx,1          ; divide cx by 2
            xor bx,bx         ; clear bx
    L1:
            rcr ax,1          ; rotate bit into carry flag
            rcl bx,1          ; rotate bit in carry to bx
            shr cx,1          ; shift bit in cx till zero
            jnz L1

            Mov AX,BX
            Mov SP,BP
            Pop BP
            Ret 4
    Reorder Endp
    Code    Ends
            end
     
     
    Last edited: Feb 22, 2005
  2. jcsd
  3. cronxeh

    cronxeh 1,232
    Gold Member

    good grief man.. who uses Pascal :/

    i have a matlab code laying around somewhere if anyone needs it
     
  4. re

    nice code, can you translate it into BASIC, its been long time since i used pascal. I'm working on a spectrum analyzer project. Thanks.
     
  5. SpaceTiger

    SpaceTiger 2,977
    Staff Emeritus
    Science Advisor
    Gold Member

    For anyone using fortran or C, I recommend Numerical Recipes for simple FFT programs (the books are online). For more advanced applications, go with FFTW
     
  6. chroot

    chroot 10,426
    Staff Emeritus
    Science Advisor
    Gold Member

    Even further, I recommend that anyone who needs a very efficient version of the FFT to look at the Intel application notes for the SSE/SSE2/SSE3 vector-processing instruction sets.

    - Warren
     
Know someone interested in this topic? Share this thead via email, Google+, Twitter, or Facebook

Have something to add?