Skip to content

Commit

Permalink
[python] Add pylint config
Browse files Browse the repository at this point in the history
  • Loading branch information
tobiashienzsch committed Sep 1, 2024
1 parent 6aea7d6 commit f33caca
Show file tree
Hide file tree
Showing 27 changed files with 1,283 additions and 648 deletions.
641 changes: 641 additions & 0 deletions .pylintrc

Large diffs are not rendered by default.

1 change: 0 additions & 1 deletion src/python/pffdtd/air_abs/get_air_absorption.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
# SPDX-License-Identifier: MIT

import numpy as np
from numpy import array as npa
from numpy import log10, exp, sqrt, log, pi

def get_air_absorption(freq_vec,temperature_celsius,rel_humidity_pnct,pressure_atmospheric_kPa=101.325):
Expand Down
58 changes: 31 additions & 27 deletions src/python/pffdtd/air_abs/modal_filter.py
Original file line number Diff line number Diff line change
@@ -1,19 +1,16 @@
# SPDX-License-Identifier: MIT

import numpy as np
import numba as nb
from numpy import array as npa
from numpy import exp, sqrt, log, pi, cos
from scipy.fft import dct,idct #default type2
from pffdtd.air_abs.get_air_absorption import get_air_absorption
from pffdtd.common.myfuncs import iround, iceil
import numpy as np
from numpy import exp, pi, cos
from scipy.fft import dct, idct # default type2
from tqdm import tqdm

#apply filter, x is (Nchannel,Nsamples) array
#see I3DA paper for details
#Tc is temperature in deg Celsius
#rh is relative humidity
def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0):
from pffdtd.air_abs.get_air_absorption import get_air_absorption
from pffdtd.common.myfuncs import iceil


def apply_modal_filter(x, Fs, Tc, rh, pad_t=0.0):
"""
This is an implementation of an air absorption filter based on a modal approach.
It solves a system of 1-d dissipative wave equations tune to air attenuation
Expand All @@ -23,13 +20,19 @@ def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0):
Hamilton, B. "Adding air attenuation to simulated room impulse responses: A
modal approach", to be presented at I3DA 2021 conference in Bologna Italy.
"""

# apply filter, x is (Nchannel,Nsamples) array
# see I3DA paper for details
# Tc is temperature in deg Celsius
# rh is relative humidity

Ts = 1/Fs

x = np.atleast_2d(x)
Nt0 = x.shape[-1]
Nt = iceil(pad_t/Ts)+Nt0
xp = np.zeros((x.shape[0],Nt))
xp[:,:Nt0] = x
xp = np.zeros((x.shape[0], Nt))
xp[:, :Nt0] = x
del x

y = np.zeros(xp.shape)
Expand All @@ -38,40 +41,41 @@ def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0):
wqTs = pi*(np.arange(Nx)/Nx)
wq = wqTs/Ts

rd = get_air_absorption(wq/2/pi,Tc,rh)
rd = get_air_absorption(wq/2/pi, Tc, rh)
alphaq = rd['absfull_Np']
c = rd['c']

P0 = np.zeros(xp.shape)
P1 = np.zeros(xp.shape)

fx = np.zeros(xp.shape)
fx[:,0] = 1
Fm = dct(fx,type=2,norm='ortho',axis=-1)
fx[:, 0] = 1
Fm = dct(fx, type=2, norm='ortho', axis=-1)

sigqTs = c*alphaq*Ts
a1 = 2*exp(-sigqTs)*cos(wqTs)
a2 = -exp(-2*sigqTs)
Fmsig1 = Fm*(1+sigqTs/2)/(1+sigqTs)
Fmsig2 = Fm*(1-sigqTs/2)/(1+sigqTs)

u = np.zeros((xp.shape[0],Nt+1))
u[:,1:] = xp[:,::-1] #flip
u = np.zeros((xp.shape[0], Nt+1))
u[:, 1:] = xp[:, ::-1] # flip

pbar = tqdm(total=Nt,desc=f'modal filter',ascii=True)
pbar = tqdm(total=Nt, desc='modal filter', ascii=True)

@nb.jit(nopython=True,parallel=True)
def _run_step(P0,P1,a1,a2,Fmsig1,Fmsig2,un1,un0):
@nb.jit(nopython=True, parallel=True)
def _run_step(P0, P1, a1, a2, Fmsig1, Fmsig2, un1, un0):
P0[:] = a1*P1 + a2*P0 + Fmsig1*un1 - Fmsig2*un0

for n in range(Nt):
#P0 = a1*P1 + a2*P0 + Fmsig1*u[:,n+1] - Fmsig2*u[:,n]
# P0 = a1*P1 + a2*P0 + Fmsig1*u[:,n+1] - Fmsig2*u[:,n]
for i in range(P0.shape[0]):
_run_step(P0[i],P1[i],a1,a2,Fmsig1[i],Fmsig2[i],u[i,n+1],u[i,n])
if n<Nt-1: #dont swap on last sample
P1,P0 = P0,P1
_run_step(P0[i], P1[i], a1, a2, Fmsig1[i],
Fmsig2[i], u[i, n+1], u[i, n])
if n < Nt-1: # dont swap on last sample
P1, P0 = P0, P1
pbar.update(1)
pbar.close()

y = idct(P0,type=2,norm='ortho',axis=-1)
return np.squeeze(y) #squeeze to 1d in case
y = idct(P0, type=2, norm='ortho', axis=-1)
return np.squeeze(y) # squeeze to 1d in case
42 changes: 21 additions & 21 deletions src/python/pffdtd/air_abs/ola_filter.py
Original file line number Diff line number Diff line change
@@ -1,18 +1,18 @@
# SPDX-License-Identifier: MIT

import numpy as np
import numba as nb
from numpy import array as npa
from numpy import exp, sqrt, log2, pi, ceil, cos
from scipy.fft import rfft,irfft
from numpy import exp, log2, pi, ceil, cos
from scipy.fft import rfft, irfft
from pffdtd.air_abs.get_air_absorption import get_air_absorption
from pffdtd.common.myfuncs import iround, iceil
from tqdm import tqdm

#apply the filter, x is (Nchannels,Nsamples) array
#Tc is temperature degrees Celsius
#rh is relative humidity
def apply_ola_filter(x,Fs,Tc,rh,Nw=1024):
# apply the filter, x is (Nchannels,Nsamples) array
# Tc is temperature degrees Celsius
# rh is relative humidity


