Skip to content

Commit

Permalink
[python] Improve RT60 tolerance report
Browse files Browse the repository at this point in the history
  • Loading branch information
@tobiashienzsch
tobiashienzsch committed Aug 7, 2024
1 parent c25eef1 commit 2ced3dd
Show file tree
Hide file tree
Showing 4 changed files with 100 additions and 77 deletions.
10 changes: 6 additions & 4 deletions run_3d.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,8 @@ sim_dir="$root_dir/data/sim_data/$sim_name/cpu"
model_dir="$root_dir/data/models/$sim_name"
materials_dir="$root_dir/data/materials"

fmax=1000
fmin=20
fmax=400

# Delete old sim
rm -rf "$sim_dir"
Expand All @@ -34,6 +35,7 @@ $engine_exe

# Post-process
cd "$python_dir"
python -m sim3d.process_outputs --data_dir="$sim_dir" --fcut_lowpass "$fmax" --N_order_lowpass=8 --symmetric --fcut_lowcut 20.0 --N_order_lowcut=4 --air_abs_filter="none" --save_wav --plot
python -m analysis.t60 --data_dir="$sim_dir" --fmin=20 --fmax="$fmax"
python -m analysis.room_modes --data_dir="$sim_dir" --fmin=10 --fmax=200 --modes=10
python -m sim3d.process_outputs --data_dir="$sim_dir" --fcut_lowpass "$fmax" --N_order_lowpass=8 --symmetric --fcut_lowcut $fmin --N_order_lowcut=4 --air_abs_filter="none" --save_wav --plot
# python -m analysis.t60 --fmin=$fmin --fmax="$fmax" --target=0.25 ../../data/sim_data/$sim_name/cpu/R001_out_normalised.wav
# python -m analysis.t60 --data_dir="$sim_dir" --fmin=$fmin --fmax="$fmax" --target=0.25
python -m analysis.room_modes --data_dir="$sim_dir" --fmin=$fmin --fmax=$fmax --modes=20
2 changes: 1 addition & 1 deletion src/python/Studio_cpu.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
'Table': 'door_wood.h5',
'Walls': 'concrete_painted.h5',
},
duration=1.5,
duration=2.0,
Tc=20,
rh=50,
fcc_flag=False,
Expand Down
10 changes: 6 additions & 4 deletions src/python/analysis/room_modes.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import glob
import os
import argparse
import pathlib

import matplotlib.pyplot as plt
import numpy as np
Expand All @@ -14,7 +15,7 @@
def collect_wav_files(directory, pattern="*.wav"):
search_pattern = os.path.join(directory, pattern)
wav_files = glob.glob(search_pattern)
return wav_files
return list(sorted(wav_files))


def hz_to_note(frequency):
Expand Down Expand Up @@ -141,6 +142,7 @@ def main():

mode_freqs = [mode['frequency'] for mode in modes][:args.modes]
for file in files:
file = pathlib.Path(file)
fs, buf = wavfile.read(file)
fmin = args.fmin if args.fmin > 0 else 10
fmax = args.fmax if args.fmax > 0 else fs/2
Expand All @@ -159,13 +161,13 @@ def main():

print(freqs[peaks][:10])

plt.plot(freqs, dB, linestyle='-', label=f'Receiver')
plt.plot(freqs, dB, linestyle='-', label="Measurement")
if args.modes > 0:
plt.vlines(mode_freqs, dB_max-80, dB_max+10,
colors='#AAAAAA', linestyles='--')
colors='#AAAAAA', linestyles='--', label="Modes")
plt.plot(freqs[peaks], dB[peaks], 'r.',
markersize=10, label='Peaks')
plt.title('Response')
plt.title(f'{file.stem[:4]}')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.xscale('log')
Expand Down
155 changes: 87 additions & 68 deletions src/python/analysis/t60.py
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@
import argparse
import glob
import os
import pathlib

import numpy as np
import scipy.io.wavfile as wavfile
import scipy.signal as signal
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from matplotlib.ticker import ScalarFormatter

from common.plot import plot_styles
Expand Down Expand Up @@ -40,30 +42,7 @@ def calculate_t60(edc_db, fs):
return t60


def ebu_3000_t60_threshold_upper(freqs):
times = np.zeros_like(freqs)
for i, freq in enumerate(freqs):
if freq < 63:
times[i] = 0.3
elif freq <= 200:
times[i] = 0.3 - (0.25 * (freq - 63) / (200 - 63))
else:
times[i] = 0.05
return times


def ebu_3000_t60_threshold_lower(freqs):
times = np.zeros_like(freqs)
for i, freq in enumerate(freqs):
if freq < 4000:
times[i] = -0.05
else:
times[i] = -0.10
return times


def run(files, fmin, fmax, show_all=False, show_tolerance=True):

