CC BY
92
УДК 681.51: 534
Н.Ф. ТОДОРОВ
БЫСТРЫЕ АЛГОРИТМЫ В МАТЕМАТИЧЕСКИХ МОДЕЛЯХ АУРАЛИЗАЦИИ В АКУСТИКЕ ПОМЕЩЕНИЙ
Разрабатываются алгоритмы аурализации для компьютерного моделирования звучания заданного звукового файла в исследуемом помещении. Методом лучевых траекторий строится импульсный отклик помещения. Слышимый звук строится как свёртка исходного сигнала с функцией импульсного отклика. Для построения свертки используется быстрое преобразование Фурье (БПФ).
Ключевые слова: аурализация, преобразование Фурье, свертка, дифракционный интеграл, акустика помещений, метод лучевых траекторий.
Введение. Субъективный характер оценки качества звучания помещений приводит к тому, что расчёт «правильно звучащих» музыкальных залов представляет собой нетривиальную задачу. С появлением современных компьютеров её решить возможно новыми методами.
Появляется новое направление в архитектурной акустике - «Аурализация» (auralization). Его определение дал Мендель Клейнер (Mendel Kleiner) в 1989 году. Оно звучит так: "Аурализация - это процесс превращения звукового поля источника в пространстве в "слышимый звук" путем физического или математического моделирования таким образом, чтобы смоделировать бинауральное слуховое ощущение на заданной позиции моделируемого пространства» [1].
Существует ряд программ акустического моделирования - CATT, EASE, Odeon, Ulysses, и другие. Лидирующую программу выделить нельзя, набор акустических оценок более продвинут у CATT, а визуальная составляющая - у EASE. Среди недостатков этих программ - наличие определённой погрешности в вычислениях, что связано с привлечением статистических методов, а также некоторыми упрощениями (например, игнорирования в ряде случаев диффузии звука на рассеивающих поверхностях). Нельзя не отметить и то, что эти программные продукты - закрытые коммерческие разработки, созданные за рубежом [2].
В данной работе предлагается компьютерный алгоритм аурализации, основанный на методе лучевых траекторий. Предложенный алгоритм проанализирован на примере тестовых комнат. Описаны основные детали разработанного алгоритма.
Процесс формирования слышимого звукового сигнала. Рассмотрим некоторый источник звука, который создаёт определенный акустический сигнал, представляющий собой звуковую волну с зависимостью звукового давления от времени p0(t) . Если рассматривать помещение как линейный фильтр, который имеет свою амплитудно-временную энергетическую характеристику Er (t) = p2r(t), то в каждой точке пространства суммарный сигнал получается как "свертка" сигнала источника и характеристик помещения:
p(t) = (Ро ® Pr)(t) = J Po(t)Pr (t -t )dt ^ p(w ) = po(w )~r(w ). (1)
Самой сложной задачей при этом является определение полной спектральной характеристики отклика помещения рг (ю ) . Значения коэффициентов поглощения и диффузии, как правило, задаются для шести несущих октавных частот: f = 125; 250; 500; 1000; 2000; 4000 Гц. Таким образом, найти импульсную характеристику можно только для этих шести частот. В дискретном виде импульсный отклик помещения для каждой из этих частот ^ (рис.1) определяется некоторой дельтаобразной последовательностью (здесь 5 - дельта - функция Дирака)
р(пЛ/Е(ПГ 8 « - tf)), (2)
j = 1
-*( п )
где - моменты прихода очередного звукового луча к прием-
тт' ( п )
нику, а Ej - соответствующие значения энергии.
В частном случае, когда отклик помещения слабо зависит от частоты, форма аудиосигнала, преобразованного помещением, может быть выписана в явном виде:
J _____
р(/) = Х л/Е^Ро(/ - *}) ■ (3)
У = 1
Представление (3) имеет прозрачный физический смысл. Пусть импульсный отклик помещения один и тот же для всех частот звучания. Тогда исходный аудиофрагмент при воспроизведении в данном зале представляет собой сумму исходных копий этого фрагмента, сдвинутых на время запаздывания ]-го пришедшего луча и имеющих амплитуду, определяемую потерей энергии данного луча при всех его отражениях от граничных поверхностей.
Стоит отметить, что суммирование по формуле (3), осуществляемое во временной области, требует квадратичного числа арифметических операций по количеству выбранных временных узлов. Для существенного сокращения времени работы программы необходимо перейти в спектральную область, воспользовавшись для соотношения (1) теоремой свертки. В этом случае функция рг (ю ) является преобразованием Фурье от (2):
J _____
Рг (ю ) = Х л/Ё" ехр^ю tj ) , (4)
У = 1
однако быстрый алгоритм вычисления такого Фурье-преобразования может быть основан не на явной формуле, а на том или ином алгоритме типа БПФ. Такой подход является линейно-логарифмическим, т.е. практически линейным по числу узлов. Перемножение в спектральной области (1) является линейным, а Фурье-обращение функции р(ю) в (1) - линейно-логарифмическим [3].
В общем случае соотношение вида (4) можно использовать лишь для шести несущих октавных частот. Однако реальная зависимость от частоты здесь будет примерно такая же, как и частотная зависимость коэффициента поглощения, т.е. достаточно плавной функцией, что подробно
Рис.1. Импульсный отклик помещения
описано в литературе. В итоге можно достаточно эффективно интерполировать спектральные свойства на всю частотную область.
Практическая реализация этой идеи такова. Если ограничиться выбранным частотным диапазоном по шести основным октавным полосам, т.е. на интервале от 125/21/2 = 88 Гц до 4000-21/2 = 5640 Гц, то исходный сигнал р0^) можно представить в виде суммы шести сигналов, отфильтрованных каждый внутри своей октавной полосы:
po(t) =I pOn)(t) . (5)
n= 1
Если в первом приближении считать отклик помещения постоянным внутри одной октавной полосы, то с учетом (5) свертка (1) принимает следующий вид:
6 ¥ 6
P(t) = £ J Р0П) (t)p(и) (t -t )dt * ~(w ) = £ ~n) (w )~n) (w ), (6)
n= 1 -¥ n=1
(n w '•J mi
где импульсный отклик pr ’(t) в каждой октавной полосе свой. В итоге с использованием (5) получаем расчётную формулу
6 Jn ,------
p(t) =££ Vjp0n)(t-1r>, (7)
n= 1 j = 1
для быстрой реализации которой также применяется алгоритм БПФ.
Метод лучевых траекторий. В настоящее время при компьютерном расчете структуры звукового поля внутри замкнутого помещения получили развитие асимптотические методы, из которых самым эффективным оказался метод лучевых траекторий (МЛТ, Ray Tracing Method). Его использование даёт результаты, хорошо согласующиеся с натурными измерениями.
Для простоты ограничимся лишь помещением с плоскими граничными поверхностями. Его модель представляет собой конечный объем, ограниченный замкнутым многогранником. Имеется замкнутый многогранник с боковой поверхностью S . Внутри S необходимо решить волновое уравнение для акустического давления p( x, y, z, t) с некоторыми начальными условиями по времени и граничными условиями на внешних гранях:
1 Э2 p
DP = (*> У’ z) е V;
с Э t
З p З n
= g.p,(x, у, z) є s; (8)
і З p
plt=o - po;
= Vo
t = О
где V - внутренний объем многогранника £, т.е. £ = д V . Здесь с -скорость звука, п - внешняя нормаль к граничной поверхности; £ -произвольная і -я грань многогранника; импеданс уі на грани £і связан с поглощающими свойствами отделочного материала этой грани.
Удивительно, что в такой постановке не удалось построить ни одного реального компьютерного алгоритма, основанного на точном решении нестационарной краевой задачи (8). Это связано с тем, что длина волны на
s
порядки меньше характерных размеров отражающих поверхностей, и в результате требуется настолько большое число узлов сетки, что точный расчет не удается реализовать в реальном масштабе времени даже на самых быстрых современных компьютерах.
В связи с этим реальные алгоритмы могут быть основаны на слежении за траекториями звуковых лучей. Данные идеи связаны с геометрической теорией дифракции и излагаются здесь для случая произвольного числа отражений от любых гладких отражателей.
Итак, пусть из точки х0 акустической среды излучается высокочастотная сферическая волна. Считаем, что далее акустическая волна рас* * * * пространяется вдоль луча хо - У1 - У2 - Уз - ••• - Ум - хм + 1, где точки
* * * *
зеркального отражения у1,у2,у3,...,ум принадлежат граничным в общем случае искривленным поверхностям различных N отражателей.
Волна принимается в точке хм +1 акустической среды, в которой требуется найти амплитуду N раз переотраженной волны. Последняя определяется отражением волны от малых окрестностей £*, £2,-,гра-у * * * * ничных поверхностей в точках зеркального отражения У1 , У2, Уз ,•••,Ум .
Давление в N раз отраженной волне в точке хм +1 будем находить в виде интеграла Кирхгоффа по окрестности £* последней точки зеркального отражения у*м лучей, полученных при однократном отраже-
С'* * ^
нии от окрестности -1 предпоследней точки зеркального отражения У N-1. Тогда давление в точке приема дается выражением
( ) ГГ 2 ( ) ^ Ф ( У N , XN + 1 )
Р( Х N + 1 )= II 2 Р ( УN ) -----------^------------- dSN
'N1 N + д П N
(9)
где функция Ф (г) = ехР(ікг)/(4р г) обозначает функцию Грина для оператора Гельмгольца [4].
В свою очередь, давление р (ум) само выражается в виде интегрального представления через падающую на окрестность £* волну, при-
^ *
шедшую после отражения на окрестности N - 1 :
( Ум -1, Ум )
Р (УN )= Я 2 Р (У N - 1 ) ^ ^ - 1
Э пм- 1
- 1
(10)
Продвигаясь последовательно в обратном направлении: хм+1 - уМ - ••• - у* - у* - х0, приходим к окончательному представлению
для Р +1) в виде 2N кратного интеграла:
р (х м +1 )= 2м Л Л • ЯЯр'" ^... _дФ ЭЭФ
■'* ** **** Э п1 Э п2 Эп м1 Эп м
о* о* о* о* 1 2 м-1 N
5М 5М - 1 ^2 51
. (11)
В случае волнового процесса на фиксированной круговой частоте а с волновым числом k = а /с в коротковолновой области можно использовать асимптотическое представление для функции Грина при
k ® ¥
Эф( ут-1> У т ) 1 (лт'У 1| Г1 'к\Ут - 1- Ут I Г, . (7-1 \1
-------------— = 1к СОБ У т-1 (4р) Ут-1 - Ут\ ? [1 +Щk )]
Э Пт-1
т = 1, 2, •• •, М + 1 , у0 = х0 , ум + 1 ХМ + 1 . (12)
Здесь ут-1 - угол между нормалью пт -1 и направлением падения луча Ут - 2 - Ут-1, т.е. угол падения в точке Ут-1, а | Ут-1 - Ут | -расстояние между точками пространства Ут-1 е 5* -1 и ут е Б*т . Заме-
т-1 т -
> * /->*
т - 1 т т ^ * *
тим, что падающий Ут-1 - Ут и отраженный Ут - Ут+1 лучи лежат в
*
одной плоскости с нормалью пт в точке ут . Обозначим расстояния
I * I т I * * \ т I* I т
| Х0 - У1 | = ^ 0 , | Ут Ут + 1 | = ^п , | УМ ХМ +1 | ^ N ,
т = 1м-1. После вынесения за знак интеграла медленно меняющихся функций имеем следующее интегральное представление для давления в точке приема:
( 'к ^М *
р (х* + 1 )= 2^ ^П ^>08/ п ЯЯ. Яе'к" ^•••^м-А
п= 1
5М 5М - 1
, (13)
Ф = | х0 - У1 | + | У1 - У2 | + ••• + | Ум - 1 - Ум К I Ум - ХМ + 1 I .
(14)
Тестирование, примеры и результаты вычислений. Решив главную задачу - построение импульсного отклика помещения - и используя вышеописанный алгоритм, можно преобразовать любой записанный (в идеале
- в безэховой камере) аудиофрагмент в форму, которую он примет при звучании в моделируемом зале. Сравнение расчетов на основе данного алгоритма с натурными измерениями даёт погрешность около 10%, что достаточно для многих случаев. На иллюстрациях показаны графики амплитудно-временных характеристик звуковых файлов, которые были обработаны на основе метода, предложенного в данной работе. На нижнем графике на рис.2 можно отчётливо увидеть несколько первых переотраже-ний, а также появляющийся реверберационный «хвост». На рис.3 показан аналогичный график для реального музыкального фрагмента.
Рис.2. Звук выстрела стартового пистолета, записанный в безэховой камере (верхний график) и его смоделированное звучание в тестовой комнате (нижний график)
Рис.3. Фрагмент музыкального произведения, записанного в безэховой камере (верхний график) и его смоделированное звучание в тестовой комнате (нижний график)
В качестве тестовой комнаты в рассмотренных двух примерах бралась простейшая модель прямоугольного помещения, одна из стен которого сделана со скосом для подавления стоячих волн. При этом коэффициент поглощения для тестирования алгоритма на звучании в гулком помещении специально выбирался малым, равным 0,1 для всех отражающих поверхностей.
Автор выражает признательность профессору М.А. Сумбатяну (Южный федеральный университет, г. Ростов-на-Дону) за ценные советы и постоянное внимание к работе.
Библиографический список
1. KuttruffK.H. Auralization of impulse responses modeled on the basis of Ray-Tracing results // J. Audio Eng. Soc.- 1993. - 41 - P. 876-880.
2. VorlanderM. Auralization. Aachen.: Springer, 2008.
3. Блейхут Р. Быстрые алгоритмы обработки сигналов. / Р. Блейхут.
- М.: Мир, 1989.
4. Макриненко Л.И. Акустика помещений общественных зданий. / Л.И. Макриненко. - М.: Стройиздат, 1986. - 173 С.
Материал поступил в редакцию 25.12.08.
N.F. TODOROV
FAST ALGORITHMS IN MATHEMATICAL AURALIZATION MODELS IN ROOM ACOUSTICS
We develop auralization algorithms for computer simulation of sounding of audio files in a certain hall. By using the Ray Tracing Method there is formed the impulse room response. The audible sound is constructed as a convolution of the initial signal with this impulse response. The computational realization is based upon FFT.
ТОДОРОВ Николай Федорович (р. 1982), аспирант кафедры теоретической гидроаэромеханики факультета «Математика, механика и компьютерные науки» Южного федерального университета. Окончил механико-математический факультет Ростовского государственного университета (2004). Область научных интересов - аурализация и акустика помещений.
Имеет 6 опубликованных научных работ. zire71@bk.ru
