Эксперимент для сотрудника с нарушением слуха, ч. 2, проверка на себе

Привет! В первой части мы рассказали, зачем вообще решили заняться этим вопросом, а также поделились переводом статьи, ставшей для нас отправной точкой для собственных изысканий. Теперь хотим рассказать, как мы доработали идею под нашего сотрудника.

Отдельное спасибо комментаторам, которые отметились в комментариях к первой части. Устройства с костной проводимостью, программные решения вроде Equalizer APO 1.2.1, слуховые устройства с поддержкой Bluetooth — мы собрали и передали все ваши идеи. Может быть, что-то из этого и выйдет. Но мы расскажем о своём варианте. Возможно, он тоже кому-то будет полезен.

Подготовка

Для проверки нужно построить свои HL-аудиграммы (‘Hearing Levels’). Сгенерируем wav-образцы на ключевых частотах с постепенно возрастающей громкостью с помощью gen_hl_samples.py, прослушаем их и запишем время, когда стало слышно каждый из образцов. Далее по «линейке» сопоставим время с уровнем громкости и получим свою аудиограмму в виде набора «частота» — «уровень слышимости частоты».

В авторский код gen_hl_samples.py были внесены некоторые изменения:

во-первых, в дополнение к wav-файлам вместо вывода в консоль параметров «время — дБ» стали выводить своеобразную «линейку»: текстовый файлик со шкалой, который помогает по моменту времени определить уровень слышимости (одинаковый для всех ключевых частот)

Пример линейки

MIN      0:00        0:01        0:02        0:03        0:04        0:05        0:06        0:07        0:08        0:09        0:10     MAX ──────────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────    volume │ -90.00 dB │ -80.00 dB │ -70.00 dB │ -60.00 dB │ -50.00 dB │ -40.00 dB │ -30.00 dB │ -20.00 dB │ -10.00 dB │   0.00 dB │    volume amplitude │   0.00003 │   0.00010 │   0.00032 │   0.00100 │   0.00316 │   0.01000 │   0.03162 │   0.10000 │   0.31623 │   1.00000 │ amplitude

во-вторых, сменили направление изменения громкости на возрастание, дабы сразу не оглохнуть, начав прослушивать образец. В оригинале громкость образцов убывала.

gen_hl_samples.py

Генератор образцов для построения своей HL-аудиограммы: генерирует для каждой частоты из списка [125, 250, 500, 1000, 2000, 4000, 6000] wav-файл, громкость нарастает от vol_min до vol_max со скоростью vol_delta дБ/сек

import numpy as np from scipy.io import wavfile   def db_to_amplitude(level_dB):  # преобразование дБ в амлитуду     return np.power(10.0, level_dB / 20.0)   def amplitude_to_db(amplitude):  # преобразование амлитуты в дБ     return 20 * np.log10(amplitude)   def line_aka(time, frequency, volume, time_str='', scale_str='', vol_str='', amp_str=''):     if len(time_str) == 0 or len(scale_str) == 0 or len(vol_str) == 0 or len(amp_str) == 0:         time_str = f"MIN      {time // 60:0>1.0f}:{time % 60:0>2.0f}"  # временные отметки для шкалы         scale_str = f"──────────┼"  # шкала         # vol_str = f"  {frequency: >4} Hz │"  # отметки уровня тона         vol_str = f"   volume │"  # отметки уровня тона         amp_str = f"amplitude │"     else:         time_str += f"{' ' * 8}{time // 60:0>1.0f}:{time % 60:0>2.0f}"         scale_str += f"{'─' * 11}┼"         vol_str += f"{volume: >7.2f} dB │"         amp_str += f" {db_to_amplitude(volume): >9.5f} │"      return time_str, scale_str, vol_str, amp_str   if __name__ == "__main__":     sampleRate = 44100     duration = 1  # сек, длительность звучания тона в образце на текущем уровнем громкости      sampleCount = sampleRate * duration      s = ''  # шкала время/уровень громкости тона     for freq in [125, 250, 500, 1000, 2000, 4000, 6000]:  # Гц, частота тонов образцов         samples = np.array([])          # amp = amp_delta = 0.01  # начальная амлитуда и шаг изменения амплитуды         # amp_max = 0.1           # конечная амплитуда         vol_min = -90           # дБ, начальный уровень громкости         vol_max = 0             # дБ, конечный уровень громкости         vol_delta = 1           # дБ, шаг         t = 0                   # сек, начальный момент времени          s_time, s_scale, s_dB, s_amp = line_aka(t, freq, None)          for vol_dB in range(vol_min, vol_max + vol_delta, vol_delta):             time_points_array = np.linspace(t, t + duration, sampleCount)              samples_new = np.sin(2 * np.pi * freq * time_points_array)             samples_new *= db_to_amplitude(vol_dB)             samples = np.append(samples, samples_new)              t += duration             s_time, s_scale, s_dB, s_amp = line_aka(t, freq, vol_dB, s_time, s_scale, s_dB, s_amp)          # сохранение частотного образца в файл         wavfile.write(f"./wav_files/sine_{freq}.wav", sampleRate, samples)         print(f"{freq: >4} Hz wav file done!")      s = f"{s_time}     MAX\n{s_scale}──────────\n{s_dB}    volume\n{s_amp} amplitude"     # print(s)      with open(f"./wav_files/sine_hearing_levels_scale.txt", 'w', encoding='utf-8') as file_handle:         file_handle.write(s)

 125 Hz wav file done!  250 Hz wav file done!  500 Hz wav file done! 1000 Hz wav file done! 2000 Hz wav file done! 4000 Hz wav file done! 6000 Hz wav file done!

