Математическое программирование звука. WAV сэмплы

 ВВЕДЕНИЕ

  Существуют различные направления синтеза электронного музыкального звука. Одни из методов используют библиотеки качественно записанных природных звуков, другие полагаются на математическое описание или математическое моделирование процесса образования звука.

   Здесь я описываю собственные эксперименты в среде Турбо бейсика с функцией FNSS2(F,K), создающей звук в простой программе изготовления звуковых сэмплов.
   Как видно, от привычной функции типа Sin(F*t*2п) эта функция отличается тем, что вместо одного параметра, фазы колебания  X= F*t*2п функции передаётся два параметра. Первый параметр  - это частота F,  второй параметр  - K  является ссылкой на таблицу с номером  K, одну из нескольких таблиц, имеющихся в теле функции. Каждая таблица описывает форму первого полупериода  колебания и изменения формы при переходе к следующему колебанию. Второй полупериод колебания получается обращением первого по знаку и в противоходе по времени.
   Естественно, такое описание с сосредоточением усилий на одном полупериоде существенно ограничивает класс получаемых форм звукового колебания. Но его можно рассматривать как следующее приближение к реальному звуку после синусоиды, дающее вполне удовлетворительные результаты. Кроме того, на такой простой модели оказывается возможным прощупать многие нюансы восприятия звука.

   Как видно, в функции совершенно не фигурирует время. Не используется также и число Пи. Вместо этого при первом обращении к функции в начале генерации звука определяется шаг изменения фазы DX=2*F/VEL , где VEL сэмпл/сек - скорость записи. Период колебания при этом изображается отрезком (0,2) на оси X, первый полупериод (0,1), второй полупериод (1,2). Переменная X изменяется при каждом обращении к функции на DX, для неё вычисляется амплитуда - значение, возвращаемое функцией, а отрезок (0,2) проходится от начала и до конца за время, равное периоду одного колебания. Не совсем точно X попадает на конец интервала, но как только X становится большим двух, оператор X=X-2 возвращает переменную X в его начало.
 

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

    Вся функция FNSS2( F,K) состоит из следующих последовательно расположенных блоков:
- инициализация, установка переменных и параметров по умолчанию (при обращении с F=0)
- вычисление шага, установка фаз, и параметров по номеру K (при обращении с новым F или K)
- делается шаг X=X+DX, если X>2 то X=X-2, и делается переустановка формы колебания
- вычисление амплитуды колебания
- добавление к амплитуде обертонных колебаний
- отдельно, как дополнительная подпрограмма ко второму блоку - вычисление параметров для произвольно заданной частоты путём интерполяции между несколькими таблицами, каждая из которых соответствует ограниченному интервалу частот. Эта подпрограмма позволяет плавно перестраивать тембр звучания при переходе звука от самых низких частот до самых высоких.
Текст функции FNSS2( F,K) вы можете видеть ниже.

ПРИМЕЧАНИЕ: Прошёл год после написания этой статьи, и в звуковой функции произошли некоторые изменения. Для формирования звука активно используется не только первый полупериод, но и весь период звукового колебания. Кроме того звуковая функция работает, как одна из трёх, внутри программы озвучивания Трио.
  О структуре функции и о структуре программы Трио можно прочесть тут -
http://www.proza.ru/2016/03/20/1794
http://www.proza.ru/2016/03/28/350
  Кроме того, я отказался от идеи тотальной интерполяции всех параметров по частоте, как от громоздкой и не очень эффективной в своём воплощении, и передал коррекцию по частоте отдельных параметров в конкретные сценарии - корректируется только то, что нужно, и только там, где это нужно.
  Изменились также некоторые алгоритмы вычислений и наименования параметров, например обозначения MU1, AU1 и DU1 используются теперь вместо MZX3, AMZ3 и DAMZ3. А некоторые параметры оказались и вовсе излишними. Но общая идеология, о которой говорится ниже, осталась.


ОСНОВНОЕ КОЛЕБАНИЕ

   Итак.  в функции имеется как основное колебание, так и дополнительные к нему обертонные колебания. Поговорим об основном колебании. Формально оно описывается 5-ю тройками параметров, задаваемых таблично. Например так:
      MZX0=67 : MZX1=32 : MZX2=0 : MZX3=.064 : MZX4=0
      AMZ0=0 : AMZ1=2 : AMZ2=0 : AMZ3=1 : AMZ4=0
      DAMZ0=1 : DAMZ1=.999 : DAMZ2=1 : DAMZ3=.9998 : DAMZ4=1
по строкам - модуляционные параметры, определяющие изменения формы колебания
- амплитудные параметры, устанавливающие начальную амплитуду
- диссипативные параметры - коэффициенты уменьшения амплитуд для следующего колебания.

