The Fast Fourier Transform


A long time ago in college, I learned a lot about signal processing. A microphone produces a signal, and a speaker processes that signal back into sound. Why is CD quality sound sampled at 44 Khz ? Because human hearing ranges from 20-20 Khz and the Nyquist theorem.  Actually, older people tend to have high frequency hearing loss. A standard NTSC TV set horizontal retrace transformer emits a 15 Khz squeal that I can’t hear any more. 15Khz is less than an octave down from 20Khz so I’m not too worried: an octave is a doubling of frequency. (For more about TV and computers, see Why was the original IBM PC 4.77 Megahertz?)


 


I worked for a company called Bolt Beranek and Newman (BBN) in Cambridge Massachusetts right after graduation from MIT in 1980.


BBN has a couple historical computer science related accomplishments:



  • 1962 Performed first public demonstration of computer time-sharing.
  • 1969 Launched the ArpaNet (a forerunner of the Internet)

 


One of the projects I worked on at BBN in the early 80’s was detecting submarines. I created a computer system that would control a sound source underwater. This source was a vertical array of cylinders of compressed air, released under computer control to make an underwater sound. The array was steerable: firing the higher ones just a little before the lower ones caused the wavefront to be steered downward.


I also created a “target simulator”, which basically consisted of an underwater microphone and speaker. My computer system would listen to the microphone for the air gun sounds, and then echo back via the speaker as if it were a submarine (at various aspects: head-on, broadside, etc). Then it would echo back the same signal 3 times: each time being 10 times louder. Thus the number of times the “submarine” was detected indicated by how much we missed.


We’d take the machines out to sea for weeks, doing sea trials, sometimes with real submarines (which are expensive to rent!). Each computer system was a PDP-11 running RT-11. My software was written mostly in Fortran, with some Macro-11 assembly language for the crucial real time processing. I also was very seasick, and never got better til well after landfall.


 


You’ve probably seen “graphic equalizers” that display the spectrum of sound coming from a music system. It’s just a bar chart that shows the bars dancing with the music. The bars range from bass frequencies on the left to treble on the right. It’s basically a spectrum of the signal currently being played. Another way to think of a spectrum is each frequency corresponds to the keys on a piano. Each “bar” could be a piano key (or range of keys). Imagine watching and listening to a piano player. The keys being pushed comprise the “spectrum” of the music.


 


A pure tone is a sine wave. This “time domain” signal can be represented as a graph with the X axis representing time. The Fourier Transform converts the signal to a “frequency domain” signal, with the X axis now representing frequency. Thus, the Fourier Transform of a simple sine wave is very simple: it consists of just a single spike at the frequency of the sine wave (a single “piano key” (minus harmonics))


 


Below is a sample of code you can run. It will allow you to choose between a few kinds of signals (2 Sine waves, Square, Impulse), display them in red, display their FFT in blue and the difference between the FFT/Inverse FFT in green. A slider and spinner allow you to vary the 2 Sine waves frequency and phase shift (theta). Theta is initially zero, and the 2 frequencies are the same initially. Thus the 2 sine waves initially appear to be one. Try varying the slider and spinner to get different signal shapes.


 


(The green signal is the difference between the original signal and the signal after being FFT’d and inverse FFT’d. Theoretically, there should be no difference, except for mathematical rounding error.)


 


I have other sample code that filters sounds: it reads a WAV file, does an FFT, filters, and inverse FFT, writes the WAV file and plays back the sound. The original WAV is my voice saying “Foxpro Rocks” with a very loud high frequency sine wave in the background. The filter removes that sine wave to recover the original voice.


Filters are quite simple. Here’s the filter to remove the loud sine wave


 


nMid=3044


nRange=150


FOR i = 1 TO np2/2


      IF BETWEEN(i,nMid-nRange, nMid+nRange)


            aSamples[i]=0


            aIm[i]=0


            aSamples[np2-i]=0


            aIm[np2-i]=0


      ENDIF


ENDFOR


 


 


BTW, a 2 dimensional FFT can be used to sharpen the focus of a photograph or for image recognition.


 


 


 


CLEAR


CLEAR ALL


*see http://www.relisoft.com/Science/Physics/fft.html and http://www.spectrumsdi.com/ch12.pdf


 


PUBLIC ox


ox=CREATEOBJECT(“signalform”)


ox.show


 


