О математическом моделировании звука

  Математическое моделирование звука можно считать неким хобби на границе программирования и музыки. И действительно, - зачем стараться, если существуют хорошие имитаторы в виде электромузыкальных клавишных инструментов с возможностями перестройки своего звучания?

  Ответ простой - всё интересно делать самому.
  Но кроме того, что красивый музыкальный звук интересно сделать самому, делая его, интересно также понять, как он устроен. Почему он звучит именно так. По каким алгоритмам можно его создать.
  В интернете нет подробной информации об этом.
  Кроме того в воздухе витает идея "фикс" - сделать идеальный звук красивого тембра. Ведь у реальных инструментов есть недостатки - и клавиши стучат, и смычок скрипит.
  Но вместе с тем, у реальных инструментов есть и колоссальные достоинства - они совершенствовались веками, и вобрали в себя лучшие достижения акустики.
  Математическое моделирование позволяет проводить экспериментальные исследования и по этой проблеме.

  В самом деле, многое ли мы знаем о физике звучания скрипки или альта?
  Нити смычка натирают канифолью для лучшего сцепления со струной, и есть представление о том, что движение смычка отклоняет струну; сцепленная канифолью, она движется в ту же сторону, что и смычок, потом срывается и ... И что - и?
Как она колеблется, что происходит дальше, как смычок вновь ловит струну, как математически имитировать созданный звук? Вспомните, что струны поддерживаются подставкой, а под поставкой и резонатор есть, и не абы какой, а некой, оптимальной для звука формы. Можно ли сделать математическую имитацию всего этого, не влезая в безумные сложности?

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


  ДЕЛАЕМ ПРОСТОЙ ЗВУК В ФОРМЕ СИНУСОИДЫ

  Как звук формы синусоиды создать программным путём? О, скажет читатель, знаю - есть такая функция Sin(t), используйте её, и делов-то!
  Нет, функция Sin вычисляется как полином, мы же, для ускорения вычислений, будем пользоваться только первым квадратичным членом этого полинома.

  Уточним задание - мы хотим создать WAV файл со звуком нескольких нот, каждая из которых должна звучать секунду или полсекунды, и иметь простейшую временную форму своего звучания - Sin(t).

  Мы будем выполнять это задание с помощью функции создания звука FNS1(F,K). К этой функции мы будем обращаться каждый раз пред записью очередного семпла в звуковой файл чтобы получить амплитуду этого семпла. Очень часто будем обращаться, поэтому надо постараться, чтобы функция наша работала быстро. Текст функции приведён в конце этой статьи.
 
  Прежде всего укажем скорость раздачи - VEL=32000 семпл/сек
Столько семплов нам нужно записывать в файл для обеспечения одной секунды звучания.
Теперь можно начинать работу.

  Чтобы сообщить функции об этом нашем желании, мы обращаемся к ней - A=FNS1(0,0)
Получив ноль в качестве первого параметра, функция понимает, что дальше нужно будет работать, и эту работу, постепенно развивая звук, следует начать с начала звука, с так называемой "атаки", если атака предусмотрена сценарием, и прежде чем начать работать, указанный сценарий следует просмотреть.

  Шаг второй - делаем первое рабочее обращение к функции A=FNS1(F,K) - здесь F - частота звука в герцах, K - номер сценария, по которому будет развиваться звук. В тексте, приведённом ниже, для получения синусоиды используется нулевой сценарий со значением K=0.
  При первом рабочем обращении функция совершает массу действий - устанавливает все параметры сценариев в значения, принятые по умолчанию, устанавливает время в ноль - T=0, по скорости раздачи определяет величину шага по времени DT=1/VEL, а по частоте -  длительность периода колебаний TT=1/FS.
  Произведя общие установки, функция  ищет и находит указанный ей сценарий K.
  Найдя его, она проводит действия, указанные в нём. Обычно в нём находятся операторы присвоения, устанавливающие численные значения параметров.
 
  Сценарии определяют всё - обращаясь к одному из них, Вы при дальнейшем развитии звука получаете звук фортепьяно, обращаясь к другому - звук альта, к третьему - звук аккордеона.
  Мы же, обращаясь сейчас A=FNS1(260.4,0) получим звук синусоиды на частоте До первой октавы, причём даже атака в сценарии на предусмотрена - MU1=0, также, как не предусмотрено затухание звука - DU1=1. Только амплитуда синусоиды указана - AU1=1.

  Установив сценарий и двигаясь дальше по программе, мы приходим на приращение времени и делаем шаг - T=T+DT, и вслед за этим шагом производим некие вычисления - X=T/TT : NT=INT(X) : X=(X-NT)*2
