Commit 4ef0081f authored by Jean Ibarz's avatar Jean Ibarz
Browse files

Added a script to generate Gammatone FIR filters in .wav format, by using the...

Added a script to generate Gammatone FIR filters in .wav format, by using the library from Bingo Todd.
Because a rectangular windowing produces ringing in the magnitude of the frequency response at low frequencies when the filter have ~512 coefficients, we apply a windowing of the filter impulse response in order to improve the filter quality. After windowing, a normalization is applied to restore the correct magnitude peak in the frequency response domain.
parent 13f8183c
import os
import GTF.GTF as gtf_proposed
import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft
from scipy.io import wavfile
from scipy.signal import windows
"""
This script generates N_BAND passband filters, in particular Gammatone filters, by using a library from Bingo Todd,
available here: https://github.com/bingo-todd/Gammatone-filters/blob/master/examples/validate.py
"""
SAMPLING_RATE = 44100
RAW_FILTER_LENGTH = 512
N_BAND = 20
APPLY_WINDOWING = True
WINDOWED_FILTER_LENGTH = 512
SHOW_PLOT = True
"""
If not already existing, create a folder to store the generated filters as .wav files.
"""
GTF_DIR = '../gtf/'
if not os.path.exists(GTF_DIR):
os.makedirs(GTF_DIR)
def generate_gtf():
fs = SAMPLING_RATE
# raw impulse length
x_len = int(fs)
x = np.zeros(shape=(RAW_FILTER_LENGTH)) # filters shape = (44100 coefficients, 1 channel)
x[0] = 1
gtf_obj = gtf_proposed(fs, cf_low=100, cf_high=20000, n_band=N_BAND)
cfs = gtf_obj.cfs
bws = gtf_obj.cal_bw(cfs)
irs = gtf_obj.filter(x)
if APPLY_WINDOWING is True:
# at low frequencies, filter length must be important
# in order to avoid side effects from rectangular truncature,
# we want to apply a Hann window on the right half side of the impulse,
filter_length = 512
hann_half_window = windows.hann(M=2 * filter_length, sym=True)[filter_length:]
# fig = gtf_obj.plot_ir_spec(irs[:, :1000])
# plt.show()
win_irs = np.array([np.multiply(ir[0:filter_length], np.expand_dims(hann_half_window, axis=-1)) for ir in irs])
# frequency magnitude of filters is affected by windowing
# we want to normalize the frequency magnitude maximum to 0dB
# we compute the mag of the fft of each ir, and normalize the
# filter impulse response accordingly.
for i in range(N_BAND):
# calculate the magnitude of the frequency response corresponding to the impulse response right zero padded.
# Zero padding in the time domain correspond to an interpolation in the frequency domain. It is done in
# order to increase the resolution of the frequency response of the filter, in order to get the best
# possible estimate of the maximal magnitude.
mag = np.abs(fft(np.concatenate([win_irs[i, :].flatten(), np.zeros(shape=(60000,))])))
mag_max = np.max(mag)
win_irs[i, :] = win_irs[i, :] / mag_max
raw_irs = np.array(win_irs)
if SHOW_PLOT:
i_filter = -1
plt.plot(irs[i_filter] / np.max(irs[i_filter]), label='Raw Filter IR')
plt.plot(hann_half_window, label='Hann right half window')
plt.plot(win_irs[i_filter] / np.max(win_irs[i_filter]), label='Windowed Filter IR')
plt.legend()
plt.xlim([0, WINDOWED_FILTER_LENGTH])
plt.yscale('log')
plt.ylim([10 ** (-80 / 20), 2]) # lower limit for amplitude value: -60dB
plt.show()
# export FIR into .wav files
for i in range(N_BAND):
filename = f"filter_{i}_CFS{int(cfs[i])}_BWS{int(bws[i])}_L{filter_length}_W{'hann'}.wav"
if APPLY_WINDOWING is True:
ir = win_irs[i, :, 0]
else:
ir = raw_irs[i, :, 0]
wavfile.write(filename=os.path.join(GTF_DIR, filename), rate=SAMPLING_RATE, data=ir)
assert os.path.exists(os.path.join(GTF_DIR, filename))
if __name__ == '__main__':
generate_gtf()
print(f'{N_BAND} Gammatone filters generated successfully in path {GTF_DIR}')
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment