20

So i recently successfully built a system which will record, plot, and playback an audio wav file entirely with python. Now, I'm trying to put some filtering and audio mixing in between the when i record and when i start plotting and outputting the file to the speakers. However, i have no idea where to start. Right now I'm to read in a the intial wav file, apply a low pass filter, and then re-pack the newly filtered data into a new wav file. Here is the code i used to plot the initial data once i recorded it.

import matplotlib.pyplot as plt
import numpy as np
import wave
import sys

spf = wave.open('wavfile.wav','r')

#Extract Raw Audio from Wav File
signal = spf.readframes(-1)
signal = np.fromstring(signal, 'Int16')

plt.figure(1)
plt.title('Signal Wave...')
plt.plot(signal)

And here is some code i used to generate a test audio file of a single tone:

import numpy as np
import wave
import struct

freq = 440.0
data_size = 40000
fname = "High_A.wav"
frate = 11025.0  
amp = 64000.0    

sine_list_x = []
for x in range(data_size):
    sine_list_x.append(np.sin(2*np.pi*freq*(x/frate)))

wav_file = wave.open(fname, "w")

nchannels = 1
sampwidth = 2
framerate = int(frate)
nframes = data_size
comptype = "NONE"
compname = "not compressed"

wav_file.setparams((nchannels, sampwidth, framerate, nframes,
comptype, compname))

for s in sine_list_x:
    wav_file.writeframes(struct.pack('h', int(s*amp/2)))

wav_file.close()

I'm not really sure how to apply said audio filter and repack it, though. Any help and/or advice you could offer would be greatly appreciated.

Philip R.
  • 345
  • 3
  • 4
  • 12
  • 2
    Have you tried looking at `scipy`'s [lfilter](http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html)? – dpwilson Jul 23 '14 at 20:35
  • 2
    Rather than the for loop to generate the sinusoid, you want something like `sine_signal = np.sin(2*np.pi*freq*(np.arange(data_size)/frate))`, then something like `wav_file.writeframes((sine_signal*amp/2).astype('h').tostring())`. – dpwe Jul 24 '14 at 14:27

2 Answers2

57

First step : What kind of audio filter do you need ?

Choose the filtered band

For the following steps, i assume you need a Low-pass Filter.

Choose your cutoff frequency

The Cutoff frequency is the frequency where your signal will be attenuated by -3dB.

Your example signal is 440Hz, so let's choose a Cutoff frequency of 400Hz. Then your 440Hz-signal is attenuated (more than -3dB), by the Low-pass 400Hz filter.

Choose your filter type

According to this other stackoverflow answer

Filter design is beyond the scope of Stack Overflow - that's a DSP problem, not a programming problem. Filter design is covered by any DSP textbook - go to your library. I like Proakis and Manolakis' Digital Signal Processing. (Ifeachor and Jervis' Digital Signal Processing isn't bad either.)

To go inside a simple example, I suggest to use a moving average filter (for a simple low-pass filter).

See Moving average

Mathematically, a moving average is a type of convolution and so it can be viewed as an example of a low-pass filter used in signal processing

This Moving average Low-pass Filter is a basic filter, and it is quite easy to use and to understand.

The parameter of the moving average is the window length.

The relationship between moving average window length and Cutoff frequency needs little bit mathematics and is explained here

The code will be

import math

sampleRate = 11025.0 
cutOffFrequency = 400.0
freqRatio = cutOffFrequency / sampleRate

N = int(math.sqrt(0.196201 + freqRatio**2) / freqRatio)

So, in the example, the window length will be 12

Second step : coding the filter

Hand-made moving average

see specific discussion on how to create a moving average in python

Solution from Alleo is

def running_mean(x, windowSize):
   cumsum = numpy.cumsum(numpy.insert(x, 0, 0)) 
   return (cumsum[windowSize:] - cumsum[:-windowSize]) / windowSize 

filtered = running_mean(signal, N)

Using lfilter

Alternatively, as suggested by dpwilson, we can also use lfilter

win = numpy.ones(N)
win *= 1.0/N
filtered = scipy.signal.lfilter(win, [1], signal).astype(channels.dtype)

Third step : Let's Put It All Together

import matplotlib.pyplot as plt
import numpy as np
import wave
import sys
import math
import contextlib

fname = 'test.wav'
outname = 'filtered.wav'

cutOffFrequency = 400.0

