Резонансные датчики в оцифрованном звуке

  Простенький физический аналог - грузик на пружинке, позволяет исследовать спектральные особенности оцифрованного звука не хуже, чем преобразования Фурье. Датчик, работающий на принципе резонанса, позволяет в общих чертах, но сравнительно просто, извлекать нотную запись из музыкального файла - http://www.proza.ru/2016/08/18/1254
  У таких датчиков, установленных в программу, по которой идёт перекачивание звука из одного звукового файла в другой, может быть масса полезных применений. И наверняка, нечто подобное применяется, но ничего конкретного по таким датчикам в интернете я найти не смог. А потому провёл собственное исследование, о результатах которого здесь рассказываю.


   КАК УСТРОЕН РЕЗОНАНСНЫЙ ДАТЧИК

  Представьте программу, которая читает музыкальный WAV файл и анализирует его. Из файла (неважно, стерео он или моно) программа получает амплитуды звука - целые числа в двухбайтовом преставлении в интервале от -32768 до +32767. В своих прикидках на возможную величину амплитуды будем ориентироваться на не очень громкий звук, полагая A=4000.
  Сопоставим это число с некой силой, действующей на груз массы m на пружинке жёсткости K. Второй закон Ньютона для этой системы запишем в виде:

  A/m - (K/m)*X = dV/dT  (1)

где dT разумеется равно 1/44100 сек - возьмём для определённости именно такую скорость раздачи, 44100 семплов в секунду.
  Резонансная частота этой системы (не в конечно разностном, а в нормальном её представлении) связана с величинами k и m следующим образом:

  K/m = (Fрез*2*Пи)*(Fрез*2*Пи)    равно 4042590 ньютон/м*кг, для Fрез=320 Гц

  Частота 320 Гц находится в первой октаве, и мы будем ориентироваться пока на такие звуки, хотя в экспериментах уйдём и в килогерцы. И пусть X в формуле (1) будет что-то в районе единицы, а раскачивающее слагаемое в левой части пусть будет на порядок меньше возвращающего фактора. Из этих соображений получаем:

  4000/m = 1/10 * 4042590 * 1    или, примерно, 1/m = 100

  Теперь, когда мы определились со всеми коэффициентами. Мы можем (обозначив через Km отношение этих величин K/m) написать кусочек программы, представляющий собой конечно разностное воплощение резонансного датчика:

   dV=(A*100-Km*X)*dT : V=V+dV 'датчик частоты          (2)
   V=V*.995 : X=X+V*dT : X=X*.995               
   Z=Z+ABS(X) : Z=Z*.9 'амплитуда датчика
   
  Поставив этот датчик в поток идущих амплитуд, по амплитуде его колебаний Z мы можем судить о том, имеются ли в звуке данная частота или нет.

  Во фрагменте (2), в дополнение к формуле (1), введён диссипативный фактор, не позволяющий датчику раскачиваться очень сильно в резонансе, и, одновременно, расширяющий полосу его чувствительности. Диссипативный фактор поровну распределён между V и X. Диссипативный фактор имеется и у интегрированной величины Z, создавая усреднение приблизительно 10-ти подряд идущих амплитуд колебаний датчика.
  Фрагмент (2) по написанию и, возможно, по скорости выполнения можно сократить:

   V=((A*100-Km*X)*dT+V)*.995 'датчик частоты            (3)
   X=(V*dT+X)*.995               
   Z=(ABS(X)+Z)*.9 'амплитуда датчика

  Если датчиков несколько, то можно заранее провести умножение A*100


   ИССЕДОВАНИЕ ПОВЕДЕНИЯ РЕЗОНАНСНОГО ДАТЧИКА

  Как выяснилось, резонансные частоты датчика (2) отличаются (особенно на больших частотах) от теоретических резонансных частот формулы (1), которые верны только в случае  dT = 1/44100 << 1/Fрез
  Поведение датчика в оцифрованном звуке изучалось с помощью программы, представленной в Приложении 1.
  Датчик настраивался на теоретическую частоту F(1). Программой создавалась синусоида оцифрованного звука частоты F и находилась резонансная частота F(2), соответствующая максимальным колебаниям датчика (в таблице указаны значения Xmax*100). Также определялась полуширина чувствительности (графики представлены на рисунке) по уровням 0.71 0.5 0.35 0.25 0.18 0.125 0.09 Эта величина совершенно не зависела от частот, вплоть до критических, 13900-14073 Гц по шкале F(1).
  Значение ZM оказывалось пропорциональным Xmax, но, естественно, оно зависит от диссипативных факторов.
  Диссипативный фактор совершенно естественным и предсказуемым образом влияет и на ширину полосу чувствительности датчика и на его чувствительность в максимуме. Произведение ширины полосы на чувствительность в максимуме для каждой частоты F(2) является постоянной величиной.

Результаты  100*   полуширина полосы чувст.
  F(1) F(2) Xmax    по разным уровням, Гц        *100
  500   500 28.73  \45 85 135 215 275 405     ZM=224
 1000  1000 14.33  \45 85 145 225 340 475 590
 2000  2000  7.234 \45 85 145 225 335 495
 4000  4045  3.736 \45 85 145 230 340 505 715
 8000  8490  2.175 \40 85 140 220 325 500 615 ZM=13
 9000  9740  2.075

10000 11095  2.036 \47 85 145 215 320 495
10250 11455  2.045
10500 11825  2.054
11000 12600  2.085 \45 90 145 230 345
12000 14335  2.285 \45 85 145 230 285 485 575
13000 16535  2.877 \45 85 145 230 245 505 745

13500 18035  3.756 \40 85 145 222 340 505 735
13775 19040  4.892 \45 85 145 225 340 510
13900 19850  6.594 \45 85 145 230 345 505
14000 20625 10.07  \40 85 145 230 340         ZM=66.9
14050 21255 17.91  \40 85 145 230 345 535 750 795
14070 21755 50.81  \45 95 160 255 275 и пошло вверх
14073 22015 286715187200 \25 30 35 и пошло вверх

14073 22085 287517849600 \50 95 155 240 335
14070 22330 50.93  \45 80 135
14050 22850 17.92  \40 85 140 220 340 495
14000 23480 10.07  \40 85 140 220 330 500
13900 24255  6.596 \ 40 85 140 220 335 500

  Интересно проследить зависимость чувствительности датчика в максимуме Xmax, и наблюдаемой частоты резонанса Fрез=F(2), от предустановленной частоты F(2). Последняя зависимость приводится на рисунке.
 
  /// В аналитическом виде обратная зависимость может быть представлена следующей приближённой формулой (частоты в ней указываются в килогерцах):

 F(1) = 14.1 - (Fрез-22.05)^2/28.171 + (Fрез-22.05)^4/73689    (4)

чтобы воспользоваться этой формулой на другой частоте раздачи, не равной 44100 семпл/сек, а равной, например 32000, нужно на входе в неё выполнить преобразование частоты Fрез=Fрез*44100/32000, а после выхода из неё выполнить FX=FX*32000/44100. \\\

  Обсудим результаты.
  На небольших частотах частота F(2) следует за частотой F(1), а чувствительность датчика оказывается обратно пропорциональной этим частотам, уменьшаясь с ростом частоты.
  Однако когда значение 1/F(1) оказывается сопоставимым со значением 1/44100, тогда уменьшение чувствительности замедляется и эта зависимость выполаживается. Одновременно значение  F(2) отрывается от значения F(1) и начинает возрастать, сначала немного, а потом всё более значительно.

  На фоне сильного роста разницы F(2)-F(1) резко возрастает и чувствительность датчика. При установке частоты выше F(1)=14073 Гц датчик переполняется.
  Критические частоты F(1)>13900 Гц, близкие к частоте переполнения датчика характеризуются наличием двух резонансных частот F(2), расположенных симметрично выше и ниже частоты переполнения.

  Резонанс на этих частотах вполне реален, он наблюдается не только в искусственных условиях возбуждения одной частотой, но и при взаимодействии датчика с реальной звуковой записью, и даже несмотря на то, что частот таких в записи в принципе не должно было бы быть - они отрезаны фильтром.


   ПОВЕДЕНИЕ РЕЗОНАНСНОГО ДАТЧИКА В РЕАЛЬНОЙ ЗВУКОВОЙ ЗАПИСИ      
 
  Слайды рисунка демонстрируют колебания датчиков, настроенных на разные частоты.
  Видно, что колебания датчиков вполне дифференцированы по частотам, да они и звучат по-разному, выделяя различные детали звука, за которыми вместе с тем узнаётся и общая мелодия.
  Спектр того фрагмента, который исследовался, показывает, что запись была подвержена воздействию фильтра, предположительно, при переводе её в mp3 формат. В записи отрезаны частоты выше 19.5 кГц.

  Постановка датчиков в частоты, присутствующие в спектре, никаких неожиданностей не даёт. Но, как уже отмечалось, можно поставить датчик и в незаполненную полосу частот, и там он тоже будет звучать. И это вполне понятно - резонансный датчик не идеальный объект, полоса его чувствительности по низкому уровню имеет широкие крылья. Однако, возбуждаясь по низким частотам, он иногда начинает резонировать и сам - на собственной резонансной частоте. Это хорошо видно по слайдам спектров в верхнем правом углу рисунка. Пики колебаний возникают там "на пустом месте". На частоте 14000 видны даже два пика колебаний. Разумеется, эти колебания не одновременны, они возникли на разных частях фрагмента, но условия для их возникновения датчик предоставил, имея две резонансные частоты, симметрично расположенные относительно частоты 22050 Гц, что тоже хорошо видно по рисунку.

  Колебания поставленных резонансных датчиков вполне можно складывать с приходящим звуком, или вычитать из него, проводя тем самым коррекцию тембра. Датчик в этом качестве будет выполнять роль эквалайзера (см. рисунок).

   ПРОБЛЕМА РЕКОНСТРУКЦИИ УТРАЧЕННЫХ ЧАСТОТ

  Проблема реставрации цифровых записей со звуком, испорченным удалением верхних полос частот существует. Её можно решать по разному, я в этих вопросах несведущ, и обсуждать их не буду. Здесь я только хочу указать на то, что резонансные датчики, поставленные в частотную полосу, в которой звук сохранился, возможно могут в какой-то степени восстановить утраченный звук.
  И в самом деле - высокие частоты обычно насыщены обертонами того, что слышится на низких частотах и создав эти обертоны в виде синусоид, управляемых по амплитуде параметром Z от датчиков, установленных ниже, мы наверное сможем утраченные обертоны воссоздать. и скорректировать спектральное содержание записи. Так, как показано на рисунке спектра красным.
  Я попробовал такое сделать, и 18 датчиков, поставленных в ряд, с этой задачей успешно справились. По крайней мере - на вид. На слух результата, естественно, никакого, поскольку эти частоты не слышны. Но по крайней мере, то, что имелось, не испортил, и это уже хорошо.
  Интересно попробовать на более испорченных записях. Это у меня в планах, но этого я ещё не делал.

  Общий алгоритм преобразования записи представляется таким:

1. определяем полосу частот F3-F4, в которую будем добавлять модулированные синусоиды, и указываем на число датчиков и на другие характеристики преобразования.

2. указываем на верхнюю границу F2 той полосы с сохранившимся звуком, в которой будут расставляются датчики.

3. запускаем программу. Она, следуя указанному числу датчиков, распределяет частоты F34 равномерно по заполняемой полосе, а частоты F(1) для расставляемых датчиков определяет по формуле (4), подставляя в неё в качестве резонансных значения
  Fрез = F34 * F2/F4
после чего, читает запись и преобразует её. Программа добавляет в запись синусоиды вида

  A34 = sin( F34*2*pi*T + FA ) * Z * AM     (5)

   с подходящей амплитудой AM, умножая их на амплитуды колебаний датчиков.

  Насколько хорошо получается, трудно проверить, поскольку высокие частоты слышны плохо. Некоторое представление об этом может дать запись фрагмента речи, находящаяся тут - http://yadi.sk/d/pKkynY6QwkVw2

  В начале идёт источник, образец речи на лёгком музыкальном фоне, а затем - его преобразование через синусы в ту же полосу частот (250-4000 Гц). Понадобилось большое число датчиков (385), чтобы перекрыть эту полосу , кроме того в формулу (5) добавлена случайная фаза FA, чтобы избежать хлопков на регулярных синусоидах.
  Результат, хоть и разборчив, но сильно хрипит и по качеству сильно хуже исходного звука. Возможно, что хрипы связаны как раз с тем, что синусоиды были взяты регулярными при случайных фазах, тогда как в естественной речи наоборот - синусоиды не регулярны, но фазы этих синусоид на отдельных участках звука не случайны, а связаны между собой.
  Как бы то ни было, но то, что речь получилась членораздельной и воспринимаемой на слух, и тихий музыкальный фон слышен, - это уже хорошо.
 ==============================
ПРИЛОЖЕНИЕ 1
 Программа для изучения поведения резонансного датчика в оцифрованном звуке.

pi=3.14159265 : dT=1/44100
OPEN "ACX.txt" FOR OUTPUT AS #1
PRINT "==== ";
 FX=10000 : FR=11000
 Km=FX*2*PI : Km=Km*Km : XMM=0 : XK=0
 FOR F=FR-500 TO FR+500 STEP 5
   XM=0 : X=0 : V=0 : Z=0 : ZM=0
   FOR T=0 TO 0.1 STEP dT
    A=SIN(F*2*pi*T)*4000
   dV=(A*100-Km*X)*dT : V=V+dV : V=V*.995
   X=X+V*dT : X=X*.995 : IF X>XM THEN XM=X
   Z=Z+ABS(X) : Z=Z*.9 : IF Z>ZM THEN ZM=Z
   NEXT T
   IF XM>XMM THEN XMM=XM : XK=XM : FM=F
   IF XM<XK THEN XK=XK/1.4142 : DF=INT(F-FM) : PRINT #1,"=====\";DF
   PRINT #1,F;XM*100;"  ZM=";ZM*100 : NEXT F
 CLOSE #1
 PRINT " Ok" : STOP