E12 Physical Systems Analysis: Lab #2 Fighting Corrupt Speech

Graphs and results from parts 1-8 due 1 week after lab date

Final report due 2 weeks after lab date

Overall Objectives

This lab will give you some hands-on experience with the Fourier transform, which is one of the most useful tools in signal processing and analysis. We will be working with 1D Fourier Transforms to clean up noisy speech.

Introduction

In order to process both speech and images on a computer, they must first be sampled at discrete intervals. The sampling frequency--the inverse of the sampling interval--is determined by the characteristics of the signal being measured.

Human speech, for example, uses a range of frequencies from approximately 200Hz to 10kHz, although the most important range is below 5kHz. The Nyquist criterion tells us that we must sample a signal at a frequency at least twice that of the frequency of the signal we want to measure. Otherwise, aliasing results, which is the signal folding back on itself in the frequency domain. Because of the Nyquist criterion, speech is usually sampled at no less than 8kHz. High quality CD music is sampled at 44kHz, which captures pretty much the whole range of human hearing abilities. For this lab, we will sample the speech at 11.025kHz.

We would like to execute a Fourier transform [FT] on our data in order to analyze and process it. Because our data is a string of samples, however, we cannot treat it as a continuous signal when executing a Fourier analysis. Instead, we must use the discrete Fourier transform [DFT]. The DFT executes the Fourier Transform of a finite-length sequence of complex points.


Discrete Fourier Transform definition for N samples

For real-valued input data like speech samples, the value of the complex part of each sample is zero.

The most important difference between the DFT and the continuous FT is that the DFT works on a finite set of samples. From the point of view of the DFT, however, the finite set of samples represents one "example" of an infinitely repeated sequence of samples. In other words, the signal you are actually analyzing is the signal you would get by repeating the sequence of samples, with x(0) coming after x(N-1), infinitely in both directions (positive and negative).

Likewise, DFT(-k) = DFT(N-k), which means the negative frequency portion of the spectrum appears in the top half of DFT(k) in locations k = N/2 to k = N-1. Normally, only the positive frequences are plotted when viewing the spectrum.

A second major difference between the DFT and the continuous FT is that, in order to extract a finite sample to analyse, we have to cut out a window of our signal of interest. For example, we may want to analyze a half second of speech sampled at 11kHz from a sample that is 3s long. In order to extract the sample we are, in effect, multiplying the 3s sample by a rectangular window function that has the value 1 within the half-second period of interest and 0 otherwise.

This doesn't seem like a real big deal until you realize that multiplication in the time domain is equal to convolution in the frequency domain. This means that the spectra of the windowing function and the signal of interest will be convolved--in other words pretty badly mixed up.

Now there is a function that acts as the identity in convolution: the delta function. If you convolve a delta function centered at 0 with another function g(n) you get back g(n). Since the inverse FT of the delta function at 0 is a uniform value over all time (a DC signal), this is the same as using a window of infinite length. This should make sense, since a window of infinite length basically gets you back to the continuous FT, which doesn't have to use windows.

So where is this going? Well, what we want to find is a windowing function that looks like a delta function in the Fourier domain so that it doesn't mess up our signal analysis too much. A square window, unfortunately, looks like a large lobe in the center and many smaller lobes to the side. While the center lobe is fairly thin and peaked--similar to a delta function--the side lobes mean that the convolution will "leak" information from one frequency to many others.

The signal processing community has studied this problem for a long time and there are several windowing functions that are in common usage. One of the most commonly used is the Hamming function, which looks like a Gaussian that goes towards zero at the edges of the window. The Hamming window (see definition below) has a spectrum that is a single lobe, 46 percent wider than a rectangular window, but with virtually no side lobes. Thus, the effect of the Hamming window is to blur the Fourier spectrum a bit, but without the wider "leakage" that occurs with a rectangular window. This is the price you pay for the uncertainty that comes with executing an FT on a finite sequence of samples.