Объясню, что они значат. Определяем X как число периодов колебаний прошедших от начала звука. Вообще-то это число не целое. Функция INT(X) находит целую часть - порядковый номер текущего периода (первый период считается нулевым, за ним следуют периоды 1, 2 ... по порядку). Величина Х перевычисляется, теперь она будет указывать на фазу колебания внутри текущего периода - число в интервале от 0 до 2-х (умножать это число на Пи для вычисления синуса не нужно).

  Если Х>1, то вычисляя последовательно V=X*2-1 и A=1-v*v,
  а иначе - V=(2-X)*2-1 и A=-(1-v*v), минус потому, что значения на 2-м полупериоде должны быть отрицательны,
  и потом - A1=AU1*A,
получаем аппроксимированное параболой значение искомой синусоиды для времени T, прошедшего с момента начала звука. 

  Аналогичные значения будет возвращать функция A=FNS1(260.4,0) при каждом обращении к ней в дальнейшем. Только начальную установку вновь делать не будет, и описание сценариев минует.
  Нам останется только умножать получаемое от функции значение A на достаточно большое число, 2000 например, и записывать полученное значение в звуковой файл в качестве очередного семпла.
  Чтобы заполнить семплами файл, звучащий одну секунду, такие действия нужно повторить 32000 раза.
В программе вычисления идут по более сложному пути, но сути вычислений это не меняет.

  Заметьте, для второго полупериода был сделан обратный ход по времени, взято (2-X)*2, а не (X-1)*2. Такой инверсный ход времени во втором полупериоде в программе используется не всегда. Но он оказывается удобен для сшивания частей синусоиды при сильном наклоне её вершины - есть такая возможность.   
  Вообще, именно в деформации синусоиды и заключается последующая работа для получения, например, фортепьянного звука.


   ОБЩАЯ СТРУКТУРА ФУНКЦИИ  FNS1

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

   К числу этих действий относится также и повторный просмотр сценариев, в котором ищется не положительный номер K, а такой же отрицательный номер. Добавка к сценарию строк CASE с отрицательным номером представляет собой онлайн коррекцию развития звука. Любой из параметров может быть с её помощью изменён. Например строка
  IF T>.5 THEN ... позволит поправить сценарий по прошествии полсекунды после начала развития звука. Или IF NT=40 THEN ... - вносит коррективы на 40-ом периоде развития звука.
  Корректировки сценариев с отрицательными CASE записываются каждая под своим сценарием, чтобы было удобно обозревать и править всё в одном месте (см. текст программы).

  Функция создания звука представлена здесь в упрощённом виде. В ней отсутствует обертонная часть и отсутствуют обращения к ревербератору, создающему эхо. Об обертонной части рассказывалось в прежних статьях, я не планирую её существенно переделывать. О ревербераторе напишу позднее. Вместе с тем сама представленная здесь функция несколько изменена и оптимизирована. Возможностей у ней стало больше, а рабочих групп стало меньше - не 5, а только 3.
  Каждая из групп представлена единообразно - модуляционным, амплитудным и диссипативным параметрами, но возможности каждой из групп различаются:

 MU1,AU1,DU1 - группа, ориентированная за звук фортепьяно
 MU2,AU2,DU2 - универсальная группа, создающая синусоиду, двугорбый профиль, треугольник или звук виолы, в зависимости от значения модификатора mod2=0,1,2,3.
 MU3,AU3,DU3 - группа, воспроизводящая форму колебания по заданному рельефу.
   
 Параметры AU представляют собой значения амплитуды колебания
 Целая часть MU - число периодов колебаний, приходящихся на один период амплитудного вибрато, а дробная часть - глубина амплитудного вибрато
 Целая часть DU3 - число периодов, во время которых происходит атака - линейное нарастание звука от нуля до амплитудного значения.
