Beginner Attempting FFT Signal Filtering With Numpy

Question

I've tried looking around for information on this, but I'm really out of my league here. I'm a guy who likes to fool around with Python, and I wanted to make a program that would filter an audio file. I'm using Python and NumPy, with the scipy.io.wavfile module for importing and exporting Wave files. I've gotten it messing around with volume, but not filtering. Here's what I have so far.

rate, data = wavfile.read('./TriLeftChannel.wav') 

filtereddata = numpy.fft.rfft(data) # FFT Filtered data
freqdata = numpy.fft.fftfreq(data.size) #Frequency data
filtereddata = AudioFunctions.Filter(filtereddata, freqdata, data, rate) # Filter the data

def Filter(filtereddata, freqdata, data, rate):

    #fftchunks = (rate / 2)# + 1

    x = freqdata[:len(data) / 2]

    for f in range(len(x)):

        if x[f] > 0.1:
            filtereddata[f] = 0.0
            filtereddata[len(data) / 2 + f] = 0.0       

    return filtereddata

In that function, filtereddata is the FFT'd data, freqdata is the frequency data that I got with fftfreq(), and data is the wave file itself, 'bare'. Rate is the sampling rate (though I don't use it). This function doesn't actually filter the frequencies (although I know it's a hard filter and no filter should really be this harsh). After I'm done, I output the file with

filteredwrite = numpy.fft.irfft(filtereddata)

filteredwrite = numpy.round(filteredwrite).astype('int16') # Round off the numbers, and get ready to save it as 16-bit depth file (?)

wavfile.write('TestFiltered.wav', rate, filteredwrite)

Since I'm kind of struggling even just to have gotten this far, I was wondering if anyone could give me any pointers or a kind of beginners' tutorial for FFT? Of course, any help would be greatly appreciated.

EDIT: Fixed indentation issues.

Answer 1

Two potential problems with your approach:

• You are computing the FFT on your whole signal, which will be terribly inefficient if your input data gets too large, and prevents you from using windowing. If you want to do frequency-domain modifications of your signal, consider using a Short-Term Fourier Transform, processing the resulting FFT frames, and resynthesizing a time-domain signal using overlap-add.
• FFT -> zeroing FFT coefficients -> IFFT, especially without windowing of the input data, is seldom used for filtering as it will yield a filter with many unwanted characteristics (side-lobes + non-causal).

A potential problem in your code:

I don't understand why you zero the filtereddata[len(data) / 2 + f] coefficients. I assume it's intended to zero the "negative frequencies" but then the array index is wrong, and you are using numpy.fft.rfft which has a different output layout from numpy.fft.fft (FFT coefficients for negative frequencies are not returned). I'm surprised this runs at all, your code should actually cause an index out of bounds error.

My suggestion would be for you to learn about FIR / IIR filters, design a suitable filter matching your constraints, and apply it to your data. There is a well-established body of work, software tools and practices on FIR / IIR filters and their design. Once you have the coefficients, for small sequences and/or if complexity is not an issue, there's scipy.signal.lfilter to do the job. As suggested by Jason R, applying a long FIR can be made very efficient computationally with FFTs + overlap-add.

Answer 2

I have had good time implementing what this guy says:

http://mpastell.com/2010/01/18/fir-with-scipy/

A typical low pass filter I use, got from the given link, is this:

def firfilt(interval, freq, sampling_rate):
    nfreq = freq/(0.5*sampling_rate)
    taps =  sampling_rate + 1
    a = 1
    b = scipy.signal.firwin(taps, cutoff=nfreq)
    return scipy.signal.lfilter(b, a, interval)

Answer 3

Linear FIR filters are applied to a signal (like your audio file) using discrete convolution. Convolution can be implemented efficiently using the FFT. Two separate schemes for doing this are called the overlap-save and overlap-add methods. I personally prefer overlap-save, as it's a bit simpler to implement. It's not clear from your question exactly what you're getting hung up on. You should have all the high-level tools with NumPy that allow you to readily put together the OLS or OLA algorithms.

While the techniques are pretty simple, it pays to try to get a basic understanding of what's happening under the hood and why they work. For a good introductory DSP book, I would recommend Lyons' DSP book. Alternatively, there is an online book called "The Scientist and Engineer's Guide to Digital Signal Processing" that I've heard good things about (you can't beat the price).