Good ol' Fourier
in [DSP]
|
From: markroworth on 15 Aug 2006 21:12 Hi there, I don't have a background in DSP or a degree in Mathematics. I'm trying to implement an FFT. I understand what it does, but I don't need to understand how it works (yet). I've spent hours scouring the net and I've found either pages with complicated mathematical equations on, which are beyond me, or undocument C and Fortran programs generally littered with gotos. Does anyone know of a webpage that contains some pseudocode that describes an FFT with comments and input and output variables documented? I'd be very grateful, for I am somewhat shooting in the dark. Many thanks, Mark Roworth
From: Randy Yates on 15 Aug 2006 21:33 "markroworth" <roworthm(a)yahoo.co.uk> writes: > Hi there, > > I don't have a background in DSP or a degree in Mathematics. I'm trying > to implement an FFT. I understand what it does, but I don't need to > understand how it works (yet). I've spent hours scouring the net and > I've found either pages with complicated mathematical equations on, > which are beyond me, or undocument C and Fortran programs generally > littered with gotos. > > Does anyone know of a webpage that contains some pseudocode that > describes an FFT with comments and input and output variables > documented? I'd be very grateful, for I am somewhat shooting in the > dark. Many thanks, > > Mark Roworth Hi Mark, Here's an old C routine I put together some 17 years ago. If it helps, go for it. --Randy #pragma title("Fast Fourier Transform Kernal") #pragma subtitle("Randy Yates") #pragma pagesize(58) #include<math.h> #include<stdio.h> #include<string.h> /* Programmer: R. Yates Date: 2/25/89 Function: Fast Fourier Transform Algorithm Kernal Description: This function computes either the direct or inverse FFT using decimation-in-time computational techniques. The calling routine passes a real array, an imaginary array, the power of two of the FFT (i.e., 2**M=length of array), and designates the forward or inverse FFT. */ /* Global Variables: */ int near fftk_i, near fftk_IP, near fftk_j, near fftk_k, near fftk_l, near fftk_N, near fftk_NV2, near fftk_LE, near fftk_LE1; double near fftk_pi, near fftk_Ux, near fftk_Uy, near fftk_Wx, near fftk_Wy, near fftk_Tx, near fftk_Ty, near fftk_tmpx, near fftk_tmpy, near fftk_tmp1x, near fftk_tmp1y; /* Setup exception handler for atan2 domain messages: */ int matherr(x) struct exception *x;{ if (x->type == DOMAIN && (strcmp(x->name,"atan2"))==0){ x->retval=0; return(1);} else{ return(0);} } void fftk(re,im,M,INV) double re[],im[]; int INV,M; { /* Calculate fftk_N = 2**M: */ fftk_N=pow(2,M); /* Reorder input data: */ fftk_j=1; fftk_NV2=fftk_N/2; for (fftk_i=1; fftk_i<fftk_N; fftk_i++){ if (fftk_i<fftk_j){ fftk_Tx=re[fftk_j-1]; fftk_Ty=im[fftk_j-1]; re[fftk_j-1]=re[fftk_i-1]; im[fftk_j-1]=im[fftk_i-1]; re[fftk_i-1]=fftk_Tx; im[fftk_i-1]=fftk_Ty;} fftk_k=fftk_NV2; while (fftk_k<fftk_j){ fftk_j=fftk_j-fftk_k; fftk_k=fftk_k/2;} fftk_j=fftk_j+fftk_k;} /* Calculate fftk_pi */ fftk_pi = 2.0 * asin(1.0); /* Compute the DFT */ for (fftk_l=1; fftk_l<=M; fftk_l++){ fftk_LE=pow(2.0,fftk_l); fftk_LE1=fftk_LE/2; fftk_Ux=1; fftk_Uy=0; fftk_Wx=cos(fftk_pi/fftk_LE1); fftk_Wy=INV * -sin(fftk_pi/fftk_LE1); for (fftk_j=1; fftk_j<=fftk_LE1; fftk_j++){ for (fftk_i=fftk_j; fftk_i<=fftk_N; fftk_i=fftk_i+fftk_LE){ fftk_IP=fftk_i+fftk_LE1; fftk_tmpx=re[fftk_IP-1]; /* save real part temp. */ fftk_tmpy=im[fftk_IP-1]; /* save imag part temp. */ fftk_Tx=(fftk_tmpx*fftk_Ux) - (fftk_tmpy*fftk_Uy); fftk_Ty=(fftk_tmpy*fftk_Ux) + (fftk_tmpx*fftk_Uy); re[fftk_IP-1]=re[fftk_i-1] - fftk_Tx; im[fftk_IP-1]=im[fftk_i-1] - fftk_Ty; re[fftk_i-1]=re[fftk_i-1] + fftk_Tx; im[fftk_i-1]=im[fftk_i-1] + fftk_Ty;} fftk_tmpx=fftk_Ux; fftk_Ux=(fftk_Wx*fftk_Ux) - (fftk_Wy*fftk_Uy); fftk_Uy=(fftk_Wy*fftk_tmpx) + (fftk_Wx*fftk_Uy);}} /* Normalize to 1/fftk_N if inverse DFT: */ if (INV==-1) { for (fftk_i=0; fftk_i<fftk_N; fftk_i++){ re[fftk_i]=re[fftk_i]/fftk_N; im[fftk_i]=im[fftk_i]/fftk_N;}} /* End FFT Kernel */ fftk_end:;} -- % Randy Yates % "I met someone who looks alot like you, %% Fuquay-Varina, NC % she does the things you do, %%% 919-577-9882 % but she is an IBM." %%%% <yates(a)ieee.org> % 'Yours Truly, 2095', *Time*, ELO http://home.earthlink.net/~yatescr
From: rtonyreeder on 16 Aug 2006 02:50 markroworth wrote: > Hi there, > > I don't have a background in DSP or a degree in Mathematics. I'm trying > to implement an FFT. I understand what it does, but I don't need to > understand how it works (yet). I've spent hours scouring the net and > I've found either pages with complicated mathematical equations on, > which are beyond me, or undocument C and Fortran programs generally > littered with gotos. > > Does anyone know of a webpage that contains some pseudocode that > describes an FFT with comments and input and output variables > documented? I'd be very grateful, for I am somewhat shooting in the > dark. Many thanks, > > Mark Roworth I struggled with understanding the FFT algorithm for years, a little at a time. I finally got it by writing an eight point DFT out, and then started factoring and regrouping terms. With a little work I had terms grouped in a method that allowed for quick in place calculation, working with two terms at a time. The order of the terms was in the bit reversed order, which explained that mystery. Give it a try. I was then able to write an FFT routine for any power of two size. MathCad and MatLab both eliminated the need for me to write such fundamental code. I haven't looked on the web for any of this, having learned it many years ago, prior to such. Brigham's "The Fast Fourier Transform" is an excellent book for understanding things about discrete Fourier transforms. Of particular import to me over the years is the part about doing non-periodic convolutions with the inherently periodic FFT. I'm think he also has a flow chart for an FFT. Good luck.
From: Andor on 16 Aug 2006 05:42 markroworth wrote: > Hi there, > > I don't have a background in DSP or a degree in Mathematics. I'm trying > to implement an FFT. I understand what it does, but I don't need to > understand how it works (yet). I've spent hours scouring the net and > I've found either pages with complicated mathematical equations on, > which are beyond me, or undocument C and Fortran programs generally > littered with gotos. > > Does anyone know of a webpage that contains some pseudocode that > describes an FFT with comments and input and output variables > documented? I'd be very grateful, for I am somewhat shooting in the > dark. Many thanks, Maybe a bit of both (FFT mathematics and example code) is the following recursive implementation of a power-of-two FFT in pseudo code (assuming a complex vector math functionality, with the symbol ".*" denoting point-wise multiplication of two complex vectors and "j" the imaginary unit): complex vector FFT(complex vector x, integer N) { if (N > 2) { complex vector Xe[N/2] = FFT( x[0, 2, ..., N-2], N/2 ); complex vector twids[N/2] = exp(-2 * pi * j * (k=0:N/2-1)/N ); complex vector Xo[N/2] = twids .* FFT( x[1, 3, ..., N-1], N/2 ); return complex vector X[N] = [Xe+Xo, Xe-Xo]; } else return complex vector X[2] = [ x[0] + x[1], x[0] - x[1] ]; } Does this get you anywhere? Regards, Andor
From: glen herrmannsfeldt on 16 Aug 2006 06:16
markroworth wrote: (snip) > Does anyone know of a webpage that contains some pseudocode that > describes an FFT with comments and input and output variables > documented? I'd be very grateful, for I am somewhat shooting in the > dark. Many thanks, I always like the Numerical Recipes (or Numerical Recipes in C) descriptions, which include both code and explanation of the theory behind the code. To me they have a good balance on the description, enough to understand it, but not so much detail that you can't find the important parts. I don't know that it is on the web, but you should be able to find it in an engineering library. -- glen |