CS 2073 Section 2, Spring 2009
Assignment 6: Arrays: Signal Processing with the Fourier Transform
piano.c

Your assignment is to write a program that will determine, given a sample of sound from a piano key, what key (1-88) and note (A-G#) it is.

Background

We often encounter data collected over time at specific intervals. For instance, the temperature each day of this year, or the number of bacteria in a dish each minute. Sometimes this data is periodic, i.e., its behavior repeats at specific intervals. An everyday example of periodic data is sound. Sound is just a signal propogating through the air in compression waves. Think of a piano key being struck; this causes the piano string to move back and forth very rapidly, causing compression in the surrounding air which moves outward as a sound wave. Your eardrum is vibrated by this back-and-forth force, and you hear the sound. If we plot the strength (or amplitude) of the compression against time, we'll get this:

Wave Plot

The x axis is the number of seconds from the time the key was struck; as you can see, we are considering only 1/20 of a second of time. This piano key happens to be the A above middle C, tuned at 440 Hertz (abbreviated Hz). That means the string vibrates 440 times per second. The y axis plots the amplitude of the compression of air around a microphone connected to my computer; you can think of it as a voltage level. As you can see from the plot, the voltage level, or signal, spikes once every 1/440 of a second.

Sound is a continous signal, so how do we reprent it inside the computer? One way is to take a sample of the amplitude of the signal once every, say, 1/8000 of a second. We can store these samples in a large array. Whenever we want to look at the signal at time t, where t is in seconds, we just multiply t by 8000, find the closest integer, and look in the array at that point. In fact, this is exactly the way sound is stored on an ordinary CD, with samples taken every 1/44,000 of a second.

How do you determine the frequency of a sampled sound? This situation is complicated by the fact that in most situations, including the piano key strike, there are many frequencies superimposed on each other. For example, the A above also vibrates at 880 Hz, 1,320 Hz, and other frequencies, although these show up at a lower amplitude. These harmonic frequencies are what give the piano, other instruments, your voice, and pretty much every other sound its characteristic texture. However, we are now only interested in determining the fundamental frequency.

The data are given to us in the time domain, i.e., we know, for any given time, what the amplitude will be. It would be nicer if the data were in the frequency domain, where we would know, for any given frequency, what the amplitude will be. Then we can just look for the frequency with the highest amplitude, and that's the fundamental frequency (usually; it works for the middle part of the piano, at least). The Fourier Transform of a function takes a function y(t) in the time domain and ``transforms'' it into a function Y(f) in the frequency domain. It works for any periodic continuous function. The Fourier Transform, for your amusement, looks like this:

Equation for Fourier Transform

If you don't know what this means, don't worry; you don't have to.

How to do the assignment

On the system in ~dj/cs2073, I have placed twenty data files called 1.dat through 20.dat. These files contain piano sounds samples at 8000 Hz. Each file contains 16384 samples; you can read them as double values using scanf in C. You should read the values from the files by redirecting standard input from them on the Unix command prompt, e.g.: 
piano < 1.dat
Your program should read the number into an array of doubles.

You can use a special version of the Fourier Transform called the Discrete Real Fourier Transform, or DFT, on the sampled data inside the computer. You give the DFT an array of n time-domain data, sampled at m samples per second, and it gives you back an array of n frequency-domain data, where each element is the amplitude of the sound at each frequency. After you take the DFT of the sample, you search the array for the point with the maximum value and convert the corresponding array index to Hz by multiplying by m/n. You should only search up to about the 2000th position in the array; anything past that is probably noise.

Once you have determined the fundamental frequency of the piano sound, you need to convert it to a number from 1-88, representing the number of the key on the piano. A formula for the conversion from a frequency f to a piano key is:

key = 1 + the closest integer to (the log base 21/12) of (f / 27.5))

Remember, to find the logarithm base b of x, you can just divide the natural logarithm of x by the natural log of b. Also, you don't necessarily get the closest integer to a double by casting to int, e.g., 3.6 is closer to 4 than it is to 3. You should write a function called closest_int that returns the closest integer to its double parameter. The log and pow functions in math.h will find the natural logarithm and number to a power, respectively. The 'key' function works since the frequency of a key doubles every twelfth note (hence taking the log base twelfth root of two) and the first note on the piano vibrates at 27.5 Hz.

Once you have the key, the final task is to decide which note it is. There are twelve symbols for notes, repeating every twelve notes. Your program should print the notes as A A# B C C# D D# E F F# G G#. (The '#' is read ``sharp''). The first note on the keyboard is an A; so is the 13th, the 25th, and so forth. Use the % (modulus) operator to help decide which key is which; don't use 88 if statements.

Using the FFT library function

To do the Fourier transform, you will use a library provided by the professor. This is similar to the standard C library or the math library, but the professor has provided it in a different directory so you will have to tell the C compiler how to find it.

The library file libfft.a on the main machines in /home/dj/cs2073 contains a precompiled version of the Discrete Real Fourier Transform as a C void function. You should declare it in your program after the #includes and before main like this: 

void rfft (double [], int);
You can apply this FFT to your double array of length 16384 using the following C function call: 
        rfft (array, 16384)
where array is the name of your array. Your time-domain data will be replaced with the frequency domain data. Compile your program with: 
cc -o piano piano.c -L/home/dj/cs2073 -lfft -lm
The output of your program should be give the maximum amplitude, the fundamental frequency, the key number, and the symbol for the note. So, your program should: Here is the output of my program run on the first three data files, 1.dat, 2.dat, and 3.dat: 
% ./piano < 1.dat
maximum amplitude is  4097.06 at 196.29 Hz, key = 35
key letter is G
% ./piano < 2.dat
maximum amplitude is  6109.61 at 622.56 Hz, key = 55
key letter is D#
% ./piano < 3.dat
maximum amplitude is  2902.94 at 154.79 Hz, key = 31
key letter is D#
Don't worry about the amplitudes matching exactly; the important thing is that the frequencies and key numbers match. To turn in this program, do exactly this from the prompt on the main machines, after having compiled your program: 
% cat piano.c > output
% ./piano < ~dj/cs2073/1.dat >> output
% ./piano < ~dj/cs2073/2.dat >> output
% ./piano < ~dj/cs2073/3.dat >> output
% ./piano < ~dj/cs2073/4.dat >> output
% ./piano < ~dj/cs2073/5.dat >> output
% ./piano < ~dj/cs2073/6.dat >> output
% ./piano < ~dj/cs2073/7.dat >> output
% ./piano < ~dj/cs2073/8.dat >> output
% ./piano < ~dj/cs2073/9.dat >> output
% ./piano < ~dj/cs2073/10.dat >> output
% ./piano < ~dj/cs2073/11.dat >> output
% ./piano < ~dj/cs2073/12.dat >> output
% ./piano < ~dj/cs2073/13.dat >> output
% ./piano < ~dj/cs2073/14.dat >> output
% ./piano < ~dj/cs2073/15.dat >> output
% ./piano < ~dj/cs2073/16.dat >> output
% ./piano < ~dj/cs2073/17.dat >> output
% ./piano < ~dj/cs2073/18.dat >> output
% ./piano < ~dj/cs2073/19.dat >> output
% ./piano < ~dj/cs2073/20.dat >> output
This will run your program on all the files and put the result into a file called output. Make sure to use the double "greater than" signs (i.e. >>) so that the output of the program will be appended to the output file. After you have done this, send the file called output to the teaching assistant.

This program can be done in about 70 lines of C code, not including comments. It is not as difficult as it sounds, but it is by no means easy, so you should get started right away.

This assignment is due Tuesday, April 28, 2009 at midnight.