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 https://www.relisoft.com/Science/Physics/fft.html and https://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