Аудиограммы готовы. Визуализируем данные с помощью matplotlib и numpy.

image_builder.py

# https://blog.demofox.org/2015/04/14/decibels-db-and-amplitude/ #   0 дБ означает полный уровень сигнала: амплитуда == 1.0, т.е. 100% громкость #   каждые 6 дБ сигнал изменяется в ~2 раза #       +n дБ - увеличение громкости #       -n дЬ - уменьшение громкости  import matplotlib.pyplot as plt from matplotlib.ticker import EngFormatter, FuncFormatter, PercentFormatter import numpy as np import matplotlib as mpl from mpl_toolkits.axes_grid1 import make_axes_locatable   def db_to_amplitude(level_db):  # преобразование дБ в амлитуду     return np.power(10.0, level_db / 20.0)   def amplitude_to_db(amplitude):  # преобразование амлитуты в дБ     return 20 * np.log10(amplitude)   audiograms = {  # decibel-frequency audiograms     1: np.array([  # An_         [125, -48.0],         [250, -42.0],         [500, -52.0],         [1000, -36.0],         [2000, -23.0],         [4000, -29.0],         [6000, -19.0],     ]),     2: np.array([  # Mi_         [125, -76.0],         [250, -68.0],         [500, -72.0],         [1000, -73.0],         [2000, -77.0],         [4000, -83.0],         [6000, -78.0],     ])}  for k, a in audiograms.items():     min_y, max_y = a[:, 1].min() * 1.1, 0     img_width, img_height = 16, 9      fig, ax = plt.subplots(1, 1, figsize=(img_width, img_height))     plt.xticks(rotation='vertical')     fig.suptitle(f"Аудиограмма {k} сотрудника Cloud4y")  # , fontsize=14)     fig.patch.set_facecolor('white')      ax.plot(a[:, 0], a[:, 1], '-o')     ax.set_xscale('log')      ax.set_xlabel(r'Частота сигнала')     ax.set_ylabel(r'лучше   <<<   Уровень слуха (слышимость сигнала)   >>>   хуже')     ax.set_ylim(min_y, max_y)     ax.xaxis.set_major_formatter(EngFormatter(unit='Гц'))     ax.yaxis.set_major_formatter(EngFormatter(unit='дБ'))     ax.xaxis.set_ticks(a[:, 0])     ax.yaxis.set_major_locator(plt.MultipleLocator(10))     ax.yaxis.set_minor_locator(plt.MultipleLocator(2))     ax.grid(axis='y')      # дополнительные вертикальные оси - способ 1     y2 = ax.secondary_yaxis('right', functions=(db_to_amplitude, amplitude_to_db))     y2.set_yscale('log')     y2.set_ylabel(r'Амплитуда сигнала')     y2.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: f"{x: .3f}"))      y3 = ax.secondary_yaxis(1.1, functions=(db_to_amplitude, amplitude_to_db))     y3.set_yscale('log')     y3.set_ylabel(r'Громкость сигнала')     # y3.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: f"{x * 100:5.1f} %"))     y3.yaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=1, symbol=' %', is_latex=False))      # # дополнительные вертикальные оси - способ 2 - происходит наложение 2го графика на 1й вместо добавления Y-оси     # #     чтобы графики совпадали, нужно подправлять пределы     # ax2 = ax.twinx()     # ax2.set_yscale('log')     # ax2.plot(df[:, 0], db_to_amplitude(df[:, 1]), 'y:x')     # ax2.set_ylabel(r'Signal Amplitude')     # ax2.set_ylim(db_to_amplitude(min_y), db_to_amplitude(max_y))     # ax2.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: f"{x:.3f}"))     # ax2.xaxis.set_ticks(df[:, 0])      ax.axhline(y=a[:, 1].max() + 1, linestyle="-", color='C2')   # demo_1_audible_normally     ax.axhline(y=np.median(a[:, 1]), linestyle="-", color='C1')  # demo_2_audible_partially     ax.axhline(y=a[:, 1].min() - 1, linestyle="-", color='C3')   # demo_3_not_audible      labels = [         "HL-Аудиограмма (Hearing Level)",         f"demo_1: {a[:, 1].max() + 1: .0f} дБ ≡ {db_to_amplitude(a[:, 1].max() + 1): .5f}, корректировка не потребуется",         f"demo_2: {np.median(a[:, 1]): .0f} дБ ≡ {db_to_amplitude(np.median(a[:, 1])): .5f}, для частичной корректировки",         f"demo_3: {a[:, 1].min() - 1: .0f} дБ ≡ {db_to_amplitude(a[:, 1].min() - 1): .5f}, для полной корректировки по аудиограмме",     ]     plt.legend(labels=labels)      divider = make_axes_locatable(ax)     cax = divider.append_axes("left", size="0.7%", pad=-.09)     cax.set_ylim(min_y, max_y)     norm = mpl.colors.Normalize(vmin=min_y, vmax=max_y)     cmap = mpl.cm.ScalarMappable(norm=norm, cmap='RdYlGn_r', )     fig.colorbar(cmap, cax=cax)     cax.yaxis.set_major_locator(plt.MultipleLocator(10))     cax.yaxis.set_minor_locator(plt.MultipleLocator(2))     cax.set_yticklabels([])      fig.tight_layout()     plt.show()

