DiscreteTransforms.mw
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
:
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:
>

ArrayTools:Copy(
data,
CdataR,
0, 2 ); 
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 inplace
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
); 
>

limv
:= floor( ( nSamples1
)/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 Corder 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
:= nSamplesLowHalfSamp;

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*( nSamples1
)4*hh,
4, Cdata2,
2*( nSamples1
)2*hh,
2 ); 
> 
ArrayTools:Copy(
hh,
Cdata,
2*( nSamples1
)4*hh+1,
4, Cdata2,
2*( nSamples1
)2*hh+1,
2 ); 
To check,
we plot a small range of one channel:
> 
tscale
:= evalf( nSamples/Freq
); 
> 
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 ); 