def run(files, fmin, fmax, show_all=False, show_tolerance=True, target=None):
# ISO 1/3 octaves
center_freqs = np.array([
20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160,
Expand All @@ -76,10 +55,12 @@ def run(files, fmin, fmax, show_all=False, show_tolerance=True):
center_freqs = center_freqs[np.where(center_freqs <= fmax)]

file_times = []
for file in files:
file_names = []
for path in files:
file = pathlib.Path(path).absolute()
fs, ir = wavfile.read(file)
t60_times = []
print(f"{file}")
print(f"---- {file.stem} ----")
for center_freq in center_freqs:
filtered_ir = third_octave_filter(ir, fs, center_freq)
edc_db = energy_decay_curve(filtered_ir)
Expand All @@ -88,56 +69,88 @@ def run(files, fmin, fmax, show_all=False, show_tolerance=True):
print(f"T60 at {center_freq} Hz: {t60:.2f} seconds")

file_times.append(np.array(t60_times))
file_names.append(file.stem[:4])

plt.rcParams.update(plot_styles)

fig, ax = plt.subplots(2, 1)
fig, axs = plt.subplots(2, 1)
formatter = ScalarFormatter()
formatter.set_scientific(False)

# T60
ax[0].margins(0, 0.1)
ax: Axes = axs[0]
ax.margins(0, 0.1)

if show_all:
for i, f in enumerate(file_times):
ax[0].semilogx(center_freqs, f, label=f"{i}")
for f, name in zip(file_times, file_names):
ax.semilogx(center_freqs, f, label=f"{name}")
else:
ax[0].semilogx(center_freqs, file_times[0])

ax[0].set_title("Reverberation Time (T60)")
ax[0].set_ylabel("Decay [s]")
ax[0].set_xlabel("Frequency [Hz]")
ax[0].xaxis.set_major_formatter(formatter)

ax[0].set_xlim((center_freqs[0], center_freqs[-1]))
ax[0].set_ylim((0.0, np.max(file_times)+0.1))

ax[0].grid(which='minor', color='#DDDDDD', linestyle=':', linewidth=0.5)
ax[0].minorticks_on()
ax[0].legend()
ax.semilogx(center_freqs, file_times[0], label="Measurement")

if target:
ax.hlines(
target,
fmin,
fmax,
color='#555555',
label=f"Target {target} s",
linestyles="dashed",
)

ax.set_title("RT60")
ax.set_ylabel("Decay [s]")
ax.set_xlabel("Frequency [Hz]")
ax.xaxis.set_major_formatter(formatter)

ax.set_xlim((fmin, fmax))
ax.set_ylim((0.0, np.max(file_times)+0.1))

ax.grid(which='minor', color='#DDDDDD', linestyle=':', linewidth=0.5)
ax.minorticks_on()
ax.legend(loc='upper right')

# Tolerance
if show_tolerance:
diff = np.insert(np.diff(file_times[0]), 0, 0.0)
diff = np.insert(file_times[0][:-1]-file_times[0][1:], 0, 0.0)
upper_threshold = ebu_3000_t60_threshold_upper(center_freqs)
lower_threshold = ebu_3000_t60_threshold_lower(center_freqs)

ax[1].margins(0, 0.1)
ax[1].semilogx(center_freqs, diff)
ax[1].semilogx(center_freqs, upper_threshold)
ax[1].semilogx(center_freqs, lower_threshold)

ax[1].set_title(f"T60 Tolerance (EBU Tech 3000)")
ax[1].set_ylabel("Difference [s]")
ax[1].set_xlabel("Frequency [Hz]")
ax[1].xaxis.set_major_formatter(formatter)

ax[1].set_xlim((center_freqs[0], center_freqs[-1]))
ax[1].set_ylim((np.min(diff)-0.1, 0.4))

ax[1].grid(which='minor', color='#DDDDDD',
linestyle=':', linewidth=0.5)
ax[1].minorticks_on()
k4 = min(fmax, 4000)
k8 = min(fmax, 8000)
k20 = min(fmax, 20000)

ax: Axes = axs[1]
ax.margins(0, 0.1)
ymin, ymax = -0.05, +0.3
if show_all:
for f, name in zip(file_times, file_names):
diff = np.insert(np.diff(f), 0, 0.0)
diff = np.insert(f[:-1]-f[1:], 0, 0.0)
ax.semilogx(center_freqs, diff, label=f"{name}")
ymin, ymax = min(ymin, np.min(diff)), max(ymax, np.max(diff))
else:
diff = np.insert(np.diff(file_times[0]), 0, 0.0)
diff = np.insert(file_times[0][:-1]-file_times[0][1:], 0, 0.0)
ax.semilogx(center_freqs, diff, label="Measurement")
ymin, ymax = np.min(diff), np.max(diff)

ax.plot([63.0, 200.0], [0.3, 0.05], color='#555555')
ax.hlines(
[+0.05, -0.05, -0.1, -0.1, +0.3, +0.05, -0.05],
[200, 100, k4, k8, fmin, k8, fmin],
[k8, k4, k8, k20, 63, k20, 100],
linestyles=["-", "-", "-", "--", "--", "--", "--"],
colors='#555555',
label="EBU Tech 3000"
)

ax.set_title(f"Tolerance")
ax.set_ylabel("Difference [s]")
ax.set_xlabel("Frequency [Hz]")
ax.xaxis.set_major_formatter(formatter)

ax.set_xlim((fmin, fmax))
ax.set_ylim((ymin-0.05, ymax+0.05))

ax.grid(which='minor', color='#DDDDDD', linestyle=':', linewidth=0.5)
ax.minorticks_on()
ax.legend(loc='upper right')

plt.show()

Expand All @@ -148,18 +161,24 @@ def main():
parser.add_argument('--data_dir', type=str)
parser.add_argument('--fmin', type=float, default=20.0)
parser.add_argument('--fmax', type=float, default=1000.0)
parser.add_argument('--target', type=float)

args = parser.parse_args()
if args.data_dir and len(args.filename) > 0:
raise RuntimeError("--data_dir not valid, when comparing IRs")

files = args.filename
if args.data_dir:
files = collect_wav_files(args.data_dir, "*_out_normalised.wav")
run(files, args.fmin, args.fmax, show_tolerance=True)

if len(args.filename) == 2:
files = [args.filename[0], args.filename[1]]
run(files, args.fmin, args.fmax, show_all=True, show_tolerance=False)
run(
list(sorted(files)),
args.fmin,
args.fmax,
show_tolerance=True,
show_all=True,
target=args.target
)


if __name__ == '__main__':
Expand Down

0 comments on commit 2ced3dd

Please sign in to comment.