Имеются два вида формы колебания:
- парабола, заменяющая на интервале -1<v<1 косинус по формуле y=1-v^2  (1)
- более тупая функция (1+v^2)(1-v^4)             (2)
Обе функции могут смешиваться в соответствии с амплитудами  AMZ3 и AMZ4, их задающими. Как видно по таблице выше, ею определяется чистая парабола с начальной амплитудой 1.  Однако устанавливается эта амплитуда не сразу. Она постепенно нарастает от нуля до значения 1 за 1/4 от числа колебаний указанных в дробной части MZX3=.064 то есть через 64/4=16 колебаний. После этого амплитуда будет уменьшаться, перевычисляясь каждый раз по формуле AMZ3=AMZ3*0.9998.
   Для чего понадобилось такое странное задание числа 64? Для того, чтобы не плодить лишних параметров. Например,  MZX1=32 задаёт большее единицы значение модуляционного параметра, меняющего амплитуду AMZ1=2. Такое задание вызовет периодическое изменение этой амплитуды по закону косинуса с периодом в 64 раза большим. чем период основного колебания.
   Тройки параметров 2 и 3 определяют уклон вершины параболы или верблюда в ту или другую сторону. Делается это добавлением множителя (см. верхнюю табл.) y=(1+AMZ1*v). Огибающая полупериода в этом случае описывается формулой y=AMZ3*(1-v^2)(1+AMZ1*v), где AMZ1 синусоидально колеблется, а AMZ3 экспоненциально уменьшается.
   Надо сказать, что AMZ1 в пике своей амплитуды  со значением большим единицы создаёт на краях полупериода выброс в отрицательную область, который тут же подхватывается положительным выбросом следующего полупериода. Такое сочетание хорошо воспринимается на слух в басах, где гладкие синусоидальные формы колебания становятся слабо слышимыми.
   Дробные указания для модуляционных параметров 1 и 2 интерпретируются немного по-другому.
Если для 3 и 4 амплитуда в начальной стадии колебания изменялась от нуля и до максимального значения, даже несколько переходя его: v =>  от 0 до 0.6 , то для 2 и 3 амплитуда проходит максимум и ослабевает: v =>  от 0.47 до 0.7 (вычисления ведутся всё по той же формуле перевёрнутой "рожками вниз" параболы). Где, и как именно это делается вы можете посмотреть по тексту программы.
    Осталось сказать о тройке параметров 0. Выше говорилось, что при каждом обращении к функции переменная X изменяется на DX. Это так, да не совсем. Если, например, AMZ0=0.1, то шаг по первому полупериоду будет на 10% короче, а сам полупериод стане фактически на 10% длиннее, зато со вторым полупериодом произойдут обратные вещи. Результат ощущается на слух как появление tremolo.
   Для модуляционного параметра  MZX0 не предусмотрено отдельной интерпретации дробной части. MZX0 задаёт период колебаний.

ДОПОЛНЕНИЕ К ОСНОВНЫМ КОЛЕБАНИЯМ

  Как показала практика, форма колебания «Верблюд» с параметрами 4 не является сильно востребованной.  Поэтому введён модификационный флаг mod4. Если он установлен в положение 1 (значение по умолчанию), то формируется колебание «Верблюд», если флаг установить в положение 2, то в качестве формы колебания будет взят равнобедренный треугольник.
Вы можете видеть это из текста программы:
  CASE 1 : AW2=AV*(1+ZZ)*(1-ZZ*ZZ)*(1+AZ2*V) 'camel
  CASE 2 : AW2=AV*(1-ABS(V))*(1+AZ2*V) 'triangle
Последние множители в этих операторах обеспечивают уклон вершины колебания в соответствии с параметрами 1 и 2, о чём было сказано выше.

  В самом конце вычислений к амплитуде основного колебания AW прибавляется амплитуда колебания обертонов AWo:
FNSS2=AW+AWo*oH
  - поскольку для обертонов не предусмотрено мер к постепенному возрастанию амплитуды каждого, то введён множитель oH (по умолчанию сразу равный 1 ), возрастающий от первоначально заданного значения до единицы с шагом doH ( doH=.01 по умолчанию) постепенно, на каждый период колебания. Такая техника даёт возможность постепенно включать звучание обертонов, и делать это независимо от основного колебания.

ОБЕРТОНЫ

   Об основном колебании сказано всё. Перейдём к обертонной части функции.
   Перед выходом из функции вы видите обращения к подпрограммам создающим три основных обертонных звука с адресами 21, 22 и 23. Три последующие обращения тоже можно использовать, но они задуманы как вспомогательные - по мере того как частота звука F будет переходить из одной области в другую, будут плавно меняются все те параметры, о которых говорилось выше.
   Для обертонов, однако, параметры являются дискретными величинами, они не могут меняться плавно. Поэтому предполагается не изменять параметры обертонной группы, а плавно заменять звучание одной группу на звучание другой, изменяя амплитуды этого звучания.

   Изменение параметров звучания в зависимости от частоты в большинстве случаев должно делаться обязательно. Звук похожий на звук струны в первой или во второй октавах, звучащий при тех же параметрах  в октавах малой или большой, вполне может походить на звук колокола. Звук фортепьяно имеет по меньшей мере пять характерных областей, для каждой из которых находятся свои оптимальные параметры - большая октава и ниже, малая, первая, вторая, третья и выше. В упрощённом виде, приведённом ниже, достаточно трёх областей - CASE 1, 2 и 3.
  После многочисленных опытов с представленной здесь функцией, мне показалось, что с ней может быть получен удовлетворительный звук фортепьяно, звук гитарной струны, звук колокола и звук варгана. Но есть звуки, которые даже отдалённо имитировать не удавалось – скрипка, аккордеон, труба.
  Для этих звуков возможно нужно другое решение, и другие функции, может быть более богатые на обертоны, или выходящие за пределы программирования одного полупериода колебания.
  И действительно, звук аккордеона удалось получить именно на этом пути -  CASE 8. Пояснения к опциям используемых там обертонов даны ниже.

  После этого отступления давайте вернёмся к обсуждению имеющейся группы трёх обертонов. Также как и для основного тона для каждого обертона имеются три параметра - модуляционный, амплитудный и диссипативный. Вот пример задания этих параметров:
      MZX5=3 : AMZ5=1 : DAMZ5=.99
 здесь 3 указывает на 3-й обертон (частота 3F), амплитуда = 1, после каждого периода основного колебания (на частоте F) амплитуда уменьшается на 100-99=1 процент. Звук этого обертона складывается со звуком основного колебания.
  А попробуем взять полуцелый М-параметр
      MZX5=1.5 : AMZ5=1 : DAMZ5=1