Hamming Window Function for N samples

One of the major uses of the Fourier Transform in signal analysis is to either detect or remove repeating signals. We will be doing the latter in this lab and the former in a follow-up lab later on.

In particular, for this lab we will be trying to remove a 400Hz signal from a captured speech signal. In the time domain, the captured signal can be represented as y(t) = s(t) + n(t), where s(t) is the speech and n(t) is the noise. In the Fourier domain, the signal is Y(w) = S(w) + N(w). While it may seem like we haven't gained much in the process, in the Fourier domain the periodic noise signal and its harmonics are compressed into a few points representing particular frequencies. If we reduce the value of these points in the Fourier domain and then take the inverse DFT, we should get a much cleaner signal.

Now to work.

Set Up

For this lab you will need the following items:

We will be using matlab for this exercise, so sharpen up your m-file skills.

Tasks

  1. Start up the sound recording tool (Start -> Programs -> Accessories -> Entertainment -> Sound Recorder). Under the Edit Menu, select Audio Properties. In the window that comes up, click on Volume under the microphone symbol. Then set the volume for Recording and for the Microphone to about half the maximum. Record a sample of you talking and make sure the volume is reasonable with no distortion. Use the Properties under the File menu in the Sound Recorder to set the recording quality and sampling rate. Use a 11.025kHz sampling rate with 8-bit mono.

  2. Record one person in your lab group saying a 3s phrase without any intentional noise. Listen to the captured phrase and make sure it is not clipped or distorted. (Clipping is hitting the boundary of the dynamic capability of the microphone.) Save this speech recording in your working directory and call it cleanphrase.

  3. Now hook up the speaker to the function generator and set the function generator to 400Hz on a reasonably loud amplitude. Place the speaker directly in front of the microphone and record a 3s sequence of sound. Save the speech recording as noise. Leave the speaker in place when you are finished.

  4. Now record yourself saying the same phrase as before, but with the speaker turned on and still placed directly in front of the microphone. Speak quietly so that the noise is at a similar level as your speech in the recording. Save this speech recording in your working directory and call it corrupt.

  5. Time to move into Matlab. Start up Matlab and set the working directory to the place you saved the speech recordings. Load each recording into Matlab using the WAVREAD() function. For example, to load the noise sound file into a variable called noise in matlab, use: noise = wavread('noise') function. If you want to see the time waveforms, plot the input variables, although you may only want to plot a thousand points or so if you really want to see any detail. Use plot(noise(1:1000)); to plot the first thousand points of the noise signal.

  6. Now we want to take the FT of each of the speech sequences. However, we need to have the same number of samples in the noise sample and the corrupted speech sample so that their FTs are directly comparable. Do this by taking the FT of only the first N samples, where N is the number of samples in the corrupt data (your noise sample should be longer than this). You can discover the size of a variable using the size() function. To identify the size of the corrupt speech sample and save the number of points in the variable N, use:

    [N, trash] = size(corrupt);

    To set the noise sample to N data points, use:

    noise = noise(1:N);

    To take the DFT of the samples use the function fft:

    noiseFFT = fft(noise);

    The first time through this process, don't worry about the Hamming window, especially if your phrases begin and end during relative silence.

  7. Look at the plots of the Fourier magnitudes with appropriate axes showing frequency. To get the magnitude of a complex number (which is what the Fourier Transform consists of) use the abs() function. In order to identify the particular frequencies you will need to create a variable to represent the frequency axis. Start by creating a variable

    DELTAF = 1 / (N * Ts)

    where N is the number of elements in the speech sample and Ts is time between samples (1 / fs). Then create a variable FREQ that goes from 0 to N*DELTAF. If DELTAF is defined as above and the sampling frequency (hopefully 11025) is defined as FS, then you can do this using the following expression:

    FREQ = 0:DELTAF:FS/2;

    FREQ = FREQ';

    The second statement is necessary to convert the FREQ matrix from a row vector to a column vector.

    Now you can plot the FTs versus FREQ and see the spectra of the three samples. For example, to plot the magnitude of the positive half of the noise spectrum, use the following:

    plot(FREQ(1:N/2), abs(noiseFFT(1:N/2)));

    You may want to make an m-file that does all of the steps so for you automatically. Look for the 400Hz signal in the noise sample and its harmonics.

  8. Write an m-file that takes as arguments the FTs of corrupt and noise and then uses noise to remove the noise in the signal. If we model the speech and noise as additive signals y(t) = s(t) + n(t), then we should be able to subtract our estimate of the noise from the speech and reduce its effect. This is generally done in the power spectrum (magnitude squared). Use the .* (element by element multiplication) and abs() operators to calculate the power spectrum as follows:

    npc = abs(noiseFFT) .* abs(noiseFFT);

    In the m-file, calculate the power spectrum and then subtract the noise from the signal. Since a power spectrum must be positive, clip any negative results to zero or invert their sign (choose one method). Use a for loop and an if statement.

    Then, using the new magnitude from the power spectrum and the phase of the corrupt phrase, recreate the complex number for each sample. The m-file should return a cleaned up signal.

    The tricky part here is recreating the Fourier transform, which is a series of complex numbers. In part one, you subtract the power spectrum of the noise from the power spectrum of the corrupted speech, which results in a new power spectrum. However, the power spectrum does not represent the phase of the signal. To get a complete cleaned up signal, you need both parts. The output signal should have the same phase as the corrupted speech signal, so you can rebuild the magnitude and phase from the new power spectrum and the old phase information.

    9. Take the inverse FT (ifft(a)) of the cleaned up spectrum and then play the real portion to see how it sounds. You can play any matrix by using the function SOUNDSC(REAL(mysound), mysamplingrate). Did this processing reduce the intensity of the noise in the signal?

  9. As a slight improvement on this method, write an m-file that blurs the noise spectrum prior to suppresssion. You can blur a signal by repeatedly convolving it with either a box filter or a Gaussian approximation (like the matrix [1 2 1]). With the blurred signal, select a threshold magnitude, like 10, and for any magnitude in the blurred noise signal greater than 10, multiply the corrupt phrase element by something like 0.001 to reduce its magnitude. Then take the IFFT of the resulting spectrum. As before, listen to the cleaned up speech to see if it made any improvement over the first method.

  10. As a final test, write an m-file that simply suppresses a particular frequency f and its harmonics by reducing the magnitude of the signal spectrum around f and its harmonics 2f, 3f, and so on. Use this technique to try and remove the 400Hz noise. Again, listen to the result to see how well it works.

Extensions

  1. Does using a Hamming window on the signal change anything or improve the results?

  2. Try other methods of removing the noise by modifying spectra. You can implement band-pass, low-pass, and high-pass filters in the Fourier domain by simply reducing the values of the frequencies you don't want.

  3. Try creating a corrupt sample where you can hardly hear the speaker and see if one of the subtraction methods will still work.

  4. Test out your m-file that suppresses a particular frequency on other frequencies of noise besides 400Hz.

  5. Try filtering out your voice and leaving just the 400Hz noise (why not?).

  6. Try executing the FTs on sub-sections of the speech rather than the whole thing. See if this produces improved results. You may want to use a Hamming window for this extension.

Lab Report

The writeup for this lab is not totally formal, but should be structured. You may turn it in either as hard copy or as a URL (web page). Please include the following sections in your writeup.

  1. Abstract: 100 words or less giving a summary of what you did.

  2. Introduction and problem definition: define the task and the overall approach.

  3. Methods and experiments: explain what you did, and what experiments you executed. Include descriptions of any extensions you made to the basic tasks.

  4. Results and discussion: give a summary of how things turned out, present graphs to show your results, and then discuss the interesting parts and answer the following questions: