Previous: 基本的な定数などの計算・変換 Up: 音の処理のためのPythonファイル Next: 信号生成 / Signal Generation
##### FILTER ##################################################################
# biquad.lpf = function(fc, samprate) {
# K = tan(pi * fc / samprate)
# b0 = K**2 / (1 + sqrt(2)*K + K**2)
# b1 = 2*K**2 / (1 + sqrt(2)*K + K**2)
# b2 = K**2 / (1 + sqrt(2)*K + K**2)
# a0 = 1
# a1 = 2*(K**2-1) / (1 + sqrt(2)*K + K**2)
# a2 = (1 - sqrt(2)*K + K**2) / (1 + sqrt(2)*K + K**2)
# filt = Arma(c(b0, b1, b2), c(a0, a1, a2))
# return(filt)
# }
# biquad.hpf = function(fc, samprate) {
# K = tan(pi * fc / samprate)
# b0 = 1 / (1 + sqrt(2)*K + K**2)
# b1 = -2 / (1 + sqrt(2)*K + K**2)
# b2 = 1 / (1 + sqrt(2)*K + K**2)
# a0 = 1
# a1 = 2*(K**2-1) / (1 + sqrt(2)*K + K**2)
# a2 = (1 - sqrt(2)*K + K**2) / (1 + sqrt(2)*K + K**2)
# filt = Arma(c(b0, b1, b2), c(a0, a1, a2))
# return filt
# }
# biquad.lowshelf = function(fc, G, samprate) {
# K = tan(pi * fc / samprate)
# if (G >= 0.0) {
# V0 = 10 ** (G/20)
# b0 = (1 + sqrt(2*V0)*K + V0*K**2) / (1 + sqrt(2)*K + K**2)
# b1 = 2*(V0*K**2 - 1) / (1 + sqrt(2)*K + K**2)
# b2 = (1 - sqrt(2*V0)*K + V0*K**2) / (1 + sqrt(2)*K + K**2)
# a0 = 1
# a1 = 2*(K**2 - 1) / (1 + sqrt(2)*K + K**2)
# a2 = (1 - sqrt(2)*K + K**2) / (1 + sqrt(2)*K + K**2)
# } else {
# V0 = 10 ** (-G/20)
# b0 = (1 + sqrt(2)*K + K**2) / (1 + sqrt(2*V0)*K + V0*K**2)
# b1 = 2*(K**2 - 1) / (1 + sqrt(2*V0)*K + V0*K**2)
# b2 = (1 - sqrt(2)*K + K**2) / (1 + sqrt(2*V0)*K + V0*K**2)
# a0 = 1
# a1 = 2*(V0*K**2 - 1) / (1 + sqrt(2*V0)*K + V0*K**2)
# a2 = (1 - sqrt(2*V0)*K + V0*K**2) / (1 + sqrt(2*V0)*K + V0*K**2)
# }
# filt = Arma(c(b0, b1, b2), c(a0, a1, a2))
# return filt
# }
# biquad.highshelf = function(fc, G, samprate) {
# K = tan(pi * fc / samprate)
# if (G >= 0.0) {
# V0 = 10 ** (G/20)
# b0 = (V0 + sqrt(2*V0)*K + K**2) / (1 + sqrt(2)*K + K**2)
# b1 = 2*(K**2 - V0) / (1 + sqrt(2)*K + K**2)
# b2 = (V0 - sqrt(2*V0)*K + K**2) / (1 + sqrt(2)*K + K**2)
# a0 = 1
# a1 = 2*(K**2 - 1) / (1 + sqrt(2)*K + K**2)
# a2 = (1 - sqrt(2)*K + K**2) / (1 + sqrt(2)*K + K**2)
# } else {
# V0 = 10 ** (-G/20)
# b0 = (1 + sqrt(2)*K + K**2) / (V0 + sqrt(2*V0)*K + K**2)
# b1 = 2*(K**2 - 1) / (V0 + sqrt(2*V0)*K + K**2)
# b2 = (1 - sqrt(2)*K + K**2) / (V0 + sqrt(2*V0)*K + K**2)
# a0 = 1
# a1 = 2*(K**2/V0 - 1) / (1 + sqrt(2/V0)*K + K**2/V0)
# a2 = (1 - sqrt(2/V0)*K + K**2/V0) / (1 + sqrt(2/V0)*K + K**2/V0)
# }
# filt = Arma(c(b0, b1, b2), c(a0, a1, a2))
# return filt
# }
# biquad.peak = function(fc, G, Q, samprate) {
# K = tan(pi * fc / samprate)
# if (G >= 0.0) {
# V0 = 10 ** (G/20)
# b0 = (1 + V0/Q*K + K**2) / (1 + 1/Q*K + K**2)
# b1 = 2*(K**2 - 1) / (1 + 1/Q*K + K**2)
# b2 = (1 - V0/Q*K + K**2) / (1 + 1/Q*K + K**2)
# a0 = 1
# a1 = 2*(K**2 - 1) / (1 + 1/Q*K + K**2)
# a2 = (1 - 1/Q*K + K**2) / (1 + 1/Q*K + K**2)
# } else {
# V0 = 10 ** (-G/20)
# b0 = (1 + 1/Q*K + K**2) / (1 + V0/Q*K + K**2)
# b1 = 2*(K**2 - 1) / (1 + V0/Q*K + K**2)
# b2 = (1 - 1/Q*K + K**2) / (1 + V0/Q*K + K**2)
# a0 = 1
# a1 = 2*(K**2 - 1) / (1 + V0/Q*K + K**2)
# a2 = (1 - V0/Q*K + K**2) / (1 + V0/Q*K + K**2)
# }
# filt = Arma(c(b0, b1, b2), c(a0, a1, a2))
# return filt
# }
def filterbank(x, samprate, oct=1):
”'
FILTERBANK octave filterbank using raised cosine frequency window
out, cfreqs = filterbank(x, samprate) create filtered signals `out' from x.
Columns of `out' is a filtered signal with corresponding center
frequency in `cfreqs'. By default, 1/1-octave band filters are used.
Filtered signals at the lowest frequency is processed with low-pass,
the highest frequency is processed with high-pass, and all frequencies
in between are processed with band-pass filters. This yields to
(almost) perfect reconstruction of the original signal, such that:
out, cfreqs = filterbank(x, samprate)
z = sum(out, 1)
creates z almost the same as x.
out, cfreqs = filterbank(x, samprate, oct) creates and filters using
1/oct-octave band filters.
”'
# prepare center frequencies
cfreqs = 1000 * 2**(np.arange(-100, 100)/oct)
cfreqs = cfreqs[np.logical_and(20<=cfreqs, cfreqs<=samprate/2)]
# create octave-band bandpass filter
if len(x) % 2 == 1:
x = np.concatenate((x, [0]))
out = np.zeros((len(x), len(cfreqs)))
xx = scipy.fft(x)
for n in range(len(cfreqs)):
fc = cfreqs[n]
fl = fc * 2**(-1/oct)
fh = fc * 2**(+1/oct)
W = np.zeros(len(xx)//2)
f = np.linspace(0, samprate/2, len(W))
ff = np.logical_and(fl<=f, f<=fh)
cc = (np.cos(np.log2(f[ff]/fc)*math.pi*oct) + 1) / 2
W[ff] = cc
if n==0:
W[f < fc] = 1
elif n==len(cfreqs)-1:
W[fc < f] = 1
W = np.concatenate((W, np.flipud(W)))
W[0] = 0
W[len(W)//2] = 0
# apply filter
y = np.real(scipy.ifft(xx * W))
out[:,n] = y
return(out, cfreqs)