Filtering of SampledSoundList Objects

In this Notebook, we'll examine the use of Fourier Transforms to construct various filters for SampledSoundList objects.

Theory of Operation

Mathematica can represent sound in several formats, as lists of samples or as functions which return sample values.  Though the functional representation is smaller, using much less memory, the list representation has some key advantages that we'll make use of in this sample Notebook.

One major advantage of using the list representation is that Mathematica provides several functions for manipulating lists such as those used to hold sampled sounds.  For example, it is trivial to reverse the order of the elements of a list (to play sounds in reverse), concatenate two lists, to multiply each element of a list by a scalar (to amplify or attenuate sounds), or even to multiply each element of a list by a corresponding element of another list (to apply an amplitude envelope).  Perhaps the most useful of these, however, is the facility to perform Fast Fourier Transforms (FFTs) on lists of data.

In this Notebook, we'll use the FFT (and a few of the other methods mentioned above) to build filters for SampledSoundList objects.  In order to understand how these filters work, it is first necessary to understand the structure of the SampledSoundList object.

A SampledSoundList object is of the form

	SampledSoundList[{num, num, num, ..., num}, srate]

The first item in the square brackets is the list of sampled sound data.  When the sound is played, these values (which must all be real numbers between -1 and +1) are converted to voltages which are applied to the speaker to make the sound.  The second parameter is the sample rate, which determines how quickly the samples are presented to the hardware—that is, srate samples are converted each second to make up the sound.  Thus, the duration of the sound is determined by the length of the list of samples and by the sample rate.

A SampledSoundList object is played using Show[Sound[sound]].

To build a filter for a sound, we can use the FFT algorithm to obtain the frequency-domain representation of a sound, adjust the frequency response to our liking, and convert the signal back to the time domain using the inverse FFT.

The Fourier function simply takes a list of data and performs an FFT on it, returning a new list of equal length that contains information about the frequency spectrum of the original.  Each element of the returned list is a complex number, indicating magnitude and phase of a specific frequency component of the original list.  We can use InverseFourier on the resultant list to obtain the original list back.  The original list is referred to as the time-domain signal, and the list as returned by Fourier is referred to as the frequency-domain signal.

To implement our filters, we'll simply find the frequency-domain representation of a signal, mask off the unwanted frequencies, and convert the result back to the time domain.

The mask will consist of a list containing 1 in positions corresponding to frequencies which are to be retained and 0 for those to be rejected, and will simply be multiplied with the frequency-domain signal to remove unwanted frequencies.

Preliminary Setup

The following functions are used to optimize the Fourier Transforms that we use in this Notebook.  Because Fourier and InverseFourier perform most efficiently when they operate on arrays whose sizes are products of perfect powers of small prime numbers, we will pad our samples out to the nearest such size, using the functions defined here.

$foo will contain a list of optimal sizes, up to a maximum of $maxsz.  Sizes larger than $maxsz will be generated using powers of two times sizes which appear in $foo.  The function

	newsize = psp[origsize]

is used to find optimal sizes.

$maxsz = 22050;

$foo = Module[
		{c2 = Log[2.], c3 = Log[3.], c5 = Log[5.],
			c1000 = Log[$maxsz + 0.5]},
		Sort[Flatten[Round[Exp[Table[i c2 + j c3 + k c5,
			{i, 0, Floor[c1000/c2]},
			{j, 0, Floor[(c1000 - i c2)/c3]},
			{k, 0, Floor[(c1000 - i c2 - j c3)/c5]}]]]]]];

