Fun with Fourier

After the success of generating sound directly in Sniff, I guess I had some signal processing ideas in my head (Also I taught a class on sampling theory yesterday so maybe that had something to do with it!). Recording audio is pretty easy - just hook a mic up to an analog input and you can measure it sufficiently fast to fill up an Arduino's memory in no time at all. We could easily record a snapshot of audio and draw it on screen, but to see what's really going on plotting amplitude against time is only of limited use.

To get a good look inside a piece of audio we need to plot amplitude against frequency. Almost all sounds are a mix of low sounds and high sounds played at the same time, and drawing them this way lets us see how a particular complex sound is made up from these simpler components.

At this point a whole load of complex maths gets involved, but the main thing is that we can convert our simple amplitude/time recording to a amplitude/frequency version using a Fourier Transform. There's a special version of this that we normally use when we implement it on a computer called the Fast Fourier Transform or FFT for short. That's not to say that the FFT is fast - it's still quite a lot of work for a small computer, but its the fastest way to do it.

There are a bunch of FFT libraries which you can download (including some good ones for Arduino), because you really don't want to write your own - its tricky, and to get good performance you need to add in a lot of tricks (check out FFTW for a serious implementation). There's also the issue that if we're running on Arduino we don't have a lot of memory. To make the FFT run quickly they use a lot of precomputed tables which are going to eat up memory.

Then I stumbled across the Goertzel algorithm - its slower than an FFT if you can't to compute the full transform, but has the neat property that it calculates one value at a time - if we're plotting the graph we can calculate each point independently rather than having to calculate all the points. It's also light on trig functions, so we don't need any pre-computed tables. In other words its slower, but doesn't use any significant memory other than the source data. In fact if we only want to measure a particular frequency (how much 440Hz is in this audio) then we don't even need to store the incoming signal!

As an extra bonus its much simpler that an FFT, and is non-recursive.

Step 1 is to make a signal for us to process:

make x list of number

.repeat 32
..add 1 to x
.repeat 32
..add 0 to x

This fills x with a 64 samples of a square wave. Next we can use the Goertzel code to show how much energy is in the signal a frequencies from 0 (any constant signal) 1 (1 wave is 64 samples) up to 32 (32 cycles within the 64 samples we have).

.set k to 0
.repeat 33
..broadcast Goertzel and wait
..say join join [ k ] ":" [ power ]
..change k by 1

The actually Goertzel code is remarkably simple. How it works on the other hand isn't:

when Goertzel
.set coeff to 2 * cos of (360*k/64)
.set sp to 0
.set sp2 to 0
.set n to 1
.repeat 64
..set s to (item n of x) + coeff*sp -sp2
..set sp2 to sp
..set sp to s
..change n by 1
.set power to (sp2*sp2+sp*sp-coeff*sp*sp2)/32

It basically goes through the list of samples calculating a value s which at each step by summing the value of that sample, and the value of s at the previous two steps (sp and sp2). The sum is weighted by coeff which is dependant on the frequency band we're interested in.

If we run this we get the results:

0.000:32.000

1.000:12.980

2.000:0.000

3.000:1.451

4.000:0.000

5.000:0.529

6.000:0.000

7.000:0.275

8.000:0.000

9.000:0.171

10.000:0.000

11.000:0.118

12.000:0.000

13.000:0.088

14.000:0.000

15.000:0.069

16.000:0.000

17.000:0.057

18.000:0.000

19.000:0.048

20.000:0.000

21.000:0.042

22.000:0.000

23.000:0.038

24.000:0.000

25.000:0.035

26.000:0.000

27.000:0.033

28.000:0.000

29.000:0.032

30.000:0.000

31.000:0.031

32.000:0.000

You'll see that the 0 has a value of 32 - the sum of all the samples.

The value of 1 is large because our square wave looks a lot like a sine wave of period 64. From there we see that all of the even values are 0, and the odd values decrease progressively -the frequency signature of a square wave!

Here's the complete code:

make k number
make n number
make x list of number

make s number
make sp number
make sp2 number
make coeff number
make power number

when Goertzel
.set coeff to 2 * cos of (360*k/64)
.set sp to 0
.set sp2 to 0
.set n to 1
.repeat 64
..set s to (item n of x) + coeff*sp -sp2
..set sp2 to sp
..set sp to s
..change n by 1
.set power to (sp2*sp2+sp*sp-coeff*sp*sp2)/32

make t number
when start
.repeat 32
..add 1 to x
.repeat 32
..add 0 to x
.
.set t to timer
.set k to 0
.repeat 33
..broadcast Goertzel and wait
..say join join [ k ] ":" [ power ]
..change k by 1
.say join "time:" [timer - t]

It's pretty slow - on an Uno its around 0.1 secs, excluding the printing, which is far slower than an FFT, but there's a big win in terms of simplicity, flexibility memory footprint.

But here's the fun bit... Scratch and Sniff are essentially the same language: If we can do it in Sniff we can do it in Scratch!

http://scratch.mit.edu/projects/19443676/

This is identical code - I simply dragged out the Scratch blocks, copying the Sniff code. Though in this case we're going the opposite direction to the way we expect students to (moving from Sniff to Scratch rather from Scratch to Sniff) the experience confirmed a couple of core Sniff ideas:

• We were able to move code directly between Scratch and Sniff
• Knowledge of one made working the other easier
• Building code in text form required more concentration
• Building complex code graphically was slower

It's also a great endorsement of Scratch - we solved a University level problem in a Kindergarten tool! It shows that as a language Scratch is sufficiently powerful to express complex ideas, though for experienced programmers Sniff is probably a better way to do Scratch!