Fast Fourier Transform

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:
{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:
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:

cronxeh
Gold Member
good grief man.. who uses Pascal :/

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

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
Staff Emeritus