1、Gammatone简介
Gammatone滤波器被广泛用于模拟人类听觉系统对信号的处理方式，作为语音信号的一类听觉分析滤波器（以下简称为GT滤波器）。GT滤波器只需要很少的参数就能很好地模拟听觉实验中的生理数据，能够体现基底膜尖锐的滤波特性，而且 GT滤波器具有简单的冲激响应函数，能够由此推导出GT函数的传递函数，进行各种滤波器性能分析，同时有利于听觉模型的电路实现。
GT滤波器的冲击响应函数定义如下：
这里n为滤波器阶数，b为滤波器的带宽，f为滤波器的中心频率，a是振幅。其等于
a=B**n 这里 B=b1ERB(f)
ERB(f)为GT滤波器的等价矩形带宽（等价矩形带宽：对于同样的白噪声输入，和指定的滤波器通过一样能量的矩形滤波器的宽度，简称ERB)，它同GT滤波器中心频率f，的关系是
ERB（f）=24.7+0.108f
b1=1.019是为了让GT函数更好地与生理数据相符而引入的参数。
GT滤波器的幅频特性如下图所示
2、python代码
以下的python代码分别实现了频域Gammtone滤波和时域Gammtone滤波器，频域滤波时参数USEFFT=1。
from Hammingstft import Hstft2,Hstft1
from ERBFilterBank import ERBFilterBank
from MakeERBFilters import MakeERBFilters
from fft2gammatonemx import fft2gammatonemx
import matplotlib.mlab
import numpy as np
import math
import scipy
import pylab
def FFtGamatone(X,SR=16000.0,TWIN=0.025,THOP=0.010,N=64.0,FMIN=50.0,FMAX=0,USEFFT=1,WIDTH=1.0):
if FMAX==0 or FMAX
其中MakeERBFilters函数产生了GT滤波器系数。输入为：采样频率、滤波器通道数、最小频率。输出为：n通道GT滤波器系数。代码如下
import numpy as np
from ERPSpace import ERBSpace
import pylab
def MakeERBFilters(fs,numChannels,lowFreq):
fs=float(fs)
lowFreq=float(lowFreq)
T = 1/fs
if np.isscalar(numChannels):
numChannels=[numChannels]
numChannels=np.array(numChannels)
if numChannels.ndim >1:
print "is not one dimision data"
return 0
if numChannels.size == 1:
cf = ERBSpace(lowFreq, fs/2.0, numChannels)
else:
cf = numChannels;
EarQ = 9.26449 # Glasberg and Moore Parameters
minBW = 24.7
order = 1.0
ERB = ((cf/EarQ)**order + minBW**order)**(1/order)
B=1.019*2.0*np.pi*ERB
A0 = T
A2 = 0.0
B0 = 1.0
B1 = -2*np.cos(2.0*cf*np.pi*T)/np.exp(B*T)
B2 = np.exp(-2.0*B*T)
A11 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) + 2.0*np.sqrt(3.0+2.0**1.5)*T*np.sin(2.0*cf*np.pi*T)/np.exp(B*T))/2.0
A12 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) - 2.0*np.sqrt(3.0+2.0**1.5)*T*np.sin(2.0*cf*np.pi*T)/np.exp(B*T))/2.0
A13 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) + 2.0*np.sqrt(3.0-2.0**1.5)*T*np.sin(2.0*cf*np.pi*T)/np.exp(B*T))/2.0
A14 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) - 2.0*np.sqrt(3.0-2.0**1.5)*T*np.sin(2.0*cf*np.pi*T)/np.exp(B*T))/2.0
gain = np.abs((-2.0*np.exp(4.0j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2.0j*cf*np.pi*T)*T* \
(np.cos(2.0*cf*np.pi*T) - np.sqrt(3.0 - 2.0**(3.0/2.0))* np.sin(2.0*cf*np.pi*T))) * \
(-2*np.exp(4.0j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2.0j*cf*np.pi*T)*T* \
(np.cos(2.0*cf*np.pi*T) + np.sqrt(3.0 - 2.0**(3.0/2.0)) * np.sin(2.0*cf*np.pi*T)))* \
(-2.0*np.exp(4.0j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2.0j*cf*np.pi*T)*T* \
(np.cos(2.0*cf*np.pi*T) - np.sqrt(3.0 + 2.0**(3.0/2.0))*np.sin(2.0*cf*np.pi*T))) * \
(-2.0*np.exp(4.0j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2.0j*cf*np.pi*T)*T* \
(np.cos(2.0*cf*np.pi*T) + np.sqrt(3.0 + 2.0**(3.0/2.0))*np.sin(2.0*cf*np.pi*T)))/(-2.0 / np.exp(2.0*B*T) \
- 2.0*np.exp(4.0j*cf*np.pi*T) + 2.0*(1.0 + np.exp(4.0j*cf*np.pi*T))/np.exp(B*T))**4.0);
allfilts = np.ones((cf.size,1),dtype=np.float64)
fcoefs=np.zeros((cf.size,10),dtype=np.float64)
fcoefs[:] = zip(A0*allfilts, A11, A12, A13, A14, A2*allfilts, B0*allfilts, B1, B2, gain)
return fcoefs
# Testing
if __name__ == '__main__':
chann=np.array([64])
cfs = MakeERBFilters(22050,chann,100);
print cfs
pylab.figure()
pylab.imshow(cfs, origin='lower', aspect='auto', interpolation='nearest')
pylab.ylabel('Frequency')
pylab.show()
代码中其中ERBSpace获得了中心频率f。输入为：最小频率、最大频率、滤波器通道数。输出为滤波器族的各个中心频率。代码如下：
import numpy as np
from matplotlib import pyplot
def ERBSpace(lowFreq = 100.0, highFreq = 44100.0/4.0, N=100.0):
earQ = 9.26449
minBW = 24.7
low = float(lowFreq)
high = float(highFreq)
N=float(N)
cf = -(earQ * minBW) + np.exp((np.arange(N+1)[1:]) * (-np.log(high + earQ * minBW) + np.log(low + earQ * minBW)) / (N)) * (high + earQ * minBW)
#cf = cf[::-1]
return cf
# Testing
if __name__ == '__main__':
cf = ERBSpace()
print cf
pyplot.plot(cf)
pyplot.show()
ERBFilterBank 函数输入分别为：原始数据和GT滤波器系数。输出为滤波后的数据。该函数实现对原始数据的时域GT滤波。代码如下：
from MakeERBFilters import MakeERBFilters
from scipy import signal
import numpy as np
import pylab
def ERBFilterBank(x, fcoefs=0.0):
if isinstance(fcoefs,np.ndarray)!=True:
fcoefs = MakeERBFilters(22050,64,100)
if fcoefs.ndim!=2:
print "fcoefs parameter passed to ERBFilterBank is the wrong size."
col=fcoefs.shape[1]
if col!=10:
print "fcoefs parameter passed to ERBFilterBank is the wrong size."
A0 = fcoefs[:,0];
A11 = fcoefs[:,1];
A12 = fcoefs[:,2];
A13 = fcoefs[:,3];
A14 = fcoefs[:,4];
A2 = fcoefs[:,5];
B0 = fcoefs[:,6];
B1 = fcoefs[:,7];
B2 = fcoefs[:,8];
gain= fcoefs[:,9];
output = np.zeros((gain.size, x.size),dtype=np.float64)
for chan in range(gain.size):
y1=signal.lfilter(np.array([A0[chan]/gain[chan],A11[chan]/gain[chan],A2[chan]/gain[chan]]), np.array([B0[chan],B1[chan],B2[chan]]), x)
y2=signal.lfilter(np.array([A0[chan],A12[chan],A2[chan]]), np.array([B0[chan],B1[chan],B2[chan]]), y1)
y3=signal.lfilter(np.array([A0[chan],A13[chan],A2[chan]]), np.array([B0[chan],B1[chan],B2[chan]]), y2)
y4=signal.lfilter(np.array([A0[chan],A14[chan],A2[chan]]),np.array([B0[chan],B1[chan],B2[chan]]), y3)
output[chan, :] = y4;
return output
# Testing
if __name__ == '__main__':
xx=np.ones((1,100))
cfs = ERBFilterBank(xx);
print cfs
pylab.figure()
pylab.imshow(cfs, origin='lower', aspect='auto', interpolation='nearest')
pylab.ylabel('Frequency')
pylab.show()
fft2gammatonemx为获得了GT滤波器的频率域系数。输入为：滤波器长度、采样频率、滤波器通道数、滤波器带宽、最小频率、最大频率、最大长度。输出为GT滤波器的频率域系数。代码如下：
import numpy as np
from ERPSpace import ERBSpace
import pylab
def fft2gammatonemx(nfft, sr=16000, nfilts=64, width=1.0, minfreq=100, *othermore):
sr=float(sr)
nfilts=float(nfilts)
minfreq=float(minfreq)
EarQ = 9.26449
minBW = 24.7
order = 1.0
if len(othermore)==0:
maxfreq=sr/2
maxlen=nfft
elif len(othermore)==1:
maxfreq=float(othermore[0])
maxlen=nfft
elif len(othermore)==2:
maxfreq=float(othermore[0])
maxlen=othermore[1]
else:
print "param error"
return 0
wts = np.zeros((nfilts, nfft),dtype=np.float64)
cfreqs = ERBSpace(minfreq, maxfreq, nfilts)
cfreqs = np.flipud(cfreqs)
GTord = 4
ucirc = np.exp(2j*np.pi*np.arange(nfft/2+1)/nfft)
justpoles = 0
for i in range(int(nfilts)):
cf = cfreqs[i];
ERB = width*((cf/EarQ)**order + minBW**order)**(1/order);
B = 1.019*2*np.pi*ERB
r = np.exp(-B/sr)
theta = 2*np.pi*cf/sr
pole = r*np.exp(1j*theta)
if justpoles == 1:
cosomegamax = (1.0+r*r)/(2.0*r)*np.cos(theta);
if np.abs(cosomegamax) > 1:
if theta < np.pi/2.0:
omegamax = 0
else:
omegamax =np.pi
else:
omegamax = np.arccos(cosomegamax)
center = np.exp(1j*omegamax)
gain = np.abs((pole-center)*(pole.conjugate()-center))**GTord
wts[i,1:(nfft/2+1)] = gain * (np.abs((pole-ucirc)*(pole.conjugate()-ucirc))**(-GTord))
else:
T = 1/sr;
A11 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) + 2.0*np.sqrt(3.0+2.0**1.5)*T*np.sin(2*cf*np.pi*T)/np.exp(B*T))/2.0
A12 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) - 2.0*np.sqrt(3.0+2.0**1.5)*T*np.sin(2*cf*np.pi*T)/np.exp(B*T))/2.0
A13 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) + 2.0*np.sqrt(3.0-2.0**1.5)*T*np.sin(2*cf*np.pi*T)/np.exp(B*T))/2.0
A14 = -(2.0*T*np.cos(2.0*cf*np.pi*T)/np.exp(B*T) - 2.0*np.sqrt(3.0-2.0**1.5)*T*np.sin(2*cf*np.pi*T)/np.exp(B*T))/2.0
zros = -np.array([A11,A12,A13,A14])/T
gain = np.abs((-2.0*np.exp(4j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2j*cf*np.pi*T)*T* \
(np.cos(2.0*cf*np.pi*T) - np.sqrt(3.0 - 2.0**(3.0/2.0))* np.sin(2.0*cf*np.pi*T)))* \
(-2.0*np.exp(4j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2j*cf*np.pi*T)*T* \
(np.cos(2.0*cf*np.pi*T) + np.sqrt(3.0 - 2.0**(3.0/2.0))*np.sin(2.0*cf*np.pi*T)))* \
(-2.0*np.exp(4j*cf*np.pi*T)*T + 2.0*np.exp(-(B*T) + 2j*cf*np.pi*T)*T*(np.cos(2*cf*np.pi*T) - \
np.sqrt(3.0 + 2.0**(3.0/2.0))*np.sin(2.0*cf*np.pi*T)))*(-2.0*np.exp(4j*cf*np.pi*T)*T + \
2.0*np.exp(-(B*T) + 2j*cf*np.pi*T)*T*(np.cos(2*cf*np.pi*T) + np.sqrt(3.0 + 2.0**(3.0/2.0))* \
np.sin(2.0*cf*np.pi*T)))/(-2.0/np.exp(2*B*T) - 2.0*np.exp(4j*cf*np.pi*T)+2.0*(1 + np.exp(4j*cf*np.pi*T))/np.exp(B*T))**4)
wts[i,np.arange(0,int(nfft/2+1))] = ((T**4)/gain)*np.abs(ucirc-zros[0])*np.abs(ucirc-zros[1])*np.abs(ucirc-zros[2])*np.abs(ucirc-zros[3]) \
*(np.abs((pole-ucirc)*(pole.conjugate()-ucirc))**(-GTord))
wts = wts[:,0:maxlen]
return wts
# Testing
if __name__ == '__main__':
cfs = fft2gammatonemx(256, 16000, 64,1,50,8000, 1024.0/2.0+1);
print cfs
pylab.figure()
pylab.imshow(cfs, origin='lower', aspect='auto', interpolation='nearest')
pylab.ylabel('Frequency')
pylab.show()
测试效果图：
360
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