На каждой аудиограмме отмечены три уровня громкости, ‘нормально слышно’, ‘частично слышно’, ‘совсем не слышно’, на которых мы сгенерируем с помощью gen_continuous_sample.py демки, которые затем обработаем их алгоритмом, чтобы понять, насколько хорошо проходит обработка.

gen_continuous_sample.py

Генерация образцов с тремя разными постоянными уровнями громкости, частота меняется от freq_start до freq_end со скоростью freq_inc/time_inc Гц/сек.

from scipy.io import wavfile import numpy as np   def db_to_amplitude(level_db):  # преобразование дБ в амлитуду     return np.power(10.0, level_db / 20.0)   def amplitude_to_db(amplitude):  # преобразование амлитуты в дБ     return 20 * np.log10(amplitude)   def generate_samples(sample_rate, level=-20.0):     amp = db_to_amplitude(level)  # коэффициент амлитуды (уровня громкости) сигнала (1 ≡ 0 dB)      freq_start = 125  # Гц     freq_end = 8000  # Гц      freq_inc = 25  # Гц, шаг увеличения тональности     time_inc = .25  # сек, длительность сигнала перед увеличением тональности      print(f"amplitude {amp: .5f} ≡ {np.round(20 * np.log10(amp), 2)} dB")      inc_sample_count = int(sample_rate * time_inc)     sample = np.array([])      freq = freq_start     t = 0.0     while freq < freq_end:         time_points_array = np.linspace(t, t + time_inc, inc_sample_count, endpoint=False)          new_samples = np.sin(2 * np.pi * freq * time_points_array)         new_samples *= amp          sample = np.append(sample, new_samples)          freq += freq_inc         t += time_inc      return sample   if __name__ == '__main__':     levels_a = {'demo_1_audible_normally': -18.0, 'demo_2_audible_partially': -36.0, 'demo_3_not_audible': -53.0}  # An_     levels_m = {'demo_1_audible_normally': -67.0, 'demo_2_audible_partially': -76.0, 'demo_3_not_audible': -84.0}  # Mi_     levels = {1: levels_a, 2: levels_m}     sampleRate = 44100      for k, l in levels.items():         for name, level in l.items():             generated_samples = generate_samples(sampleRate, level=level)             wavfile.write(f"./wav_files/{k}_{name}.wav", sampleRate, generated_samples)

amplitude  0.12589 ≡ -18.0 dB amplitude  0.01585 ≡ -36.0 dB amplitude  0.00224 ≡ -53.0 dB amplitude  0.00045 ≡ -67.0 dB amplitude  0.00016 ≡ -76.0 dB amplitude  0.00006 ≡ -84.0 dB

Демо-wav готовы для обработки с помощью gain.py.

gain.py