psp[x_] := Module[{p, q},
			If[x < $maxsz,
				First[Select[$foo, (# >= x)&]],
			(* else *)
				p = Ceiling[Log[2, N[x]/$maxsz]];
				q = 2^-p;
				2^p First[Select[$foo, (# >= q x)&]]
			]
		]

Band Pass Filter

First, we define a simple band pass filter.  It takes a SampledSoundList and two corner frequencies as arguments, and returns a SampledSoundList as its result.

The two corner frequencies define a "pass band", or a group of frequencies which will be retained from the original sound in the result.  Frequencies outside the pass band (lower than the lower corner and higher than the upper corner) will be rejected.

BandPassFilter[sound_SampledSoundList, lcorner_, hcorner_] :=
	Block[{samps, srate, origsize, newsize, kern,
		 	nlcorn, nhcorn, n},
		
		(* Separate the sound into samples and sample rate *)
		samps = sound[[1]];
		srate = sound[[2]];
		
		(* Because the FFT algorithm is much more efficient
		   when it operates on tables whose sizes are products
		   of perfect powers of small prime numbers, we will
		   pad the original sample with silence to be of such
		   a size.  After we perform the filtering, we will
		   return the result table to the correct size. *)
		   
		origsize = Length[samps];
		(* Calculate optimum size for FFT and pad sound *)
		newsize = psp[origsize];
		samps = Join[samps, Table[0.0, {newsize - origsize}]];
		
		(* The Fourier[] function generates a table of 
		   frequencies relative to the size of the original
		   table.  Therefore, frequencies in the frequency
		   domain must be scaled to relate to a "real" time-
		   domain signal. *)
		   
		nlcorn = (lcorner * newsize) / srate;
		nhcorn = (hcorner * newsize) / srate;
		
		(* Generate the frequency response table.  This
		   is a simple table in the frequency-domain table
		   which will be multiplied with the frequency spectrum
		   of the original sound.  Thus, it contains 1 for
		   frequencies which are to be kept and 0 for those
		   that are to be rejected.  *)
		
		kern = Table[If[n >= nlcorn && n <= nhcorn, 1.0, 0.0],
				{n, Ceiling[newsize/2.0]}];
		(* The table is symmetric due to aliasing. *)
		kern = Join[kern,
					Take[Reverse[kern], Floor[newsize/2.0]]];
		
		(* Generate the new samples by conversion to the
		   frequency domain, performing multiplication and
		   conversion back to the time domain. *)
		   
		samps = Take[Chop[Re[
					InverseFourier[kern Fourier[samps]]]],
					origsize];
		
		(* Make sure all the samples are still within range,
		   because some destructive interference may have
		   been lost in the filter. *)
		
		samps /= Max[Abs[samps]];
		
		(* Build the SampledSoundList back from the samples *)
		Return[SampledSoundList[samps, srate]]
	]

Demonstration of the Band Pass Filter

In order to demonstrate the band pass filter, we'll construct a sound that is composed of three frequencies, at 400 Hz, 800 Hz and 1600 Hz.  Just to keep things small, it will be only 1/20 second long.

sines = SampledSoundList[
		Table[N[0.3 (Sin[2 Pi 400. x] + Sin[2 Pi 800. x] + Sin[2 Pi 1600. x])],
		{x, 0, 0.05, 1/22050.0}], 22050];

Show[Sound[sines]]

Now, we'll filter to retain only the lowest of the frequencies using our band pass filter we constructed above.  In this case, we're actually using it as a low pass filter, by specifiying the lower edge of the band to be at 0 Hz.

We'll keep the 400 Hz tone by passing frequencies between 0 and 600 Hz.

lsine = BandPassFilter[sines, 0, 600];

Show[Sound[lsine]]

Next, we'll extract the 800 Hz tone, using a passband of 600-1200 Hz.

msine = BandPassFilter[sines, 600, 1200];

Show[Sound[msine]]

And finally, we extract the high tone, using a passband of 1200-2400 Hz.

hsine = BandPassFilter[sines, 1200, 2400];

Show[Sound[hsine]]

Just to make sure we've done everything correctly, we add the three extracted tones together.  We should get the same sound we started with.

ssine = SampledSoundList[lsine[[1]] + msine[[1]] + hsine[[1]], 22050];

Show[Sound[ssine]]

Low Pass Filter

The low pass filter is simply a band pass filter that sets the lower edge of the pass band at 0 Hz.

LowPassFilter[sound_SampledSoundList, cutoff_] :=
		BandPassFilter[sound, 0, cutoff]

Demonstration of the Low Pass Filter

To demonstrate the low pass filter, we'll use a sound that sweeps from high frequencies to low frequencies.  Here is the unfiltered sound:

sweep = SampledSoundList[N[Table[Exp[-50 t] Sin[10/(t + 0.01)],
	{t, 0.0, 0.05, 1./22050}]], 22050];

Show[Sound[sweep]]

Now, we will filter out all frequencies above 1800 Hz.

lpsweep = LowPassFilter[sweep, 1800];

Show[Sound[lpsweep]]

High Pass Filter

The high pass filter is simply a band pass filter that sets the upper edge of the pass band at Infinity.

HighPassFilter[sound_SampledSoundList, cutoff_] :=
		BandPassFilter[sound, cutoff, Infinity]

Demonstration of the High Pass Filter

To demonstrate the high pass filter, we'll use the same sound we used for the low pass filter, one that sweeps from high to low frequencies.

hpsweep = HighPassFilter[sweep, 1800];

Show[Sound[hpsweep]]

Notch Reject Filter

The Notch Reject filter is almost exactly the opposite of the band pass filter: instead of retaining all frequencies between two "corner" frequencies and rejecting all others, it rejects frequencies in a given band and retains all others.  Its implementation is very nearly identical to the band pass filter, but with one important change: the sense of the "If" statement in building the filter.

NotchRejectFilter[sound_SampledSoundList, lcorner_, hcorner_] :=
	Block[
		{samps, srate, origsize, newsize, kern,
		 nlcorn, nhcorn, n},
		samps = sound[[1]];
		srate = sound[[2]];
		origsize = Length[samps];
		newsize = psp[origsize];
		samps = Join[samps, Table[0.0, {newsize - origsize}]];
		nlcorn = (lcorner * newsize) / srate;
		nhcorn = (hcorner * newsize) / srate;
		kern = Table[If[n >= nlcorn && n <= nhcorn, 0.0, 1.0],
				{n, Ceiling[newsize/2.0]}];
		kern = Join[kern, Take[Reverse[kern],
				Floor[newsize/2.0]]];
		samps = Take[
					Chop[Re[
						InverseFourier[kern Fourier[samps]]]],
					origsize];
		samps /= Max[Abs[samps]];
		Return[SampledSoundList[samps, srate]]
	]

Demonstration of the Notch Reject Filter

We'll use the same sound to demonstrate the notch reject filter that we used to demonstrate the low pass and high pass filters.  This time we'll remove everything between 800 and 2500 Hz.

nrsweep = NotchRejectFilter[sweep, 800, 2500];

Show[Sound[nrsweep]]