Physics Forums

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 02:39 PM.

Powered by vBulletin Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.
© 2014 Physics Forums