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