def apply_ola_filter(x, Fs, Tc, rh, Nw=1024):
"""
This is an implementation of overlap-add (STFT/iSTFT) air absorption filtering.
Tuned for 75% overlap and 1024-sample Hann window at 48kHz.
Expand All @@ -35,35 +35,35 @@ def apply_ola_filter(x,Fs,Tc,rh,Nw=1024):
assert Np < Nw
Nfft_h = np.int_(Nfft//2+1)

xp = np.zeros((x.shape[0],Nw+Nt0+Np))
xp[:,Nw:Nw+Nt0] = x
y = np.zeros((x.shape[0],Nt0+Np))
xp = np.zeros((x.shape[0], Nw+Nt0+Np))
xp[:, Nw:Nw+Nt0] = x
y = np.zeros((x.shape[0], Nt0+Np))
del x

wa = 0.5*(1-cos(2*pi*np.arange(Nw)/Nw)) #hann window
ws = wa/(3/8*Nw/Ha) #scaled for COLA
wa = 0.5*(1-cos(2*pi*np.arange(Nw)/Nw)) # hann window
ws = wa/(3/8*Nw/Ha) # scaled for COLA

fv = np.arange(Nfft_h)/Nfft*Fs
rd = get_air_absorption(fv,Tc,rh)
rd = get_air_absorption(fv, Tc, rh)
c = rd['c']
absNp = rd['absfull_Np']

for i in range(xp.shape[0]):
pbar = tqdm(total=NF,desc=f'OLA filter channel {i}',ascii=True)
pbar = tqdm(total=NF, desc=f'OLA filter channel {i}', ascii=True)
yp = np.zeros((xp.shape[-1],))
for m in range(NF):
na0 = m*Ha
assert na0+Nw<=Nw+Nt0+Np
assert na0+Nw <= Nw+Nt0+Np
dist = c*Ts*(na0-Nw/2)
xf = xp[i,na0:na0+Nw]
if dist<0: #dont apply gain (negative times - pre-padding)
xf = xp[i, na0:na0+Nw]
if dist < 0: # dont apply gain (negative times - pre-padding)
yp[na0:na0+Nw] += ws*xf
else:
Yf = rfft(wa*xf,Nfft)*exp(-absNp*dist)
yf = irfft(Yf,Nfft)[:Nw]
Yf = rfft(wa*xf, Nfft)*exp(-absNp*dist)
yf = irfft(Yf, Nfft)[:Nw]
yp[na0:na0+Nw] += ws*yf

pbar.update(1)
y[i] = yp[Nw:]
pbar.close()
return np.squeeze(y) #squeeze to 1d in case
return np.squeeze(y) # squeeze to 1d in case
44 changes: 21 additions & 23 deletions src/python/pffdtd/air_abs/test_air_abs_filters.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,58 +13,56 @@
##############################################################################

import numpy as np
import numba as nb
from numpy import array as npa
from numpy import exp, sqrt, log, pi
from numpy import exp, log
from numpy.random import random_sample
import matplotlib.pyplot as plt
from pffdtd.air_abs.get_air_absorption import get_air_absorption

from pffdtd.air_abs.visco_filter import apply_visco_filter
from pffdtd.air_abs.modal_filter import apply_modal_filter
from pffdtd.air_abs.ola_filter import apply_ola_filter
from numpy.random import random_sample
from pffdtd.common.myfuncs import iround, iceil
from pffdtd.common.myfuncs import iround

Tc = 20
rh = 60

SR = 48e3

#single channel test
# single channel test
t_end = 1
Nt0 = iround(t_end*SR)
tx = np.arange(Nt0)/SR
td = 0.02
tau = t_end/(6*log(10))
x = exp(-(tx-td)/tau)*(2*random_sample(Nt0)-1)
x[:int(td*SR)]=0
x[:int(td*SR)] = 0

#multi channel test
#x = np.zeros((2,iround(1.2*SR)))
#nstart = iround(0.01*SR)
#nskip = iround(0.005*SR)
#x[0,nstart::nskip] = 1.0
#x[1,nstart::nskip] = -1.0
#tx = np.arange(x.shape[-1])/SR
# multi channel test
# x = np.zeros((2,iround(1.2*SR)))
# nstart = iround(0.01*SR)
# nskip = iround(0.005*SR)
# x[0,nstart::nskip] = 1.0
# x[1,nstart::nskip] = -1.0
# tx = np.arange(x.shape[-1])/SR

y1 = apply_visco_filter(x, SR, Tc, rh)
y2 = apply_modal_filter(x, SR, Tc, rh)
y3 = apply_ola_filter(x, SR, Tc, rh)

ty1 = np.arange(0,y1.shape[-1])/SR
ty2 = np.arange(0,y2.shape[-1])/SR
ty3 = np.arange(0,y3.shape[-1])/SR
ty1 = np.arange(0, y1.shape[-1])/SR
ty2 = np.arange(0, y2.shape[-1])/SR
ty3 = np.arange(0, y3.shape[-1])/SR

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(tx,x.T,linestyle='-',color='b',label='orig')
ax.plot(ty1,y1.T,linestyle='-',color='g',label='stokes')
ax.plot(ty2,y2.T,linestyle='-',color='r',label='modal')
ax.plot(ty3,y3.T,linestyle='-',color='y',label='OLA')
ax.plot(tx, x.T, linestyle='-', color='b', label='orig')
ax.plot(ty1, y1.T, linestyle='-', color='g', label='stokes')
ax.plot(ty2, y2.T, linestyle='-', color='r', label='modal')
ax.plot(ty3, y3.T, linestyle='-', color='y', label='OLA')
ax.margins(0, 0.1)
ax.grid(which='both', axis='both')
ax.legend()
plt.show()


#multi channel test
# multi channel test
SR = 48e3
40 changes: 20 additions & 20 deletions src/python/pffdtd/air_abs/visco_filter.py
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
# SPDX-License-Identifier: MIT

import numpy as np
import numba as nb
from numpy import array as npa
from numpy import exp, sqrt, log, pi
from pffdtd.common.myfuncs import iround, iceil
from pffdtd.common.myfuncs import iceil
from pffdtd.air_abs.get_air_absorption import get_air_absorption
from tqdm import tqdm

#function to apply filter, main input being x, np.ndarray (Nchannels,Nsamples)
#enter temperature (Tc) and relative humidity (rh)
#NdB should be above 60dB, that is for truncation of Gaussian kernel
def apply_visco_filter(x,Fs,Tc,rh,NdB=120,t_start=None):
# function to apply filter, main input being x, np.ndarray (Nchannels,Nsamples)
# enter temperature (Tc) and relative humidity (rh)
# NdB should be above 60dB, that is for truncation of Gaussian kernel


def apply_visco_filter(x, Fs, Tc, rh, NdB=120, t_start=None):
"""
This is an implementation of an air absorption filter based on approximate
Green's function Stoke's equation (viscothermal wave equation)
Expand All @@ -20,14 +20,13 @@ def apply_visco_filter(x,Fs,Tc,rh,NdB=120,t_start=None):
Hamilton, B. "Air absorption filtering method based on approximate Green's
function for Stokes' equation", to be presented at the DAFx2021 e-conference.
"""
rd = get_air_absorption(1,Tc,rh)
c = rd['c']
rd = get_air_absorption(1, Tc, rh)
g = rd['gamma_p']

Ts = 1/Fs
if t_start is None:
t_start = Ts**2/(2*pi*g)
print(f'{t_start=}')
t_start = Ts**2/(2*pi*g)
print(f'{t_start=}')

x = np.atleast_2d(x)
Nt0 = x.shape[-1]
Expand All @@ -36,23 +35,24 @@ def apply_visco_filter(x,Fs,Tc,rh,NdB=120,t_start=None):
dt_end = Fs*sqrt(0.1*log(10)*NdB*n_last*Ts*g)
Nt = Nt0+iceil(dt_end)

y = np.zeros((x.shape[0],Nt))
y = np.zeros((x.shape[0], Nt))
n_start = iceil(t_start*Fs)
assert(n_start>0)
assert n_start > 0

y[:,:n_start] = x[:,:n_start]
y[:, :n_start] = x[:, :n_start]
Tsg2 = 2*Ts*g
Tsg2pi = 2*Ts*g*pi
gTs = g*Ts
dt_fac = 0.1*log(10)*NdB*gTs
pbar = tqdm(total=Nt,desc=f'visco filter',ascii=True)
for n in range(n_start,Nt0):
pbar = tqdm(total=Nt, desc='visco filter', ascii=True)
for n in range(n_start, Nt0):
dt = sqrt(dt_fac*n)/Ts
dt_int = iceil(dt)
nv = np.arange(n-dt_int,n+dt_int+1)
assert n>=dt_int
y[:,nv] += (Ts/sqrt(n*Tsg2pi))*x[:,n][:,None] * exp(-((n-nv)*Ts)**2/(n*Tsg2))[None,:]
nv = np.arange(n-dt_int, n+dt_int+1)
assert n >= dt_int
y[:, nv] += (Ts/sqrt(n*Tsg2pi))*x[:, n][:, None] * \
exp(-((n-nv)*Ts)**2/(n*Tsg2))[None, :]
pbar.update(1)
pbar.close()

return np.squeeze(y) #squeeze to 1d in case
return np.squeeze(y) # squeeze to 1d in case
6 changes: 1 addition & 5 deletions src/python/pffdtd/analysis/response.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import click
import scipy.io.wavfile as wavfile
from scipy.io import wavfile
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import ScalarFormatter
Expand Down Expand Up @@ -48,9 +48,6 @@ def main(filename, fmin, fmax, label_a, label_b, smoothing, target):
file_a = filename[0]
file_b = filename[1]

label_a = label_a
label_b = label_b

fs_a, buf_a = wavfile.read(file_a)
fs_b, buf_b = wavfile.read(file_b)

Expand All @@ -75,7 +72,6 @@ def main(filename, fmin, fmax, label_a, label_b, smoothing, target):
dB_b += 75.0

if smoothing > 0.0:
smoothing = smoothing
dB_a = fractional_octave_smoothing(dB_a, fs_a, nfft, smoothing)
dB_b = fractional_octave_smoothing(dB_b, fs_b, nfft, smoothing)

Expand Down
6 changes: 3 additions & 3 deletions src/python/pffdtd/analysis/room_modes.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import find_peaks
import scipy.io.wavfile as wavfile
from scipy.io import wavfile

from pffdtd.common.myfuncs import iceil
from pffdtd.common.plot import plot_styles
Expand All @@ -18,8 +18,8 @@ def find_nearest(array, value):
return array[idx]


def collect_wav_paths(dir, pattern="*.wav"):
return list(sorted(glob.glob(os.path.join(dir, pattern))))
def collect_wav_paths(folder, pattern="*.wav"):
return list(sorted(glob.glob(os.path.join(folder, pattern))))


def hz_to_note(frequency):
Expand Down
Loading

0 comments on commit f33caca

Please sign in to comment.