что получится, вы видите на рисунке "обертон 1.5" над заголовком статьи.
  Отрицательные М-параметры имеют особую интерпретацию. Если в строке выше указать  MZX5=-2 , то получим колебание "обертон -2". Для чего нужен такой обертон? Несмотря на то. что пики в нём расположены вдвое чаще основного колебания, никакого высокого звучания он не создаёт. Он звучит по-прежнему на частоте основного звука, но с другим тембром. Дело в том, что те пички, которые идут в нижней части колебания, действуют в точности в противофазе к верхним пичкам, таким образом, если ощущение чего-то колебательного с вдвое большей частотой и возникает внутри нашего уха, так оно тут же и гасится. Назовём подобные обертоны квазиобертонами.
  А вот если при такой же частоте пичков задать отрицательным Д-параметр:
       MZX5=-2 : AMZ5=1 : DAMZ5=-1      см. "обертон -2 Д -1"
мы услышим звук вдвое меньшей частоты, на октаву ниже того звука, который задан частотой  F.

   Наверное, плохо задавать М-параметры квазиобертонов дробными числами? Оказывается, совсем не плохо - надо только предусмотреть, что делать в том случае когда кривая квазиобертона, начавшись с нуля в начале полупериода не доходит до нуля в его конце.
Кстати, чем плохо, что не доходит? Плохо тем, что при переходе на второй полупериод образуется резкий разрыв в плавном ходе кривой, а значит в звуке будут слышны треск или шипение.
   Не будем доходить до этого, на последнем "приземлении в ноль" сделаем так. чтобы до конца полупериода квазиобертон не звучал совсем.
       MZX5=-1.99 : AMZ5=1 : DAMZ5=1  дают интересную кривую "обертон -1.99"
   Комбинация  MZX5=-1.99 ,  MZX5=-2.7 ,  MZX5=-3.4  даёт пилообразную кривую с глубоким звуком пароходного гудка (см. рисунок).

   Достаточно плавные и похожие на синусоиду кривые обертонов и квазиобертонов могут быть обострены, Если вместо MZX5=2 указывать 102, 202 и 302. Подобные тонкие пички квазиобертонов  можно видеть на нижней правой кривой рисунка, где они совмещаются с меняющими уклон параболой и верблюдом, которые нарастают по амплитуде, и смешаны между собой в одинаковой пропорции.
    
   Дополнительные обертонные возможности дают М-параметры вида 1002, 1301.03, 1101.25 и подобные. Единица в четвёртом порядке указывает на особую привязку вычислений, при которой не следует опасаться разрывов кривой при использовании не целочисленных значений. Вычисления обертонов привязываются не к переменной X, меняющейся с частотой F, а к переменной Xo, которая проходит свой период за 10 секунд. Потому при основном тоне До можно смело указывать не только на октавные обертоны высоких До, но и на близкое Ми -  MZX5=1001.25.
Хотите указать на Соль - задайте MZX5=1001.5
Не задавайте для обертонов такого вида отрицательных значений. Красивые картинки вы получите (см. самую верхнюю), но звучать они не будут. :))

  Имеется ещё два интересных параметра  D56 и M67. Первый параметр перекачивает амплитуду от обертона 6 к обертону 5, а второй периодически меняет амплитуды обертонов 6 и 7 местами. Они могут быть использованы, например, по такой схеме:
      MZX5=1 : AMZ5=0 : DAMZ5=.998
      MZX6=301 : AMZ6=14 : DAMZ6=1 : D56=98
      MZX7=0 : AMZ7=14 : DAMZ7=1 : M67=8
начинает звучать довольно громко и обострённо обертон 6. Обертон 5 начинает звучать более глубоким звуком по мере поступления энергии от 6-го (2% за каждый период). Но и 6-й не оскудевает - через каждые 8 периодов он обменяется энергией с 7-м. В результате образуется каплеобразный звук (по форме, а не по звучанию). Он сперва нарастает, слегка пульсируя, а потом постепенно ослабевает и меняет при этом тембр.

