• The WAV package
• Analyzing the Windows "Ding" sound file
• Altering the sound file

Discrete Transforms in Maple 9

by Maplesoft

The new DiscreteTransforms package in Maple 9 provides routines for Fast Fourier Transforms that are many times faster and more flexible than the older FFT  routines of previous releases.  In this demonstration, we use the DiscreteTransforms and ArrayTools packages to analyze a .wav sound file ( the Windows "Ding" ) , alter the signal, and write out a new sound file.

For manipulating the .wav files, we use an enhanced version of the WAV package from the Application Center.

The WAV package

Analyzing the Windows "Ding" sound file

As a sound file, we select "ding.wav", which is the sound that Windows produces when it wants to call your attention to something.

First we read in the data from the WAV file:

>	data := WAV:-ReadWAV( "C:/ding.wav", 'Freq' );

We see from the data storage that this sound sample is in stereo and has a total of 20191 samples. Thus, the duration of the sound in seconds is:  

>	nSamples := op( [1, 2], [rtable_dims( data )] ):

>	evalf( nSamples/Freq );

So a little less than 1 second.

We can plot a small sample of the sound data as follows:

>	plot( [seq( [i, data[i, 1]], i=5000..5500 )] );

Let's look at the data in the frequency domain using a discrete fourier transform.

We also want to look at the two different stereo channels at the same time. We can do this efficiently as follows:

Construct a complex Array for the transform data ( used for input and output ).

>	Cdata := Array( 1..nSamples, 1..2, datatype=complex[8], order=C_order );

Now create an alias to real data for copying the WAV data:

>	CdataR := ArrayTools:-ComplexAsFloat( Cdata );

The following step is not necessary, but we do it for clarity.  Create an alias of the real data to an Array of the same form as the WAV data for the copy:

>	CdataRA := ArrayTools:-Alias( CdataR, [1..nSamples, 1..2], C_order );

Copy the data over

>	ArrayTools:-Copy( data, CdataRA, 0, 2 );

Check :

>	data[5000, 1..2];

>	Cdata[5000, 1..2];

Exactly what is wanted.

Note: as mentioned, the Alias step is not needed prior to the copy - we could have used CdataR instead of CdataRA. To demonstrate this:

>	Cdata[5000, 1]:=0:

>	ArrayTools:-Copy( data, CdataR, 0, 2 );

>	Cdata[5000, 1];

Exactly as before. This is because the copy routine acts directly on the raw storage of the Vector, Matrix or Array it is working with. This has consequences for the data ordering of any Array or Matrix data being used. For more information, see C_order, Fortran_order .

Now the data is in the desired form, so we transform it using an in-place fast discrete fourier transform. We specify that the dimension is the first ( 1..20191 ):

>	DiscreteTransforms:-FourierTransform( Cdata, 1, inplace=true );

We can now examine the frequency domain data in a plot.  

Since we want the plot frequency in Hz, we need to adjust based on the number of samples and the sample frequency, that is, the time duration of the signal ( the ratio of the two ). In addition, array element 1 is actually the frequency 0 component, so we shift by 1.

Finally, there is far too much data to effectively plot, so we will only plot every 10th value, and only over the first half range ( as this is all that is needed when the data is real ).

>	tscale := evalf( nSamples/Freq );

>	skip := 10:

>	limv := floor( ( nSamples-1 )/skip/2 ):

>	freq1 := [seq( [i*skip/tscale, abs( Cdata[i*skip+1, 1] )], i=1..limv )]:

>	freq2 := [seq( [i*skip/tscale, abs( Cdata[i*skip+1, 2] )], i=1..limv )]:

>	plots[display]( {plot( freq1 ), plot( freq2 )}, labels=["Frequency ( Hz )", "A"] );

We see a dominant frequency somewhere around 800Hz, a smaller one around 3200 Hz, and an even smaller one around 6900 Hz.  We can extract these from our data:

>	select( a->a[2]>0.01, freq1 );

>

The first is around 786.3 Hz, the second around 3170 Hz and the third around 6935 Hz.

On the musical scale, the first spike is roughly at a G ( 783.99 Hz ), and the second is approximately two octaves higher at around 3135.96 Hz.

Altering the sound file

Let's drop the 'ding' by one octave, keeping the duration the same. We do this by halving the frequency of the data, or in transform terms, map freq 2*i to i. This also means that the top half of the frequency will be zero, as we have no data for it. We can do this with ArrayTools:-Copy.

>	Cdata2 := Array( 1..nSamples, 1..2, datatype=complex[8], order=C_order, fill=0 );

We are dealing with a nSamples x 2 C-order array, so the left and right channels of the data are interleaved.
In the first step we will copy the relevant portion of the first half of the left data to the new data array.

The data is symmetric, so we need to copy the first half to the start of the Array, and the last half to the end, to get the correct result:

>	LowHalfSamp := iquo( nSamples+1, 2 );

>	HighHalfSamp := nSamples-LowHalfSamp;

Since we are halving the frequency, we only want to copy every second sample.

Taking into account the interleaving of the left and right channels, we do this for the low frequency portion of the left channel as follows:

>	ArrayTools:-Copy( iquo( LowHalfSamp, 2 ), Cdata, 0, 4, Cdata2, 0, 2 );

This command copies the first iquo( LowHalfSamp )entries from Cdata starting from 0, skips by 4 at each step to Cdata2 starting from 0, and skips by 2 in each step.

It is equivalent to the loop:

for i from 0 to iquo( LowHalfSamp, 2 )-1 do

 Cdata2[0+i+1, 1] := Cdata2[0+2*i+1, 1]

end do

The right channel is similar:

>	ArrayTools:-Copy( iquo( LowHalfSamp, 2 ), Cdata, 1, 4, Cdata2, 1, 2 );

The upper frequencies are a little trickier ( with respect to the 'offset' ), but can be done as follows:

>	hh := iquo( HighHalfSamp, 2 ):

>	ArrayTools:-Copy( hh, Cdata, 2*( nSamples-1 )-4*hh, 4, Cdata2, 2*( nSamples-1 )-2*hh, 2 );

>	ArrayTools:-Copy( hh, Cdata, 2*( nSamples-1 )-4*hh+1, 4, Cdata2, 2*( nSamples-1 )-2*hh+1, 2 );

To check, we plot a small range of one channel:

>	tscale := evalf( nSamples/Freq );

>	skip := 10:

>	limv := 50:

>	plot( [seq( [i*skip/tscale, abs( Cdata2[i*skip+1, 1] )], i=1..limv )] );

We now have a spike at around 400 Hz, which is what we wanted.
Now we apply the inverse transform to map to the time domain:

>	DiscreteTransforms:-InverseFourierTransform( Cdata2, 1, inplace=true );

Use ArrayTools to ComplexAsFloat and copy the data to a real array, which we will use to construct our new WAV file:

>	Cdata2R := ArrayTools:-ComplexAsFloat( Cdata2 );

>	data2 := hfarray( 1..nSamples, 1..2 ):

>	ArrayTools:-Copy( nSamples, Cdata2R, 0, 4, data2, 0, 2 );

>	ArrayTools:-Copy( nSamples, Cdata2R, 2, 4, data2, 1, 2 );

And now construct the WAV file with the WAV package( and no normalization ):

>	WAV:-WriteWAV( "C:/ding_low.wav", data2, Freq, 8, false );

>