Дробная часть DU3 - коэффициент затухания, начинающегося после окончания атаки.

  Каждая рабочая группа создаёт свою отдельную амплитуду A1,A2 и A3, перед выходом из программы эти амплитуды складываются.


  РАБОТА С ГРУППОЙ ПАРАМЕТРОВ MU1,AU1,DU1

  Звук фортепьяно создаётся по сценарию 2 с онлайн коррекцией -2:

 CASE 2
    MU1=64.1 : AU1=2 : DU1=3.995 : EA=.6 : E0=1.25 : E1=1.3 : E2=1.3
    Apm$="505050707050505560657035504744433857" : Dpm=.95
  CASE -2 : E0=1.25*(1+.2*FNTIM(4,T))
    E1=1.3*(1+.2*FNTIM(10,T)) : E2=1.3*(1+.2*FNTIM(7,T))

  Группа параметров указывает на лёгкую неторопливую вибрацию - MU1=64.1, резкую атаку и медленное затухание - DU1=3.995 и организуют вычисление синусоиды с амплитудой AU1=2.

  Далее идут коэффициенты пропорций EA, E0, E1 и E2, поправляющие форму вершины синусоиды, а также длительности полу и четверть периодов. Если синусоида нужна в чистом виде, то все E равны 1.
  EA поправляет вершину синусоиды - для звука фортепьяно нам нужна более заострённая вершина (о всех этих изменениях см. предыдущую статью "Создание звука фортепьяно").
Задание EA=.6 показывает, что только 60% звука берётся от синусоиды, а 40% звука добавляет заострённая форма колебания, для получения которой значения синусоиды возводятся в 4-ю степень. Всё сказанное организуется строкой:

ZZ=V*V : A1=(1-ZZ) : A1=EA*A1+(1-EA)*A1*A1*A1*A1 : A1=AA1*A1*(1+Eo*V)

в последнем операторе которой присутствует ещё один коэффициент - Eo. Этот коэффициент организует наклон вершины колебания вправо или влево. Если наклон не нужен, то полагается Eo=0. В сценарии этого задания нет, поскольку оно уже было выполнено по умолчанию при нашем первом рабочем входе в функцию. Очень сильный наклон Eo=4 используется для получения переколебания - резкого и короткого выброса, такая форма может быть успешно использована для получения фортепьянного звука на низких октавах.

  Коэффициент E0>1 растягивает первый период колебания, перенося границу между полупериодами с координаты X=1 на координату E0
  Коэффициенты E1 и E2 делают аналогичные вещи с первым и со вторым полупериодами колебания, сказанное реализуется в предваряющих вычисления строках:

 IF X<E0 THEN V=X*2/E0-1 : B=1-E1 : ELSE V=(2-X)*2/(2-E0)-1 : B=1-E2
 IF V<-B THEN V=(V+B)/(1-B)  ELSE V=(B+V)/(1+B)

  В результате мы вместо чистой синусоиды получаем более заострённую и несколько изменённую форму колебания (см. рисунок). Звучит она вполне соответствуя звуку фортепьяно в первой октаве.
 
  Вместе с тем, полученный звук нельзя назвать удовлетворительным. Он слишком монотонный, от впечатления монотонности не избавляет даже лёгкое вибрато.
  В чём же дело? Проделаем эксперимент.
  Загрузим в редактор "Wave edit" естественный звук фортепьяно. На первый взгляд  колебания похожи на синусоиду, форма их м.б. чуть другая.
  Вырежем из звука десяток периодов и аккуратно растиражируем их. Затем прослушаем. И не услышим живого звука. Звук, составленный из "клонов", окажется монотонным.
  Следовательно, форма реального звука не только отличается от формы синусоиды, но она ещё и эволюционирует, медленно меняясь от периода к периоду.

  Но ведь то же самое можем сделать и мы. Вы помните про участок программы в слежах? Он как раз и предназначен для осуществления подобных медленных изменений. От периода к периоду. Из него мы обращается на поправочную часть сценария CASE -2 и с помощью периодической функции FNTIM проводим вариацию параметров E0,E1 и E2.
Функция FNTIM имеет первым параметром частоту, а вторым параметром - время. Она возвращает значения колебания треугольной с амплитудой от -1 до +1
  Мы используем разные частоты - 4, 10 и 7 Гц.
  В отношении выбора частот - более 2-х Гц, иначе эффекта на фоне длительности звучания ноты не заметим, но менее 15 Гц, иначе будет слышаться вибрация наподобие звучащей буквы "Р". Кроме того частоты не должны быть кратны друг другу, иначе будут заметны биения звука. всё сказанное легко проверяется экспериментально.
Экспериментально подбираются и амплитуды отклонения параметров.
  У нас взяты одинаковые амплитуды .2, специально утрированно большие, чем это требуется - для того, чтобы продемонстрировать создаваемый эффект.

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

  Но не только деревянные стуки составляют феномен удара. Многочисленные и быстро затухающие обертоны слышатся в первый момент. Наиболее сильны 4-й, 5-й и 6-й обертоны. Эти номера обусловлены расположением точки удара молоточка по струне.
Короткий тон, на октаву или квинту более низкий, чем основная частота в сочетании с высокими обертонами привносит в звук удара металлический призвук.

  Но не только стуки и обертоны знаменуют начальную фазу звучания струн после удара.
  Бывают заметны и не гармонические колебания. В окне редактора"Wave edit" можно видеть, что иногда несколько периодов колебания развивающегося звука без понятных причин смещаются то вверх, то вниз.
  Ну хорошо, раз такой феномен замечен, то почему бы его не воспроизвести?
  Специальная строка Apm$ (толкай в плюс-p, толкай в минус-m) позволяет это сделать.

   Apm$="505050707050505560657035504744433857" : толкай в плюс-p, толкай в минус-m

  Запись в строке расшифровывается просто - представьте звуковое колебание человеком, имеющим положительную верхнюю часть - голову и руки, и отрицательную нижнюю часть - ноги. И вот этот человек стоит по пояс в воде и наружу торчат 50% всей его длины - видите число 50 в начале строки? Так он стоит первый, второй и третий период. Потом залезает на камешек, и наружу торчат уже 70% всей его длины. Постоял 2 периода на камешке и снова прыгнул на дно. Думаю, принцип построения строки "толкай в плюс, толкай в минус" Вам понятен.
  В действительности всё происходит несколько не так. Чтобы избежать разрывов на нулевой линии, положительная и отрицательная части колебания не перемещаются вверх вниз, а изменяют свой вертикальный размер.

  Я попробовал - даже простой тембр синусоиды с введением процесса "толкай вверх вниз" приобрёл рояльную живость.
  Диссипативный фактор Dpm=.95 обеспечивает быстрое затухание процесса - вызванный процессом посторонний звук низкого тона должен быть коротким, чтобы не стать заметным.
  Строка Apm$ может использоваться не только при демонстрации начальной фазы удара. Незатухающий процесс "толкай вверх вниз" при небольших его амплитудах очень хорошо оживляет синтезированный звук аккордеона.

  Прослушать постепенное становление звука фортепьяно можно на этой записи –
  https://yadi.sk/d/u8Sg4vB7pvkEP  а характер колебаний увидеть на рисунке.

Первые 4 звука - чистая синусоида с короткой атакой и долгим затуханием
Вторые звуки сделаны с привлечением параметров E, то есть с изменением формы колебания
В третьей группе звуков добавлены вибрации параметров E
В четвёртой группе ко всему этому добавлена и амплитудная вибрация MU1=64.1
В следующих 4-х звуках добавлен процесс "толкай вверх вниз"

  И, наконец, последние 4 звука сделаны с добавлением обертонных колебаний обертонов 3 и 4. Сделано это добавлением в сценарий CASE 2 строк, использующих 3-ю группу параметров (об этой группе речь пойдёт ниже):

     MU3=0 : AU3=1.5 : DU3=2.99
    '>>>>="_1_2_3_4_5_6_7_8_9_0_1_2"
     A3r$="905030307050505070303050"

  Обертоны заметно начинают звучать при достаточно сильном ударе. Последние 4 звука записи производят именно такое впечатление.
  В полном своём варианте функция создания звука включает в себя обертонную часть. С её помощью обертоны создаются проще. Здесь же я просто хотел показать возможности третьей группы параметров.


  РАБОТА С ТРЕТЬЕЙ ГРУППОЙ ПАРАМЕТРОВ MU3,AU3,DU3

  Группа параметров MU3,AU3,DU3 работает со строкой A3r$, задающей рельеф звукового колебания внутри периода. Мнемоника этой строки сходна с мнемоникой строки "толкай вверх вниз". Взгляните на сценарий CASE 3. Числа 70 99 50 00 в строке A3r$ задают внутри периода 4 точки, по которым подпрограммой M1A3 делается сплайн интерполяция. На краях сплайна слева и справа предполагаются значения 50, соответствующие нулевой амплитуде.
  Помните? Это как уровень воды, в которой человек стоит по пояс.

  CASE 3
    MU3=64.1 : AU3=2 : DU3=3.995
    A3r$="70995000"
  Apm$="505050707050505560657035504744433857" : Dpm=.95

  Строка A3r$="70995000" в этом сценарии почти точно воспроизводит форму фортепьянного звука, о котором речь шла выше. Но там эта форма создавалась параметрами группы 1.
Если несколько варьировать числа в строке перед началом каждого периода, то можно сделать звук более реалистичным и не монотонным. В принципе хорошо, что программа предоставляет различные возможности для создания формы колебаний.

  Сценарий CASE 4 создаёт звук аккордеона. Форма полученного звука показана на рисунке.

  CASE 4
    MU2=0 : AU2=1 : DU2=30.999 : mod2=2
     MU3=0 : AU3=3 : DU3=30.999
             Apm$="505152504546504651555047504751535057" : Dpm=1

  Здесь использованы одновременно две группы параметров. Группа 2 создаёт треугольную форму, а группа 3 добавляет пички. Чтобы пички получались острыми, использовалась длинная строка A3r$, содержащая 24 числа. А чтобы не запутаться в расстановке чисел, сверху к строке приставлена разметка в форме комментария - всё, что находится правее апострофа в Турбо Бейсике рассматривается как комментарий.

  Сделаю небольшое отступление общего плана. Как Вы понимаете, чтобы практически озвучивать мелодии, этих особых изысков со строками и с различными группами параметров не нужно. Нужно иметь большой числовой массив, в котором записаны формы звука на все случаи жизни, ходить по этому массиву и делать сплайн интерполяцию. Но это, во-первых, - скучно, а во вторых не позволяет экспериментально исследовать звук. Чтобы продвигать наши знания о музыкальном звуке вперёд, нужны гибкие инструменты. 

  Группа параметров MU2,AU2,DU2 помимо универсальной поддержки других групп, имеет и свою специализацию. Она "заточена" на создание звука смычковых инструментов. В этой статье группа Виолы показана не полностью, и о ней только упоминается. Широкие и интересные возможности для смычкового звука предоставляет ревербератор. О нём и о группе Виолы речь пойдёт в следующих статьях. Отдельно поговорим о флажолетах.
============
Дополнение.
Клип о создании звука фортепьяно см. тут - http://youtu.be/RF_JHuuOqck


  ТЕКСТ ФУНКЦИИ СОЗДАНИЯ ЗВУКА

'=============================
DEF FNS1(F,K)
STATIC FS,KS,T,DT,TT,PRENT
STATIC MU1,AU1,AA1,DU1,MU2,AU2,AA2,DU2,MU3,AU3,AA3,DU3
STATIC AK,DK1,mod2,Vi,Vo,EA,E0,E1,E2,Eo
STATIC Ap,Am,Apm$,Npm,Dpm,Mpm,A3r$
LOCAL X,NT,NK,A,A1,A2,A3,V,ZZ,B,B1,B2,L,Nr

 IF F=0 THEN FS=F : KS=K : FNS1=0 : EXIT DEF 'first setup
 IF F=FS AND K=KS THEN
       ELSE 'parameter reset
   FS=F : KS=K : T=0 : PRENT=-1 : DT=1/VEL : TT=1/FS
   MU1=0 : AU1=0 : DU1=1 : MU2=0 : AU2=0 : DU2=1 : MU3=0 : AU3=0 : DU3=1
   AK=0 : DK1=1 : mod2=0 : Vi=.35 : Vo=0
   Eo=0 : EA=1 : E0=1 : E1=1 : E2=1 'proportion for A1
   Apm$="50" : Npm=0 : Dpm=1 : Mpm=1 : A3r$="50"
   NK=KS
M1S1:
   SELECT CASE NK
  CASE 0
    MU1=0 : AU1=1 : DU1=1
  CASE 1
    MU1=32.2 : AU1=2 : DU1=INT(FS/40)+.995
    MU2=24.3 : AU2=.7 : DU2=100.999 : mod2=3
  CASE 2
    MU1=64.1 : AU1=2 : DU1=3.995 : EA=.6 : E0=1.25 : E1=1.3 : E2=1.3
    Apm$="505050707050505560657035504744433857" : Dpm=.95
   CASE -2 : E0=1.25*(1+.2*FNTIM(4,T))
    E1=1.3*(1+.2*FNTIM(10,T)) : E2=1.3*(1+.2*FNTIM(7,T))

  CASE 3
    MU3=64.1 : AU3=2 : DU3=3.995
    A3r$="70995000"
  Apm$="505050707050505560657035504744433857" : Dpm=.95

   END SELECT : IF NK<0 THEN RETURN
 END IF

  T=T+DT : X=T/TT : NT=INT(X) : X=(X-NT)*2
  IF PRENT<NT THEN '////////////////////////once at begin of period
   PRENT=NT
   AA1=AU1 : V=INT(MU1) : IF V>1 THEN AA1=AU1*(1+FNTIM(F/V,T)*(MU1-V))
   AA2=AU2 : V=INT(MU2) : IF V>1 THEN AA2=AU2*(1+FNTIM(F/V,T)*(MU2-V))
   AA3=AU3 : V=INT(MU3) : IF V>1 THEN AA3=AU3*(1+FNTIM(F/V,T)*(MU3-V))
   IF DU1>1 THEN V=INT(DU1) : AA1=(NT+1)*AA1/V : IF NT+1>=V THEN DU1=DU1-V
   IF DU2>1 THEN V=INT(DU2) : AA2=(NT+1)*AA2/V : IF NT+1>=V THEN DU2=DU2-V
   IF DU3>1 THEN V=INT(DU3) : AA3=(NT+1)*AA3/V : IF NT+1>=V THEN DU3=DU3-V
   IF DU1<1 THEN AU1=AU1*DU1
   IF DU2<1 THEN AU2=AU2*DU2
   IF DU3<1 THEN AU3=AU3*DU3
    Ap=Mpm*(VAL(MID$(Apm$,Npm*2+1,2))/50-1)+1 : Am=2-Ap
    Mpm=Mpm*Dpm : Npm=Npm+1 : IF Npm*2+1>LEN(Apm$) THEN Npm=0
   IF KS>0 THEN NK=-KS : GOSUB M1S1
  END IF '\\\\\\\\\\\\\\\\\\\\\\\

 IF X<E0 THEN V=X*2/E0-1 : B=1-E1 : ELSE V=(2-X)*2/(2-E0)-1 : B=1-E2
 IF V<-B THEN V=(V+B)/(1-B)  ELSE V=(B+V)/(1+B)
 ZZ=V*V : A1=(1-ZZ) : A1=EA*A1+(1-EA)*A1*A1*A1*A1 : A1=AA1*A1*(1+Eo*V)

 SELECT CASE mod2
  CASE 0 : A2=(1-ZZ) ' ~ sin  .^.
  CASE 1 : A2=(1+ZZ)*(1-ZZ*ZZ) 'camel  .^^.
  CASE 2 : A2=(1-ABS(V)) 'triangle   /\

  CASE 3  ' viola  /\___/\___/\___
     IF X<Vi THEN
       IF X<Vi/2 THEN A2=2*X/Vi ELSE A2=2*(Vi-X)/Vi
       ELSE : A2=0
     END IF
        A2=A2-Vi/4

 END SELECT : A2=A2*AA2
 GOSUB M1A3 'A3 take
 IF A3<0 THEN A3=A3*Am ELSE A3=A3*Ap
 IF X>E0 THEN A1=-A1 : IF mod2<3 THEN A2=-A2
 IF X>E0 THEN A1=A1*Am ELSE A1=A1*Ap

A=A1+A2+A3 : A=A+AK : AK=AK-(1-DK1)*A
FNS1=A
EXIT DEF
'=====
M1A3: 'spline interpolation
  IF AU3=0 THEN A3=0 : RETURN
  L=LEN(A3r$)/2 : IF L<1 THEN A3=0 : RETURN
  L=X*L/2 : Nr=INT(L) : L=(L-Nr)*2-1 : Nr=Nr*2+1
  B=VAL(MID$(A3r$,Nr,2))
  IF Nr=1 THEN B1=50 ELSE B1=(VAL(MID$(A3r$,Nr-2,2))+B)/2
  IF Nr=LEN(A3r$)-1 THEN B2=50 ELSE B2=(VAL(MID$(A3r$,Nr+2,2))+B)/2
  B2=(B2-B1)/2 : B1=B2+B1 'B2=(-)/2 B1=(+)/2
  ZZ=B1+B2*L : ZZ=ZZ+(B-B1)*(1-L*L)
  A3=AA3*(ZZ/50-1) : RETURN

END DEF
'=========
DEF FNTIM(F,Time) ' \/\/\/ low frequency vibration +-1 with F on Time
LOCAL T
  T=Time*F : T=T-INT(T) : IF T<.25 THEN T=.5-T
  IF T>.75 THEN T=1.5-T
  T=(T-.5)*4 : FNTIM=T
END DEF
'=============