ОБЕРТОНЫ С ДОПОЛНИТЕЛЬНЫМИ ОПЦИЯМИ

  Выше говорилось о трудностях получения некоторых сложных по гармоническому составу звуков, таких как звук аккордеона. Тем не менее, введение обертонов с опциями 2000 позволило получить звук, похожий на звук аккордеона. Как это сделано, вы можете посмотреть по таблице CASE 8.
  Обертон опции 2000 ( например, 2002, 2301.5, 2202.5 ) заполняет обертонными колебаниями только первый, положительный полупериод основного колебания, тогда как обертон с отрицательной опцией -2000 заполняет только второй полупериод.
  И этого оказалось достаточно. Подробности можно посмотреть тут - http://www.proza.ru/2015/05/13/876

МНЕМОНИЧЕСКИЕ ОБОЗНАЧЕНИЯ

  Получаемые сэмплы звуков, чтобы они не перепутались, следует как-то маркировать. Удобно давать названия WAV файлам следующим образом:
   2CП0 1Д997-2-4Д99-104 3Д985-206-2Д98И_1 16Д992
   2RC - до большой октавы, частота 65.185 Гц
   П - парабола амплитудой 1, затухающая с фактором 0.997
к основному тону дополнительно прибавляются колебания - квазиобертоны:
   -2-4Д99 - 2 плавных мах на полупериод амплитуды -4 (вниз) с отражением их вниз на 2-й полупериод, затухание 0.99 на каждый период основного колебания
   -104 3Д985 - 3 обострённых мах амплитуды 3 (вверх) с затуханием 0.985
  -206-2Д98 - 6  острых мах амплитуды -2 (вниз) с затуханием 0.98
   И_1 16Д992 интервал 1-го полупериода шире номинального в 1.1, а ширина второго полупериода составляет 1-0.1=0.9  через 8 периодов соотношение интервалов обратное - (0.9)(1.1)  а ещё через 8 снова (1.1)(0.9). Здесь _1=0.1 - амплитуда расширения полупериода, 16=8+8 - периодичность процесса, Д992 = 0.992 - множитель уменьшающий амплитуду на каждый период основного колебания.

4CТ_26 1Д998о2301.5-2Д998-2302 1.5Д998оН=0 в соответствии с CASE 8 создаёт звук аккордеона.



4CDCcП_240 2Д9998 1003 1Д1 1003_05 1Д1 1003_03 1Д1
означает 4-я (первая) октава, ноты До,Ре,До и верхнее "до" исполняются со слегка, но не критично расстроенными друг относительно друга обертонами (звучит нота соль на октаву выше).
Основной тон - нарастающая парабола, амплитудой 2, затем медленно затухающая. Амплитуды обертонов равны 1.
Вместо десятичной точки в мнемонической записи иногда используется подчёркивание.

АЛГОРИТМ СОЗДАНИЯ ЗВУКА ПРОИЗВОЛЬНОЙ ЧАСТОТЫ

С помощью программы ZAW43.BAS создаётся набор звуков - сэмплы, наиболее подходящие для различных частот. Параметры, характеризующие отдельные частоты, в порядке повышения частоты указываются в наборах 1, 2, 3 и 4. В этих же наборах указывается и частоты им соответствующие. А именно -
F1- нижняя, F2-верхняя (F2=9999 указывает на все разумные частоты выше).
  В промежуточных случаях параметры основного колебания для Z0 и Z1 получаются линейной интерполяцией в соответствии с положением частоты F между F2 предыдущего набора и F1 следующего набора .
  Для трёх обертонов создаётся дополнительная тройка, абсолютные значения амплитуд в тройках получаются линейной интерполяцией (с сохранением знака амплитуд и с сохранением вида обертонов и сохранением затухания обертонов). Таким образом тройка предыдущих обертонов будет по мере повышения частоты постепенно переходить в тройку следующих обертонов, смешиваясь в совместном звучании.

  Программа  ZAW43.BAS при обращении к ней FNSS2(F,K) если K>=0 работает в штатном режиме - создаёт звук по указанному набору K. Однако, если K будет отрицательным, то программа обратится к набору -K, как первому из последовательности ряда описанных выше наборов. Если значения F попадут в интервал указанных там частот (F1,F2), то берутся указанные параметры. Если нет, то набор параметром запоминается и считывается следующий набор. В промежуточных случаях для определения параметров проводится интерполяция. 
  Текст функции, приведённый ниже содержит звуки, похожие на звуки фортепьяно. Чтобы получить звук частоты F нужно обратиться к функции, указав K=-1. Обращение с  K=6 делает красивую картинку бабочек с крылышками, которая к сожалению совершенно не звучит.
'====================
DEF  FNSS2( F,K)
STATIC 
FS,KS,X,DX,DAMZ0,AMZ0,DAMZ1,AMZ1,DAMZ2,AMZ2,DAMZ3,AMZ3,DAMZ4,AM
Z4
STATIC MZX0,MZX1,MZX2,MZX3,MZX4,AZ0,AZ1,AZ2,AZ3,AZ4
STATIC 
STATIC DZ0,DZ1,DZ2,DZ3,DZ4,Z0,Z1,Z2,Z3,Z4,AP,AV,F1,F2,PF1,PF2
STATIC PM0,PA0,PD0,PM1,PA1,PD1,PM2,PA2,PD2,PM3,PA3,PD3,PM4,PA4,PD4
STATIC PM5,PA5,PD5,PM6,PA6,PD6,PM7,PA7,PD7,PD56,PM67,Pn67
STATIC AWo1,AWo2,AWo3,AWo,Xo,DXo,BB,VV,PA,PB,PPA,PPB,mod4,oH,doH
LOCAL V,AW1,AW2,AW,DXM,XK,NK,B,BF,Ob,Obb

   IF F=0 THEN 'first setup
     FS=F : KS=K
     FNSS2=0 : EXIT DEF
   END IF

   IF F=FS AND K=KS THEN
           ELSE 'reset frequency, step, phases, amplitudes
       FS=F : X=2 : DX=2*FS/VEL : Z0=0 : Z1=0 : Z2=0 : Z3=1/2 : Z4=1/2
       Xo=0 : DXo=0.2/VEL ' 0.1 per second for obertones
       AZ=0 : AZ0=0 : AZ1=0 : AZ2=0 : AZ3=0 : AZ4=0
    MZX0=0 :  MZX1=0 :  MZX2=0 :  MZX3=0 :  MZX4=0
    AMZ0=0 : AMZ1=0 : AMZ2=0 : AMZ3=1 : AMZ4=0 : mod4=1 '1=> camel
    DAMZ0=1 : DAMZ1=1 : DAMZ2=1 : DAMZ3=1 : DAMZ4=1
    MZX5=0 : AMZ5=0 : DAMZ5=1 : MZX6=0 : AMZ6=0 : DAMZ6=1
    MZX7=0 : AMZ7=0 : DAMZ7=1 : D56=1 : M67=0 : n67=0 '3 obertones - 5,6,7
          AWo1=0 : AWo2=0 : AWo3=0 'bufer
    F1=0 : F2=0 :  PM0=0 : PA0=0 : PD0=1 : PM1=0 : PA1=0 : PD1=1
    PM2=0 : PA2=0 : PD2=1 : PM3=0 : PA3=0 : PD3=1 : PM4=0 : PA4=0 : PD4=1
    PM5=0 : PA5=0 : PD5=1 : PM6=0 : PA6=0 : PD6=1 : oH=1 : doH=.01
    PM7=0 : PA7=0 : PD7=1 : PD56=1 : PM67=0 : BB=0 : VV=-1 : PA=0 : PB=0
             'reset parameters
       KS=K  : IF KS>-.5 THEN NK=K : ELSE NK=-K