DEFINE CLASS signalform as form


      width=2^10  && must be power of 2


      height=600


      left=300


      backcolor=0xffffff


      allowoutput=.f.


      ADD OBJECT spn as spinner WITH left=130,;


            Width=60,;


            SelectOnEntry=.t.,;


            SpinnerLowValue=0,;


            SpinnerHighValue=360,;


            Increment=10,;


            value=0


      ADD OBJECT cboInput as Combobox WITH left=200


      PROCEDURE init


            this.AddObject(“oSlider”,”MySlider”)


            WITH this.oSlider


                  .min=-100


                  .max=100


                  .SmallChange=5    && a penny


                  .Value=0


                  .visible=1


            ENDWITH


            WITH this.cboInput


                  .AddItem(“2 Sine Waves”)


                  .AddItem(“Square Waves”)


                  .AddItem(“Impulse”)


                  .AddItem(“DC Input”)               


                  .Value=1


            ENDWITH


      PROCEDURE keypress(nKeyCode, nShiftAltCtrl)


            DO case


            CASE nKeyCode=27  && escape


                  thisform.release


            ENDCASE


      PROCEDURE activate


            thisform.DrawWaves


      PROCEDURE DrawWaves


            nTheta= thisform.spn.Value*PI()/180


            nFactor=thisform.oSlider.Value/100


            thisform.Cls


            thisform.Caption=TRANSFORM(nFactor)+” “+TRANSFORM(nTheta)


            yMax=thisform.Height/4


            yOffset=thisform.Height/2     && x axis halfway down


            thisform.Line(0,yOffset,thisform.Width,yOffset) && draw x axis


            xMax=thisform.Width     && must be power of 2


            DIMENSION aInput[xMax],aOutput[xMax]


            aOutput=0


            this.ForeColor=RGB(255,0,0)


            x0=0


            y0=yOffset


            cMode=thisform.cboInput.Text


            FOR x=1 TO xMax


                  DO CASE


                  CASE cMode=”2 Sine Waves”


                        aInput[x]=yMax*   (SIN(2*PI()*(1+nFactor)*x/100) + SIN(2*PI()*(1+nFactor*10.5)*x/100 + nTheta))


                  CASE cMode=”Square Waves”


                        aInput[x]=yMax * SIGN(SIN(2*PI()*(1+nFactor)*x/100) + SIN(2*PI()*(1+nFactor*10.5)*x/100 + nTheta))


                  CASE cMode=”DC Input”


                        aInput[x]=10


                  CASE cMode=”Impulse”


                        aInput[x]=IIF(x<100,-100,0)


                  ENDCASE


                  y1=aInput[x]+yOffset


                  thisform.Line(x0,y0,x,y1)


                  x0=x


                  y0=y1


            ENDFOR


            ACOPY(aInput,a2)


            FFT(@aInput,@aOutput,.f.)


            this.ForeColor=RGB(0,0,255)


            x0=0


            y0=yOffset


            FOR x =1 TO xMax/2


                  y1=yOffset-SQRT(aInput[x]^2+aOutput[x]^2)/32


                  thisform.Line(x0,y0,x,y1)


                  x0=x


                  y0=y1


            ENDFOR


            FFT(@aInput,@aOutput,.t.) && now calculate inverse FFT, display difference in green


            this.ForeColor=RGB(0,255,0)


            x0=0


            y0=yOffset


            FOR x =1 TO xMax


                  y1=aOutput[x]+yOffset


                  thisform.Line(x0,y0,x,y1)


                  x0=x


                  y0=y1


            ENDFOR


      PROCEDURE cboInput.Valid


            thisform.spn.value=0    && reset phase and freq factor to 0


            thisform.oSlider.value=0


            thisform.DrawWaves


      PROCEDURE spn.InteractiveChange


            thisform.DrawWaves


ENDDEFINE


 


DEFINE CLASS MySlider as olecontrol


      oleclass=”MSComctllib.slider.2″


      PROCEDURE Change


            thisform.DrawWaves


      PROCEDURE keypress(nKeyCode, nShiftAltCtrl)


            thisform.keypress(nKeyCode, nShiftAltCtrl)


ENDDEFINE


 


PROCEDURE FFT(aRe,aIm,fInverse as Logical )     && Fast Fourier Transform


      n=ALEN(aRe)


      nlg2 = INT(LOG(n)/LOG(2))


      n2 = n/2


      j=n2+1


      IF fInverse


            FOR i = 1 TO n


                  aIm[i]=-aIm[i]


            ENDFOR


      ENDIF


      FOR i = 2 TO n-2  && Bit Reversal order


            IF i<j


                  tr=aRe[j]


                  ti=aIm[j]


                  aRe[j]=aRe[i]


                  aIm[j]=aIm[i]


                  aRe[i]=tr


                  aIm[i]=ti


            ENDIF


            k=n2


            DO WHILE k<j


                  j=j-k


                  k=k/2


            ENDDO


            j=j+k


      ENDFOR


      le2=1


      FOR lp = 1 TO nlg2


            le=2*le2


            ur=1


            ui=0


            sr=COS(PI()/le2)


            si=-SIN(PI()/le2)


            FOR j = 1 TO le2  && each sub DFT


                  FOR i = j TO n STEP le  && butterfly loop


                        ip=i+le2


                        tr=aRe[ip]*ur – aIm[ip]*ui


                        ti=aRe[ip]*ui + aIm[ip]*ur


                        aRe[ip] = aRe[i] – tr


                        aIm[ip] = aIm[i] – ti


                        aRe[i] = aRe[i] + tr


                        aIm[i] = aIm[i] + ti


                  ENDFOR


                  tr = ur


                  ur = tr * sr – ui * si


                  ui = tr * si + ui * sr


            ENDFOR


            le2=le2*2


      ENDFOR


      IF fInverse


            FOR i = 1 TO n


                  aRe[i]=aRe[i]/n


                  aIm[i]=-aIm[i]/n


            ENDFOR


      ENDIF


RETURN


 

Comments (8)

  1. Hi Calvin,

    Always enjoy reading your blog. This reminded me of the book (and movie) "The Hunt for Red October" in which FFT technology was used and explained (somewhat) to detect enemy (Russian) submaries.

    I have running on my computer as a screen saver the SETI project which uses the FFT to detect possible intelligent signals from space. http://setiathome.ssl.berkeley.edu/

  2. I took my family to Whistler for some great skiing. Because our cell phones didn’t work too reliably…

  3. In this post about The Fast Fourier Transform , I described how simple it was to create a filter to remove

  4. bob handler says:

    Does anyone here know anything about the Temperton FFTs :

    specifically RFFTMLT and CFFTMLT ?

  5. Abey says:

    Hi, I wonder if you met Dr. John Makhoul