from scipy.io import wavfile  import numpy as np import time import os.path   # находит и возвращает главную частоту в заданном окне и уровень сигнала def get_dominant_freq(sample_rate, window):     # область частот окна     yf = np.fft.fft(window)     yf = np.abs(2 * yf / len(window))      window_size = len(window)     window_half = len(window) // 2      max_amp = yf[:window_half].max()     max_amp_idx = yf[:window_half].argmax()      # поиск частоты по её индексу     frequency = sample_rate * max_amp_idx / window_size  # frequency = (sample_rate/2) * max_amp_idx / (window_size/2)      return frequency, amplitude_to_db(max_amp)   def get_gain(freq, dB, audiogram):  # возвращает необходимый сдвиг уровня для сигнала для указанной частоты согласно аудиограмме     threshold = None      if freq < audiogram[0][0]:         threshold = audiogram[0][1]     else:         for i in range(len(audiogram) - 1):             if freq >= audiogram[i + 1][0]:                 continue             threshold = (audiogram[i + 1][1] - audiogram[i][1]) * \                         (freq - audiogram[i][0]) / (audiogram[i + 1][0] - audiogram[i][0]) + audiogram[i][1]             break          if threshold is None:             threshold = audiogram[-1][1]      if threshold <= dB:  # уровень сигнала уже в зоне слышимости         return 0.0     else:  # сдвиг уровня сигнала в зону слышимости         return threshold - dB + 3   def db_to_amplitude(dB):  # преобразование дБ в амлитуду     return np.power(10.0, dB / 20.0)   def amplitude_to_db(amp):  # преобразование амлитуты в дБ     return 20 * np.log10(amp)   if __name__ == "__main__":     start_time = time.time()      df_1 = np.array([  # An_         [125, -48.0],         [250, -42.0],         [500, -52.0],         [1000, -36.0],         [2000, -23.0],         [4000, -29.0],         [6000, -19.0],     ])     df_2 = np.array([  # Mi_         [125, -76.0],         [250, -68.0],         [500, -72.0],         [1000, -73.0],         [2000, -77.0],         [4000, -83.0],         [6000, -78.0],     ])      # считывание демо-файлов     wavfiles = [         ['1_demo_1_audible_normally.wav', '1_demo_2_audible_partially.wav', '1_demo_3_not_audible.wav'],         ['2_demo_1_audible_normally.wav', '2_demo_2_audible_partially.wav', '2_demo_3_not_audible.wav']     ]      for files, audiogram in zip(wavfiles, [df_1, df_2]):         for file in files:             file_name = file.split('.')[0]             sampleRate, samples = wavfile.read(os.path.join(os.getcwd(), 'wav_files', file), mmap=True)             outSamples = np.array([])  # набор обработанных участков для нового файла             windowSize = 512             sampleIdx = 0             while sampleIdx < len(samples):                 window = samples[sampleIdx:sampleIdx + windowSize]                  freq, dB = get_dominant_freq(sampleRate, window)                 gain = get_gain(freq, dB, audiogram=audiogram)                  startSec = np.round(sampleIdx / sampleRate, 2)                 endSec = np.round((sampleIdx + windowSize) / sampleRate, 2)                  if gain > 0.0:                     # print(f"{file_name}\t{startSec: >6.3f}..{endSec: <6.3f}s\tGAINED\t{freq: >8.2f} Hz\tvol {dB: >6.2f} dB\tgain: {gain: >6.2f} dB")                     amp = db_to_amplitude(gain)                     amplified = window * amp                 else:                     # print(f"{file_name}\t{startSec: >6.3f}..{endSec: <6.3f}s\tSTOCK\t{freq: >8.2f} Hz\tvol {dB: >6.2f} dB\tgain: {gain: >6.2f} dB")                     amplified = window                  outSamples = np.append(outSamples, amplified)                  sampleIdx += windowSize              processed_file = os.path.join(os.getcwd(), 'wav_files', f"{file_name}_processed.wav")             wavfile.write(processed_file, sampleRate, outSamples)             print('ready', os.path.relpath(processed_file))      print(f"{(time.time() - start_time)} seconds")

ready wav_files\1_demo_1_audible_normally_processed.wav ready wav_files\1_demo_2_audible_partially_processed.wav ready wav_files\1_demo_3_not_audible_processed.wav ready wav_files\2_demo_1_audible_normally_processed.wav ready wav_files\2_demo_2_audible_partially_processed.wav ready wav_files\2_demo_3_not_audible_processed.wav 180.6359121799469 seconds

Всё готово! Можно послушать демки, оценить обработку, попробовать обработать реальные аудиофайлы, наконец. Результат не идеален, но вполне рабочий, демки стало слышно.

Ради интереса и объективного сравнения можно даже посмотреть файлы «до» и «после» в Praat. Становится видно изъян алгоритма: усиливаются те частоты, которые заведомо находятся в зоне слышимости и усиливаться не должны. Это особенно заметно на первой демке, которая была сгенерирована заведомо слышимой на всем частотном диапазоне: