CS 1713 Section 1, Summer 1997
Assignment 6: Final Program: 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 runner in /home/staff3/djimenez/cs1713, 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. Note: Your program should get the name of the file to process from the command line, i.e., argc and argv; it should process only one file at a time. Your program should read the number into an array of complex numbers represented in C as an array of complexs defined: 
typedef struct _complex {
	double	real,	/* real part */
		imag;	/* imaginary part */
} complex;
Initially, you will just set the imaginary parts equal to 0, placing the values read in into the real parts.

You can use a special version of the Fourier Transform called the Discrete 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 complex-valued frequency-domain data, where the magnitude of each element represents 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 magnitude 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 complex 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 simple formula for the conversion from a frequency f to a piano key is:

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

where 2^(1/12) means "two to the one twelfth power." 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. The files complex.o and complex.h in the same directory are object code and a header file for some complex number functions including a fast version of the DFT, called the complex Fast Fourier Transform, or FFT, as a C void function. You can apply this FFT to your complex array of length 16384 using the following C function call: 
        cfft (complex_array, 16384)
where complex_array is the name of your array. Your time-domain data will be replaced with the frequency domain data from the Fourier transform. Copy complex.o and complex.h to your directory. Call your program file piano.c. Use the following Makefile: 
CC      =       gcc
CFLAGS  =       -Wall
LDLIBS  =       -lm

all:    piano

piano:  piano.c
        $(CC) -o piano piano.c complex.o $(LDLIBS)
The complex.o object file also contains code for two functions: 
complex make_complex (double, double);
double magnitude_complex (complex);
make_complex accepts two doubles and returns a complex with the first number as the real part, and the second as the imaginary part. magnitude_complex returns the magnitude of a complex number.

You should do all your work in a subdirectory called assign5. 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: 
runner% ./piano 1.dat
maximum amplitude is  4097.06 at 196.29 Hz, key = 35
key letter is G
runner% ./piano 2.dat
maximum amplitude is  6109.61 at 622.56 Hz, key = 55
key letter is D#
runner% ./piano 3.dat
maximum amplitude is  2902.94 at 154.79 Hz, key = 31
key letter is D#
To turn in this program, do exactly this from the prompt on runner, after having compiled your program: 
% cat piano.c > output
% piano ~djimenez/cs1713/1.dat >> output
% piano ~djimenez/cs1713/2.dat >> output
% piano ~djimenez/cs1713/3.dat >> output
% piano ~djimenez/cs1713/4.dat >> output
% piano ~djimenez/cs1713/5.dat >> output
% piano ~djimenez/cs1713/6.dat >> output
% piano ~djimenez/cs1713/7.dat >> output
% piano ~djimenez/cs1713/8.dat >> output
% piano ~djimenez/cs1713/9.dat >> output
% piano ~djimenez/cs1713/10.dat >> output
% piano ~djimenez/cs1713/11.dat >> output
% piano ~djimenez/cs1713/12.dat >> output
% piano ~djimenez/cs1713/13.dat >> output
% piano ~djimenez/cs1713/14.dat >> output
% piano ~djimenez/cs1713/15.dat >> output
% piano ~djimenez/cs1713/16.dat >> output
% piano ~djimenez/cs1713/17.dat >> output
% piano ~djimenez/cs1713/18.dat >> output
% piano ~djimenez/cs1713/19.dat >> output
% piano ~djimenez/cs1713/20.dat >> output
% Mail -s "piano program" djimenez@ringer.cs.utsa.edu < output
This will run your program on all the files and send me the output. (If you have written other files beside piano.c, e-mail me them, too, separately.) Expect the run time of this program to be rather long, on the order of a few seconds; an FFT is a computationally expensive operation, especially when there are 30 other students on runner running the same program.

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 Friday, August 1, 1997 at midnight.