# from http://stackoverflow.com/questions/13728392/moving-average-or-running-mean
def running_mean(x, windowSize):
  cumsum = np.cumsum(np.insert(x, 0, 0)) 
  return (cumsum[windowSize:] - cumsum[:-windowSize]) / windowSize

# from http://stackoverflow.com/questions/2226853/interpreting-wav-data/2227174#2227174
def interpret_wav(raw_bytes, n_frames, n_channels, sample_width, interleaved = True):

    if sample_width == 1:
        dtype = np.uint8 # unsigned char
    elif sample_width == 2:
        dtype = np.int16 # signed 2-byte short
    else:
        raise ValueError("Only supports 8 and 16 bit audio formats.")

    channels = np.fromstring(raw_bytes, dtype=dtype)

    if interleaved:
        # channels are interleaved, i.e. sample N of channel M follows sample N of channel M-1 in raw data
        channels.shape = (n_frames, n_channels)
        channels = channels.T
    else:
        # channels are not interleaved. All samples from channel M occur before all samples from channel M-1
        channels.shape = (n_channels, n_frames)

    return channels

with contextlib.closing(wave.open(fname,'rb')) as spf:
    sampleRate = spf.getframerate()
    ampWidth = spf.getsampwidth()
    nChannels = spf.getnchannels()
    nFrames = spf.getnframes()

    # Extract Raw Audio from multi-channel Wav File
    signal = spf.readframes(nFrames*nChannels)
    spf.close()
    channels = interpret_wav(signal, nFrames, nChannels, ampWidth, True)

    # get window size
    # from http://dsp.stackexchange.com/questions/9966/what-is-the-cut-off-frequency-of-a-moving-average-filter
    freqRatio = (cutOffFrequency/sampleRate)
    N = int(math.sqrt(0.196196 + freqRatio**2)/freqRatio)

    # Use moviung average (only on first channel)
    filtered = running_mean(channels[0], N).astype(channels.dtype)

    wav_file = wave.open(outname, "w")
    wav_file.setparams((1, ampWidth, sampleRate, nFrames, spf.getcomptype(), spf.getcompname()))
    wav_file.writeframes(filtered.tobytes('C'))
    wav_file.close()
ankostis
  • 8,579
  • 3
  • 47
  • 61
piercus
  • 1,136
  • 9
  • 17
  • Hi, I think your code has accidentally become a high-pass filter. I used it on a wav audio and ended up removing the bass. I need both low and high pass filters, so could you point me to the part that needs to be changed? – naman1901 May 28 '17 at 17:39
  • I confirm it is working as a low pass on my end. Can you share the wav file you are using ? – piercus Jun 02 '17 at 17:04
  • Sorry, my bad. It is working as an LPF. Apparently the bass was a higher frequency than the foreground. – naman1901 Jun 27 '17 at 06:28
  • What is the meaning of `[1]` in `win, [1], signal`? – user13107 Feb 05 '18 at 10:05
  • @user13107 Doesn't GitHub's terms conflict with Stack Overflow's license? – BSMP Feb 05 '18 at 16:34
  • @BSMP I tried to search for Meta question where this issue (`copying code from github to SO`) is handled; but I couldnt find it. Note that the github code is also written by OP. Can we just mention that this code is taken from GitHub and same license is applicable here? – user13107 Feb 06 '18 at 02:28
  • 2
    @user13107 see https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html#scipy-signal-lfilter in `nfilter`, we use a more generic concept of IIR or FIR filter. This is a polynomial defined-filter per the equation `a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M] - a[1]*y[n-1] - ... - a[N]*y[n-N]` We use `b = win, a = [1], x = signal` then equation is `1*y[n] = win[0]*signal[n] + win[1]*signal[n-1] + ... + win[M]*signal[n-M]` It's like using a sledgehammer to crack a nut but it works ! – piercus Feb 06 '18 at 11:41
  • I use this to filter certain instruments from wav-files and it works like a charm. Only thing: it gives me a depraction warning: DeprecationWarning: The binary mode of fromstring is deprecated, as it behaves surprisingly on unicode inputs. Use frombuffer instead channels = np.fromstring(raw_bytes, dtype=dtype). Anyone knows how to resolve this? – Philip F. Dec 24 '20 at 11:02
2

sox library can be used for static noise removal.

I found this gist which has some useful commands as examples

penduDev
  • 4,743
  • 35
  • 37