2

I have been working on audio digital signal processing. I wish to design a digital filter. I have attached the screenshot. Of course, I can design band pass or bandstop filter using FIR/IIR filters, but what is special about this filter is that it controls low, mid and high frequencies. So my question, how can I design such filter in python enter image description here manually with this magnitude response as shown. Please help. I dont know how do I limit the filter to those values.

Another alternartive could be to such filter in MATLAB's filterdesigner tool and export the coefficients to python. But again, I dont know how to design such filter there. So please help.

Scott Stensland
  • 26,870
  • 12
  • 93
  • 104
Sanket Jain
  • 21
  • 2

1 Answers1

1

There are a number of design options to match more arbitrary frequency band specifications than the simple lowpass/bandpass/bandstop/hipass classification. The most common ones are :

The first two are readily available from scipy.signal as scipy.firls and scipy.remez respectively, whereas the third can be done by taking the inverse FFT of the sampled frequency response curve.

Sample implementation using the least-square method:

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

bands = [0, fL1, fL2, fH1, fH2, fs/2]
desired = [Bl, Bl, 1, 1, Bh, Bh]
y = signal.firls(numtaps, bands, desired, fs=fs)
f,h = signal.freqz(y, fs=fs)

plt.plot(bands, desired, 'k', linewidth=2)
plt.plot(f, np.abs(h), 'r')

Sample implementation using the Parks-McClellan algorithm:

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

bands = [0, fL1, fL2, fH1, fH2, fs/2]
desired = [Bl, 1, Bh]
y = signal.remez(numtaps, bands, desired, fs=fs)
f,h = signal.freqz(y, fs=fs)

plt.plot(bands, desired, 'k', linewidth=2)
plt.plot(f, np.abs(h), 'b')

Sample implementation using the frequency sampling method:

from scipy import interpolate, fft
import matplotlib.pyplot as plt
import numpy as np

# interpolate desired response at equally spaced frequency points
interpolator = interpolate.interp1d(bands, desired)
N = 1024
fsampling = np.linspace(0, fs/2, N)
sampled = interpolator(fsampling)
# take the inverse FFT
y = fft.fftshift(fft.irfft(sampled))
# truncate response to keep numtaps coefficients
n1 = N-numtaps//2
n2 = n1 + numtaps
y = y[n1:n2]

For sake of illustration, using the following set of parameters

numtaps = 19
fL1 = 100
fL2 = 200
fH1 = 500
fH2 = 600
fs  = 1600
Bl  = 2
Bh  = 2.3

would yield the corresponding designs:

Frequency responses using least-squares, Parks-McClellan and frequency sampling design methods

Note that if your frequency specification curves are do-not-exceed values (rather than nominal ones), you can still use the above methods in an iterative fashion each time tweaking the bands and desired parameters until you converge to an acceptable response. With the initial set of parameters above, using the least-squares method, you could end up with tweaked parameters such as

bands   = [0, 100, 180, 520, 600, fs/2]
desired = [1.966,1.966,0.980,0.980,2.225,2.225]

with the corresponding response:

enter image description here

SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • Thank you for your detailed answer. I wish to implement this filter in an equation. For example, I have equation( B'(k)*B(k) ), where ' is conjugate for example and B(k) is above such filter calculated for k in range(0,4096) ( 4096 point FFT suppose). Can you please tell me then how would I implement this in python with above filter which you just answered? – Sanket Jain Dec 06 '21 at 10:58
  • Specifically for that equation you could try using [`numpy.conjugate`](https://numpy.org/doc/stable/reference/generated/numpy.conjugate.html) and [`numpy.multiply`](https://numpy.org/doc/stable/reference/generated/numpy.multiply.html)... – SleuthEye Dec 07 '21 at 12:20
  • Could you please also share the code that iteratively tweaks the bands and desired parameters? Many thanks! – John Apr 18 '23 at 08:06
  • @John I didn't use code to iteratively tweak the bands at the time. It was more of a try something, look at the result and manually tweak. If a response in a band was too big, I then adjusted the `desired` value to something a lower by the band error. The width of the middle section was also manually adjusted to get the ripples at the corners of the transition band. Presumably that could be done in code as well, but I haven't considered any implications of an 'optimal' tweaking sequence between band width and desired band values. – SleuthEye Apr 20 '23 at 00:21