Fast Fourier Transform

In summary: End; {--------------------------------------------- DIF FFT ---------------------------------------------} In summary, 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.
  • #1
willib
227
0
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 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: 
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:
Engineering news on Phys.org
  • #2
good grief man.. who uses Pascal :/

i have a MATLAB code laying around somewhere if anyone needs it
 
  • #3
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.
 
  • #4
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
 
  • #5
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
 

What is Fast Fourier Transform (FFT)?

Fast Fourier Transform (FFT) is an algorithm used to efficiently compute the discrete Fourier transform (DFT) of a sequence or signal. It converts a signal from its original representation in the time or spatial domain to a representation in the frequency domain.

What are the advantages of using FFT?

FFT is significantly faster than the traditional DFT algorithm, making it ideal for real-time applications. It also requires less computational resources, making it more efficient for processing large data sets. Finally, FFT can be easily implemented on modern computer architectures, making it a popular choice for many applications.

How does FFT work?

FFT works by breaking down a signal into smaller, simpler components and computing the DFT of each component. Then, the results are combined to produce the final DFT of the entire signal. This is achieved through a process of recursive decomposition and recombination, making FFT much more efficient than the traditional DFT algorithm.

What are some common applications of FFT?

FFT has a wide range of applications in various fields such as signal processing, image processing, data compression, and audio and video encoding. It is also used in scientific computing, particularly in solving partial differential equations and analyzing complex systems.

Are there any limitations to using FFT?

While FFT is a powerful algorithm, it does have some limitations. It can only be applied to signals that are discrete and periodic. It also assumes that the signal is stationary, meaning that it does not change over time. Additionally, FFT is not suitable for signals with rapid changes, as this can lead to errors in the computed DFT.

Similar threads

  • Linear and Abstract Algebra
Replies
11
Views
1K
Replies
1
Views
759
  • General Math
Replies
7
Views
1K
Replies
4
Views
212
  • Advanced Physics Homework Help
Replies
8
Views
622
  • Engineering and Comp Sci Homework Help
Replies
6
Views
2K
Replies
27
Views
742
  • General Engineering
Replies
2
Views
2K
  • Programming and Computer Science
Replies
1
Views
916
Replies
8
Views
2K
Back
Top