Physics Forums (http://www.physicsforums.com/index.php)
-   General Engineering (http://www.physicsforums.com/forumdisplay.php?f=113)
-   -   Fast Fourier Transform (http://www.physicsforums.com/showthread.php?t=64548)

 willib Feb22-05 07:58 AM

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

 cronxeh Feb25-05 09:56 PM

good grief man.. who uses Pascal :/

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

 waht Mar3-05 11:00 AM

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 Mar5-05 03:28 AM

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 Mar5-05 04:36 AM

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

 All times are GMT -5. The time now is 11:04 PM.