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.

• 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

WITH this.oSlider

.min=-100

.max=100

.SmallChange=5    && a penny

.Value=0

.visible=1

ENDWITH

WITH this.cboInput

.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

Tags

1. Tech Guru says:
2. Ryan Jentzsch says:

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/

3. Calvin Hsia's WebLog says:

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

4. Calvin Hsia's WebLog says:

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

5. bob handler says:

Does anyone here know anything about the Temperton FFTs :

specifically RFFTMLT and CFFTMLT ?

6. Calvin Hsia s WebLog The Fast Fourier Transform | Wood TV Stand says:

PingBack from http://woodtvstand.info/story.php?id=14072

7. Calvin Hsia s WebLog The Fast Fourier Transform | bar stools says:

PingBack from http://barstoolsite.info/story.php?id=3204

8. Abey says:

Hi, I wonder if you met Dr. John Makhoul