[SciPy-User] butterworth filter on .WAV file

[Peter Howard]

I'm trying to write a very simple example of applying a band pass filter to a .WAV music file.
I'm distinctly rusty on DSP and inexperienced with SciPy/NumPy so apologies if I've made a dumb mistake.

It executes without errors or warnings.
It produces the output file, but this is twice the size of the input file, which is clearly wrong.
I'm most uncertain about casting the filtered data back to integers and thus being suitable for writing back out to .WAV.
I'm also bit uncertain about my interpretation / understanding of the frequency and gain specifications.

Any help and advice very much appreciated.

Pete


<snip>

from scipy.io.wavfile import read, write
from scipy.signal.filter_design import butter, buttord
from scipy.signal import lfilter
from numpy import asarray

def convert_hertz(freq):
    # convert frequency in hz to units of pi rad/sample
    # (our .WAV is sampled at 44.1KHz)
    return freq * 2.0 / 44100.0

rate, sound_samples = read('monty.wav')
pass_freq = convert_hertz(440.0) # pass up to 'middle C'
stop_freq = convert_hertz(440.0 * 4) # max attenuation from 3 octaves higher
pass_gain = 3.0 # tolerable loss in passband (dB)
stop_gain = 60.0 # required attenuation in stopband (dB)
ord, wn = buttord(pass_freq, stop_freq, pass_gain, stop_gain)
b, a = butter(ord, wn, btype = 'low')
filtered = lfilter(b, a, sound_samples)
integerised_filtered = asarray(filtered, int)
write('monty-filtered.wav', rate, integerised_filtered)

[Joshua Holbrook]

> _______________________________________________
> SciPy-User mailing list
> SciPy...@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>

What does the output file sound like?

--Josh
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[Peter Howard]

I forgot to mention that part - it plays without complaint in Window's Media Player - but is completely silent.
Would it help if I provide the input .WAV file?
Can you attach things to posts via email?
(It is 13MB)

[Joshua Holbrook]

> _______________________________________________
> SciPy-User mailing list
> SciPy...@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>

That's huge! Nevermind that.

Interestingly enough, I tried messing with the problem using a random
wav file I found from the internet (that plays in vlc just fine) and
got an ugly error. So, I guess I can't help you in this way :S

That said, here's what I would do: First, plot the array that results
from reading the wav file to get an idea of what it looks like. Then,
after you butterworth it up, plot the result of THAT to see what you
have. That could give you an idea as to whether the problem is in the
reading, the butterworthing, the output, or whatever.

Just my $0.02. Good luck!

--Josh

p.s:

If anyone cares about the error I got, here you go:

In [12]: a,b=wavfile.read('viacom2.wav')
Reading fmt chunk
Reading data chunk
Warning: %s chunk not understood
---------------------------------------------------------------------------
error Traceback (most recent call last)

/home/josh/<ipython console> in <module>()

/usr/lib/python2.6/site-packages/scipy/io/wavfile.pyc in read(file)
65 else:
66 print "Warning: %s chunk not understood"
---> 67 size = struct.unpack('I',fid.read(4))[0]
68 bytes = fid.read(size)
69 fid.close()

error: unpack requires a string argument of length 4

The file could be borked for all I know, and maybe this was fixed a
long time ago and my scipy package (probably stock fedora) could be a
bit stale, but the %s bit makes me think that maybe that part wasn't
completely finished? Anyways.

[Benjamin Root]

Well, to give you hope, it is feasible as some friends of mine did a notch filter on some .wav files at the conference last week.  They had a bunch of issues as well.  I will see if he is actively reading this list.

Ben Root

[David Baddeley]

In addition to the comments that have already been made (suggesting you plot the result etc ...), I think you might have problems with your final cast (the asarray(filtered, int) bit). On my computer that results in a 64 bit integer type (might be 32 bit on some platforms). I guess you want a 16 bit integer, as this is what most sound programs will be expecting. You could try:

asarray(filtered, 'int16)

or 

filtered.astype('int16')

cheers,

David

---

From: Peter Howard <peterh...@gmail.com>
To: scipy...@scipy.org
Sent: Thu, 8 July, 2010 6:40:31 AM
Subject: [SciPy-User] butterworth filter on .WAV file

[Christopher Brown]

I've had good luck with:

# Read
fs,data = read(infilename)
data = np.float64(data/32768.)

# ... process ...

# Write
write(outfilename, fs, np.int16(data*32768))

--
Christopher Brown, Ph.D.
Associate Research Professor
Department of Speech and Hearing Science
Arizona State University
http://pal.asu.edu

[Peter Howard]

This tip from Chris, has nearly solved the problem.
Thanks to everybody that has contributed.
I now get sound out that sounds similar to but different from the input - so the fundamentals are sorted thank you.
However I can only get it to work by experimentation with the frequencies and gains - it doesn't work with my understanding of what they theoretically should be. Either the filter design module objects to the coefficients or I get silence.
Also, (and almost certainly significantly) the filtered output is only on the right stereo channel.
I wonder if each integer sound sample in a WAV file is split bitwise into left and right sectors and the type conversion is corrupting half of each sample?

Here's the state of play:

<snip>
from scipy.io.wavfile import read, write
from scipy.signal.filter_design import butter, buttord

from scipy.signal import lfilter, lfiltic
import numpy as np
from math import log

rate, sound_samples = read('monty.wav')

sound_samples = np.float64(sound_samples / 32768.0)
pass_freq = 0.2
stop_freq = 0.3
pass_gain = 0.5 # permissible loss (ripple) in passband (dB)
stop_gain = 10.0 # attenuation required in stopband (dB)

ord, wn = buttord(pass_freq, stop_freq, pass_gain, stop_gain)
b, a = butter(ord, wn, btype = 'low')
filtered = lfilter(b, a, sound_samples)

filtered = np.int16(filtered * 32768 * 10)
write('monty-filtered.wav', rate, filtered)





Pete

[Sebastian Haase]

You could "post" (upload) large files maybe to http://drop.io for
everyone interested to listen in (100MB file limit)
-Sebastian

[Peter Howard]

Hang on a minute - my brain is working backwards... am I wasting everyone's time here...

If a .WAV file contains stereo sound data - then it's never going to work sucking it in and simply applying a filter to the numbers is it?

So how come it sort of works?

Or is the entire NumPy/SciPy suite so clever about working in N-dimensions that it treats the two sound channels as multi-dimensioned arrays right from the .WAV read() call? - and the filtering is separate for each dimension?

Pete

[Fabrice Silva]

So you might consider NumPy/SciPy clever! When reading a wav file, the
output array is a 2D array with shape (M,N) where M is the number of
samples of each channel (time range*sampling frequency) and N is the
number of channels. Each channel is stored in a column of the output
array.

And it is so clever that it handles to apply a filter (with
scipy.signal.lfilter) on each of the columns of the array.

Example
In [2]: import scipy.io.wavfile as wv

In [3]: Fs,Sig = wv.read("STE-023.wav")

Warning: %s chunk not understood
Reading fmt chunk
Reading data chunk

In [4]: Fs
Out[4]: 44100

In [5]: Sig.shape, Sig.dtype
Out[5]: ((4434112, 2), dtype('int16'))

In [6]: import scipy.signal as ss

In [8]: SigFilt = ss.lfilter([1],[1, .1], Sig)
Out[8]:
array([[ 0. , 0. ],
[ 4. , -9.4],
[ 7. , -20.7],
...,
[ 0. , 0. ],
[ 0. , 0. ],
[ 0. , 0. ]])

In [9]: SigFilt.shape, SigFilt.dtype
Out[9]: ((4434112, 2), dtype('float64'))