M1S2: ' <=========== for next NK
  SELECT CASE NK
  CASE 0
      MZX0=0 : MZX1=0 : MZX2=0 : MZX3=0 : MZX4=0
      AMZ0=0 : AMZ1=0 : AMZ2=0 : AMZ3=1 : AMZ4=0
      DAMZ0=1 : DAMZ1=1 : DAMZ2=1 : DAMZ3=1 : DAMZ4=1

   CASE 1
      MZX0=16 : MZX1=.256 : MZX2=.256 : MZX3=.04 : MZX4=.04 : F1=0 : F2=66
      AMZ0=.03 : AMZ1=4 : AMZ2=4 : AMZ3=.5 : AMZ4=.5
      DAMZ0=1 : DAMZ1=.98 : DAMZ2=.98 : DAMZ3=.9997 : DAMZ4=.9997
      MZX5=-304 : AMZ5=.8 : DAMZ5=.998 : VV=-1 : PA=0 : PB=0
      MZX6=0 : AMZ6=0 : DAMZ6=1 : D56=1
      MZX7=0 : AMZ7=0 : DAMZ7=1  : M67=0 : n67=0
   CASE 2
  MZX0=53 : MZX1=0 : MZX2=0 : MZX3=.16 : MZX4=0 : Z0=1/2 : F1=135 : F2=520
  AMZ0=.12 : AMZ1=0 : AMZ2=0 : AMZ3=2 : AMZ4=0 : Z1=1/2 : Z4=0
  DAMZ0=.998 : DAMZ1=1 : DAMZ2=.95 : DAMZ3=.998 : DAMZ4=.998
      MZX5=201 : AMZ5=1.6 : DAMZ5=.998 : PA=.27 : PB=-.1
      MZX6=1004 : AMZ6=.4 : DAMZ6=.99  ': VV=-1
      MZX7=1005 : AMZ7=.4 : DAMZ7=.99
   CASE 3 : AMZ0=.2 : MZX5=0 : MZX3=.02 : F1=700 : F2=999999

  CASE 4 'save f-no sound - звук ф-но в среднем диапазоне
      MZX0=57 : MZX1=0 : MZX2=0 : MZX3=.16 : MZX4=69 : Z0=1/2
      AMZ0=.12 : AMZ1=0 : AMZ2=0 : AMZ3=1 : AMZ4=0 : Z1=1/2 : Z4=0
      DAMZ0=.998 : DAMZ1=1 : DAMZ2=.95 : DAMZ3=.998 : DAMZ4=.998
      MZX5=201 : AMZ5=.8 : DAMZ5=.998 : PA=.27 : PB=-.1
      MZX6=1004 : AMZ6=.2 : DAMZ6=.99 : VV=-1
      MZX7=1005 : AMZ7=.2 : DAMZ7=.99
   CASE 5 'красивый звук с обертонами для "Полета шмеля" хотя и не скрипка
      MZX0=53 : MZX1=0 : MZX2=0 : MZX3=.06 : MZX4=69 : Z0=1/2
      AMZ0=.2 : AMZ1=0 : AMZ2=0 : AMZ3=1 : AMZ4=0 : Z1=1/2 : Z4=0
      DAMZ0=.998 : DAMZ1=1 : DAMZ2=.95 : DAMZ3=1 : DAMZ4=.998
      MZX5=2008 : AMZ5=.02 : DAMZ5=.998   
      MZX6=2006 : AMZ6=.02 : DAMZ6=.998 : VV=-1
      MZX7=2004 : AMZ7=.02 : DAMZ7=.998
 
  CASE 6 ' рисует красивую верхнюю картинку, но не звучит :))
      MZX0=67 : MZX1=.04 : MZX2=0 : MZX3=4 : MZX4=8
      AMZ0=0 : AMZ1=0 : AMZ2=0 : AMZ3=0 : AMZ4=0
      DAMZ0=1 : DAMZ1=.7 : DAMZ2=1 : DAMZ3=1 : DAMZ4=1
      MZX5=-1002 : AMZ5=1 : DAMZ5=1
      MZX6=1312 : AMZ6=1 : DAMZ6=1
  CASE 7 ' свистящий звук
      MZX0=0 : MZX1=0 : MZX2=0 : MZX3=.26 : MZX4=.06 : Z0=1/2
      AMZ0=.4 : AMZ1=0 : AMZ2=0 : AMZ3=1 : AMZ4=0 : Z1=1/2 : Z4=0
      DAMZ0=.998 : DAMZ1=1 : DAMZ2=.95 : DAMZ3=.998 : DAMZ4=.998
      MZX5=-4 : AMZ5=2 : DAMZ5=.998 : PA=-.428 : PB=0
      MZX6=-30 : AMZ6=.7 : DAMZ6=.998 : VV=-1
      MZX7=0 : AMZ7=.2 : DAMZ7=.998
  CASE 8  'звук аккордеона
      MZX0=0 : MZX1=0 : MZX2=0 : MZX3=.26 : MZX4=.26 : Z0=1/2 : Z2=1/2
      AMZ0=0 : AMZ1=0 : AMZ2=0 : AMZ3=0 : AMZ4=1 : Z1=1/2 : Z4=0 : mod4=2
      DAMZ0=1 : DAMZ1=1 : DAMZ2=1 : DAMZ3=.998 : DAMZ4=.998
      MZX5=2301.5 : AMZ5=-2 : DAMZ5=.998
      MZX6=-2302 : AMZ6=1.5 : DAMZ6=.998 : VV=-1 : oH=0
      MZX7=0 : AMZ7=0 : DAMZ7=.998

 END SELECT
 NK=-NK : IF KS<0 THEN GOSUB 29
 IF NK>0 THEN GOTO M1S2

      DZ0=0 : IF MZX0>0 THEN DZ0=2/MZX0
  DZ1=0 : IF MZX1>0 THEN DZ1=2/MZX1 : IF MZX1<1 THEN Z1=.45 : DZ1=DZ1/1000
  DZ2=0 : IF MZX2>0 THEN DZ2=2/MZX2 : IF MZX2<1 THEN Z2=.45 : DZ2=DZ2/1000
  DZ3=0 : IF MZX3>0 THEN DZ3=2/MZX3 : IF MZX3<1 THEN Z3=0 : DZ3=DZ3/1000
  DZ4=0 : IF MZX4>0 THEN DZ4=2/MZX4 : IF MZX4<1 THEN Z4=0 : DZ4=DZ4/1000
      END IF

  IF X<1 THEN
    V=DX/(1+AZ0) : IF X<.5 THEN DXM=V/(1+PA) ELSE DXM=V/(1-PA)
         ELSE
    V=DX/(1-AZ0) : IF X<1.5 THEN DXM=V/(1+PB) ELSE DXM=V/(1-PB)
  END IF
   Xo=Xo+DXo : IF Xo>2 THEN Xo=Xo-2
   X=X+DXM : IF X>2 THEN ' \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
   X=X-2
Z0=Z0+DZ0 : IF Z0>2 THEN Z0=Z0-2
IF Z0<1 THEN ZZ=Z0*2-1 : AZ0=(1-ZZ*ZZ)*AMZ0 ELSE ZZ=(Z0-1)*2-1 : AZ0=
(ZZ*ZZ-1)*AMZ0
  IF MZX1<1 THEN IF Z1>.7 THEN DZ1=0
Z1=Z1+DZ1 : IF Z1>2 THEN Z1=Z1-2
IF Z1<1 THEN ZZ=Z1*2-1 : AZ1=(1-ZZ*ZZ)*AMZ1 ELSE ZZ=(Z1-1)*2-1 : AZ1=
(ZZ*ZZ-1)*AMZ1
  IF MZX2<1 THEN IF Z2>.7 THEN DZ2=0
Z2=Z2+DZ2 : IF Z2>2 THEN Z2=Z2-2
IF Z2<1 THEN ZZ=Z2*2-1 : AZ2=(1-ZZ*ZZ)*AMZ2 ELSE ZZ=(Z2-1)*2-1 : AZ2=
(ZZ*ZZ-1)*AMZ2
  IF MZX3<1 THEN IF Z3>.6 THEN DZ3=0
Z3=Z3+DZ3 : IF Z3>2 THEN Z3=Z3-2
IF Z3<1 THEN ZZ=Z3*2-1 : AZ3=(1-ZZ*ZZ)*AMZ3 ELSE ZZ=(Z3-1)*2-1 : AZ3=
(ZZ*ZZ-1)*AMZ3
  IF MZX4<1 THEN IF Z4>.6 THEN DZ4=0
Z4=Z4+DZ4 : IF Z4>2 THEN Z4=Z4-2
IF Z4<1 THEN ZZ=Z4*2-1 : AZ4=(1-ZZ*ZZ)*AMZ4 ELSE ZZ=(Z4-1)*2-1 : AZ4=
(ZZ*ZZ-1)*AMZ4
 IF NOT(MZX1<1 AND DZ1>0) THEN AMZ1=AMZ1*DAMZ1
 IF NOT(MZX2<1 AND DZ2>0) THEN AMZ2=AMZ2*DAMZ2
 IF NOT(MZX3<1 AND DZ3>0) THEN AMZ3=AMZ3*DAMZ3
 IF NOT(MZX4<1 AND DZ4>0) THEN AMZ4=AMZ4*DAMZ4
 AMZ0=AMZ0*DAMZ0 : AMZ5=AMZ5*DAMZ5 : AMZ6=AMZ6*DAMZ6 : 
AMZ7=AMZ7*DAMZ7
IF D56<1 THEN AMZ5=AMZ5+AMZ6*(1-D56) : AMZ6=AMZ6*D56
IF M67>0 THEN n67=n67+1 : IF n67>M67 THEN n67=0 : SWAP AMZ6,AMZ7
IF PD56<1 THEN PA5=PA5+PA6*(1-PD56) : PA6=PA6*PD56
IF PM67>0 THEN Pn67=Pn67+1 : IF Pn67>PM67 THEN Pn67=0 : SWAP PA6,PA7
AP=AZ3 : AV=AZ4'*0.673077  '(4/3)/(208/105)
IF oH<1 THEN oH=oH+doH 'head of obertones
                END IF ' \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
IF X<1 THEN V=X*2-1  ELSE V=(X-1)*2-1 : V=VV*V
ZZ=V*V : AW1=AP*(1-ZZ)*(1+AZ1*V)
SELECT CASE mod4
  CASE 1 : AW2=AV*(1+ZZ)*(1-ZZ*ZZ)*(1+AZ2*V) 'camel
  CASE 2 : AW2=AV*(1-ABS(V))*(1+AZ2*V) 'triangle
END SELECT
AW=AW1+AW2 : IF X>1 THEN AW=-AW
IF X>1 THEN AW=AW/(1-AZ0) ELSE AW=AW/(1+AZ0)
GOSUB 21 : AWo=AZ*AMZ5
GOSUB 22 : AWo=AWo+AZ*AMZ6
GOSUB 23 : AWo=AWo+AZ*AMZ7
GOSUB 24 : AWo=AWo+AZ*PA5
GOSUB 25 : AWo=AWo+AZ*PA6
GOSUB 26 : AWo=AWo+AZ*PA7
'GOSUB 20 'подключение сглаживающего буфера
 FNSS2=AW+AWo*oH : EXIT DEF
'=======
20 AWo4=AWo3 : AWo3=AWo2 : AWo2=AWo1 : AWo1=AWo :
AWo=(AWo1+AWo2+AWo3+AWo4)/4 : RETURN
'=========
21 NK=MZX5 : IF NK=0 THEN AZ=0 : RETURN
      BF=1 : IF DAMZ5<0 THEN BF=0
      GOSUB 28 : RETURN
'=========
22 NK=MZX6 : IF NK=0 THEN AZ=0 : RETURN
      BF=1 : IF DAMZ6<0 THEN BF=0
      GOSUB 28 : RETURN
'=========
23 NK=MZX7 : IF NK=0 THEN AZ=0 : RETURN
      BF=1 : IF DAMZ7<0 THEN BF=0
      GOSUB 28 : RETURN
'=========
24 NK=PM5 : IF NK=0 THEN AZ=0 : RETURN
      BF=1 : IF PD5<0 THEN BF=0
      GOSUB 28 : RETURN
'=========
25 NK=PM6 : IF NK=0 THEN AZ=0 : RETURN
      BF=1 : IF PD6<0 THEN BF=0
      GOSUB 28 : RETURN
'=========
26 NK=PM7 : IF NK=0 THEN AZ=0 : RETURN
      BF=1 : IF PD7<0 THEN BF=0
      GOSUB 28 : RETURN
'=============
28    IF NK<0 THEN BF=0 : NK=-NK
      B=1 : Ob=0 : Obb=0 : IF NK>2000 THEN NK=NK-2000 : Obb=1
      IF NK>1000 THEN NK=NK-1000 : Ob=1
      IF NK>100 THEN NK=NK-100 : B=2 '=> _/""\_
      IF NK>100 THEN NK=NK-100 : B=3       '=> __/~\__
      IF NK>100 THEN NK=NK-100 : B=4       '=> ___/\__
      XK=X*NK/2
      IF Ob=1 THEN NK=F*NK : XK=Xo*NK*5 : IF BB>0 THEN B=BB
   IF Obb=1 THEN
     AZ=0
     IF X>1 AND BF=1 THEN B=0 'second 1/2period
     IF X<1 AND BF=0 THEN B=0 'first 1/2period
     IF X<1 AND NK*X>FIX(NK) THEN B=0
     IF X>1 AND NK*(X-1)>FIX(NK) THEN B=0
     BF=1 : IF X>1 THEN XK=(X-1)*NK/2
   END IF
   IF BF=0 THEN
     IF X<1 THEN V=X ELSE V=2-X
     XK=V*NK : XK=XK-FIX(XK) : ZZ=XK*2-1
   END IF
  IF BF=1 THEN XK=XK-FIX(XK) : ZZ=XK*4-1 : IF ZZ>1 THEN ZZ=ZZ-2
   ZZ=(1-ZZ*ZZ)
  SELECT CASE B
     CASE 1 : AZ=ZZ*ZZ : IF BF=1 THEN AZ=(ZZ*3+AZ)/4 '=> sin /"~"\
     CASE 2 : AZ=ZZ*ZZ*ZZ*ZZ
     CASE 3 : AZ=ZZ*ZZ*ZZ*ZZ*ZZ*ZZ
     CASE 4 : AZ=ZZ*ZZ*ZZ*ZZ*ZZ : AZ=AZ*AZ
  END SELECT
  IF BF=1 THEN IF XK>.5 THEN AZ=-AZ
  IF BF=1 THEN RETURN
  IF V>FIX(NK)/NK THEN AZ=0
  IF X>1 THEN AZ=-AZ
RETURN
'=====
29 IF F1<F AND F<F2 THEN PA5=0 : PA6=0 : PA7=0 : PD56=1 : PM67=0 : RETURN
     IF F2=0 THEN RETURN
     IF F>F2 THEN
   PM0=MZX0 : PA0=AMZ0 : PD0=DAMZ0 : PM1=MZX1 : PA1=AMZ1 : PD1=DAMZ1
   PM2=MZX2 : PA2=AMZ2 : PD2=DAMZ2 : PM3=MZX3 : PA3=AMZ3 : PD3=DAMZ3
   PM4=MZX4 : PA4=AMZ4 : PD4=DAMZ4 : PPA=PA : PPB=PB
   PM5=MZX5 : PA5=AMZ5 : PD5=DAMZ5 : PM6=MZX6 : PA6=AMZ6 : PD6=DAMZ6
   PM7=MZX7 : PA7=AMZ7 : PD7=DAMZ7 : PD56=D56 : PM67=M67
   PF1=F1 : PF2=F2 : NK=1-NK : F2=0 : RETURN
      END IF
   V=(F-PF2)/(F1-PF2) : B=1-V
MZX0=PM0+(MZX0-PM0)*V : AMZ0=PA0+(AMZ0-PA0)*V : DAMZ0=PD0+(DAMZ0
-PD0)*V
MZX1=PM1+(MZX1-PM1)*V : AMZ1=PA1+(AMZ1-PA1)*V : DAMZ1=PD1+(DAMZ1
-PD1)*V
MZX2=PM2+(MZX2-PM2)*V : AMZ2=PA2+(AMZ2-PA2)*V : DAMZ2=PD2+(DAMZ2
-PD2)*V
MZX3=PM3+(MZX3-PM3)*V : AMZ3=PA3+(AMZ3-PA3)*V : DAMZ3=PD3+(DAMZ3
-PD3)*V
MZX4=PM4+(MZX4-PM4)*V : AMZ4=PA4+(AMZ4-PA4)*V : DAMZ4=PD4+(DAMZ4
-PD4)*V
PA=PPA+(PA-PPA)*V  :  PB=PPB+(PB-PPB)*V
'MZX0=CEIL(MZX0) : MZX1=CEIL(MZX1)
' D56=PD56+(D56-PD56)*V : M67=PM67+(M67-PM67)*V : M67=CEIL(M67) ' ?
AMZ5=AMZ5*V : AMZ6=AMZ6*V : AMZ7=AMZ7*V
PA5=PA5*B : PA6=PA6*B : PA7=PA7*B
RETURN
        ' test
PRINT "==== F=";F
PRINT "0= ";MZX0;AMZ0;DAMZ0
PRINT "1= ";MZX1;AMZ1;DAMZ1
PRINT "2= ";MZX2;AMZ2;DAMZ2
PRINT "3= ";MZX3;AMZ3;DAMZ3
PRINT "4= ";MZX4;AMZ4;DAMZ4
PRINT "5= ";MZX5;AMZ5;DAMZ5
PRINT "6= ";MZX6;AMZ6;DAMZ6;" d56=";D56
PRINT "7= ";MZX7;AMZ7;DAMZ7;" m67=";M67
PRINT "P5= ";PM5;PA5;PD5;" PA=";PA;" PB=";PB
PRINT "P6= ";PM6;PA6;PD6;" Pd56=";PD56
PRINT "P7= ";PM7;PA7;PD7;" Pm67=";PM67
STOP
'=======
 END DEF
'=============

Текст списан с рабочей версии. Функция работает внутри программы создания WAV сэмплов, которая описана в предыдущей статье.  Вот ссылка на неё - http://www.proza.ru/2015/04/11/1943
Спасибо за внимание.