How to Implement a Fast Fourier Transform Program in Pascal?

[willib]

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..

{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 doesn't 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: 
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

 

[cronxeh]

good grief man.. who uses Pascal :/

i have a MATLAB code laying around somewhere if anyone needs it

 

[waht]

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.

 

[SpaceTiger]

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

 

[chroot]

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