Využití interpolace pro simulaci elektroakustických měničů v reálném čase

Rok / Year: 2011

Svazek / Volume: 13

Číslo / Number: 3

Využití interpolace pro simulaci elektroakustických měničů v reálném čase Real-time simulation of electroacoustic devices using interpolation Jaromír Mačák12 , Vladimír Tichý1 [email protected],[email protected],[email protected] 1

2

Vědecko-výzkumné centrum JAMU Fakulta elektrotechniky a komunikačních technologií VUT v Brně.

Abstrakt: Simulace elektroakustických měničů nachází stále větší uplatnění při číslicovém zpracování hudebních signálů. Obvykle se simulace sestává z jednoduchého číslicového filtru navrženého podle změřené kmitočtové charakteristiky. Při simulaci elektroakustických měničů je však třeba brát v potaz i směrové charakteristiky měničů a to obzvlášť v blízkém poli, kde dochází ke snímání nástroje či reproduktoru. V tomto článku je prezentována metoda interpolace kmitočtové charakteristiky mezi jednotlivými měřenými polohami. Výsledkem interpolace je impulsní charakteristika, která je použita pro filtraci pomocí algoritmu rychlé konvoluce.

Abstract: Simulation of electroacoustic devices can be used in many applications for real-time digital audio processing. Usually, simulation is based on usage of simle digital filter designed according measured frequency response. However, during the simulation. directional characteristics have to be taken account, especially when simulating the near field of the electroacoustic device. The method based on usage of interpolation between measured frequency responces in near field of the electroacoustic device is presented in this paper.

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

Využití interpolace pro simulaci elektroakustických měničů v reálném čase Jaromír Mačák1,2 , Vladimír Tichý1 1 Vědecko-výzkumné centrum JAMU Email: [email protected],[email protected] 2 Fakulta

elektrotechniky a komunikačních technologií VUT v Brně Email: [email protected]

Abstrakt – Simulace elektroakustických měničů nachází stále větší uplatnění při číslicovém zpracování hudebních signálů. Obvykle se simulace sestává z jednoduchého číslicového filtru navrženého podle změřené kmitočtové charakteristiky. Při simulaci elektroakustických měničů je však třeba brát v potaz i směrové charakteristiky měničů a to obzvlášť v blízkém poli, kde dochází ke snímání nástroje či reproduktoru. V tomto článku je prezentována metoda interpolace kmitočtové charakteristiky mezi jednotlivými měřenými polohami. Výsledkem interpolace je impulsní charakteristika, která je použita pro filtraci pomocí algoritmu rychlé konvoluce.

1

Úvod

Frekvenční vlastnosti elektroakustických měničů zásadním vlivem ovlivňují výslednou barvu zvuku při procesu snímání a nahrávání hudebních nástrojů. Proto je volba správného typu a polohy mikrofonu klíčová vzhledem k odlišným charakteristikám jednotlivých typů. Občas však nastane situace, kdy se zjistí, že pro snímání nebyl použit správný mikrofon. V praxi to znamená opětovné pořízení nahrávky s novým mikrofonem nebo lze již pořízenou nahrávku zpracovat pomocí algoritmu, který simuluje chování daného mikrofonu. Nahrávku je tak možné upravit bez nutnosti nového nahrávání. Ještě komplikovanější situace nastává u snímání zvuku elektrické kytary, jak ukazují praktické zkušenosti zvukových mistrů, protože zde se na výsledné barvě zvuku podílejí dva elektroakustické měniče – mikrofon a reproduktor kytarového komba. Mikrofon je umístěn v blízkém poli reproduktoru, kde je akustické vyzařování reproduktoru velmi proměnné a frekvenční charakteristika reproduktoru je proto v každé pozici mikrofonu, kde je zvuk snímán, odlišná [5]. Obvykle se simulace kmitočtových charakteristik reproduktorů s ozvučnicí provádí pomocí elektroakustické a elektromechanické analogie na základě náhradního schematu. Parametry lze určit měřením [9]. Směrovou závislost kmitočtové charakteristiky v blízkém akustickém poli lze získat řešením rovnice kmitání kruhového pístově kmitajícího zařiče [5]. Tento způsob simulace, pokud má dávat

opravdu přesné výsledky, je výpočetně náročný a proto se příliš nehodí pro použití pro zpracování zvukových signálů, jelikož zpracování obvykle probíhá v reálném čase. Z pohledu zpracování zvukových signálů se elektroakustický měnič, pokud tedy simulujeme pouze lineární vlastnosti daného měniče, jeví jako prostý filtr s určitou frekvenční charakteristikou. Při číslicovém zpracování signálů lze linearizovaný elektroakustický měnič simulovat FIR nebo IIR filtrem s vhodným koeficienty [15]. Koeficienty lze získat měřením daného elektroakustického měniče, kdy ze změřené frekvenční charakteristiky lze vypočítat impulsní odezvu, jejíž vzorky jsou přímo koeficienty FIR filtru [8]. Pro návrh IIR filtru je třeba využít složitějších algoritmů, např. Yule-Walker algoritmus [2] nebo lze IIR filtr navrhnout tak, že aproximuje průběh frekvenční charakteristiky FIR filtru [3]. Výpočet a změření impulsní odezvy a následná filtrace je tedy poměrně snadnou záležitostí a tento postup se běžně využívá v komerčně využívaných algoritmech. Problémem však zůstává získání impulsní odezvy pro libovolnou polohu mikrofonu vůči snímanému reproduktoru, případně pro jakékoliv natočení mikrofonu, uvažujeme-li pouze snímaní akustického nástroje nebo hlasu mikrofonem. Pro zajištění plynulé změny frekvenční charakteristiky by tedy muselo být naměřeno velké množství impulsních odezev, teoreticky až nekonečně mnoho. Vyjdeme-li ovšem z hypotézy, že naměřené frekvenční charakteristiky a tedy i impulsní odezvy v sousedních polohách mají témět stejný průběh, lze naměřit pouze impulsní odezvy ve více vzdálených polohách a impulsní odezvy odpovídající polohám mezi těmito polohami dopočítat pomocí interpolace.

2

Měření kmitočtových charakteristik elektroakustických měničů

Při měření elektroakustických měničů je nutno dodržet několik podmínek. V případě měření frekvenčních charakteristik jsou největší požadavky kladeny na akustické prostředí, ve kterém bude elektroakustický měnič měřen. Prostředí musí být bezodrazové a musí mít minimální, ideálně nulovou, zbytkovou hladinu hluku. Měření může probíhat ve volném nebo difuzním poli, přičemž realizace difuzního pole bývá obvykle finančně nákladná a proto se častěji volí

25 – 1

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

varianta volného pole. Volným polem může být otevřený prostor, což je bezodrazové prostředí, nevýhodou je však vyšší úroveň zbytkových hladin hluku. Tento nedostatek odstraňuje tzv. bezodrazová komora, což je také bezodrazové prostředí a je zde i minimální úroveň zbytkových hladin hluku. Nevýhodou je nepřesnost měření na nízkých frekvencích díky nedokonalé absorpci způsobené omezenou velikostí absorpčních bloků [1]. Měřící soustava pro měření kmitočtových charakteristik reproduktorových boxů je poměrně jednoduchá. Sestává se z budicího zesilovače měřeného reproduktorového boxu, měrného (kalibrovaného) mikrofonu a mikrofonního předzesilovače. Roli generátoru budicího signálu a zároveň záznamového zařízení může hrát osobní počítač s vhodným software, samozřejmě lze využít i specializovaných měřících zařízení. Na budicí zesilovač a předzesilovač jsou kladeny vysoké nároky, co se týče nelineárního zkreslení a vyrovnanosti frekvenční charakteristiky. Celá měřící soustava je zobrazena na obrázku 1.

vůči nelineárnímu zkreslení, které se, třebaže v minimální míře, vždy v měřícím řetězci objevuje. Technika využívající SWEEP signálu dokáže z naměřené impulsní odezvy oddělit nelineární zkreslení vyššími harmonickými složkami [7]. Pokud definujeme odstup signálu od šumu jako poměr průměrného výkonu signálu zaznamenaného měřícím mikrofonem a průměrného výkonu šumu a zkreslení obsaženého ve vypočtené impulsní odezvě (zejména v části odpovídající mnohonásobným odrazům [12]), zjistíme, že použitím techniky se SWEEP signálem dosáhneme až o 20 dB lepších výsledků než u MLS [11]. Výhodou MLS oproti SWEEP je pak odolnost vůči ruchům zaznamenaným během měření, např. v koncertním sále s obecenstvem. Odstup signálu od šumu lze také ovlivnit volbou délky budicího signálu (zdvojnásobením délky MLS sekvence dosáhneme zlepšení o 3 dB) nebo zprůměrováním výsledků vícenásobného měření [11]. Pro výpočet dekonvoluce se používá algoritmus rychlé dekonvoluce, který se skládá z následujících kroků: 1. prodloužení vstupního signálu nulovými vzorky, 2. výpočet spektra vstupního signálu pomocí FFT, 3. prodloužení výstupního signálu nulovými vzorky, 4. výpočet spektra výstupního signálu pomocí FFT, 5. podělení spekter, 6. výpočet impulsní odezvy pomocí zpětné FFT.

Obrázek 1: Blokové schéma měřící soustavy.

Při měření kmitočtové charakteristiky hledáme funkci H(ω), kde ω je spojitý úhlový kmitočet. Podle [6] platí Y (ω) = X(ω)H(ω),

(1)

kde X(ω) je spektrum budicího signálu a Y (ω) je spektrum změřeného signálu. Operací dekonvoluce ve spektrální oblasti danou vztahem H(ω) =

Y (ω) X(ω)

(2)

je možné vypočítat kmitočtovou charakteristiku H(ω). Aby byla funkce H(ω) definovaná pro všechny kmitočty, je nutné, aby budicí signál také obsahoval všechny kmitočty. Vzhledem k tomu, že je navrhován číslicový filtr, stačí, aby budicí signál obsahoval všechny kmitočty v intervalu < 0, fvz /2 >, kde fvz je vzorkovací kmitočet. Této podmínce vyhovuje hned několik signálů, např. pseudonáhodný signál MLS (maximum length sequence) nebo logaritmicky přelaďovaný signál SWEEP [8]. Oba typy signálů mají své výhody a nevýhody. Největší výhodou signálu SWEEP oproti MLS je vysoká míra odolnosti

Tímto postupem lze vypočítat signál s dvojnásobnou délkou než vstupní resp. výstupní signál, který obsahuje impulsní odezvu doplněnou nulovými vzorky. Naměřená a vypočítaná impulsní odezva kytarového komba ENGL je zobrazena na obrázku 2 Jsou zde patrné parazitní impulsní odezvy, které vznikly vlivem nelineárního zkreslení měřícího řetězce a reproduktoru. Jak je vidět, jsou zcela odděleny od hledané impulsní odezvy (leží před ní). Z těchto parazitních impulsních odezev je např. možné určit Volterova jádra [7] a pomocí nich pak simulovat nelineární zkreslení reproduktoru, pokud bychom vycházeli z předpokladu, že nelineární zkreslení budicího zesilovače bylo nulové. Pro návrh lineárního číslicového filtru je však nutné tyto parazitní impulsní odezvy odstranit a samozřejmě i zkrátit délku výsledné impulsní odezvy tak, aby ji bylo možné použít jako koeficienty FIR filtru. Tuto operaci lze provést váhováním oknem, nejčastěji pravoúhlým. Délku okna je nutno volit s ohledem na kmitočtové rozlišení u nízkých kmitočtů, nejlépe pomocí poslechových testů. Pokud není z nějakého důvodu dostupný budicí signál obsahující všechny kmitočty z intervalu < 0, fvz /2 >, je možné použít i techniku s postupným manuálním, či automatizovaným přelaďováním harmonického signálu s určitým krokem. Měření pak probíhá postupným přehráním několika harmonických signálů a současným záznamem jejich odezev. Následně je možné pomocí efektivní hodnoty zaznamenaných signálů rekonstruovat předpokládaný průběh kmitočtové charakteristiky. Tento způsob lze použít i pro více nelineární systémy, kdy by se frekvenční cha-

25 – 2

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

−4

1

kde n je časový index. Na obrázku 4 je provedeno srovnání mezi naměřeným a interpolovaným průběhem impulsní charakteristiky a na obrázku 5 je zobrazeno srovnání mezi kmitočtovými charakteristikami. Délka impulsní charakteristiky je 4096 vzorků. Impulsní charakteristika byla měřena středové pozici s krokem 2cm směrem k okraji reproduktoru.

x 10

h (t)

0.5

0

−0.5

−1 0

0.5

1

1.5

2

2.5

t [s]

0

Obrázek 2: Změřená impulsní charakteristika. V levé části jsou patrné parazitní impulsní charakteristiky vyšších řádů.

rakteristika určovala na základě hodnot základních harmonických složek signálů. Během měření se obvykle přelaďuje kmitočet harmonického signálu logaritmicky, neboť to lépe odpovídá vlastnostem lidského sluchu [14]. Pro výpočet impulsní odezvy je však nutné zvolit lineární krok přelaďování, aby bylo možné např. využít zpětné FFT. Pro návrh koeficientů číslicového filtru typu FIR pak lze využít metodu vzorkování kmitočtové charakteristiky.

M (f) [dB]

−40 −60 −80 20000 −100 5

4

3

2

1

Pozice

20

f [Hz]

Obrázek 3: Modulová frekvenční charakteristika naměřená v sousedních pozicích.

Simulace směrových charakteristik lineární interpolací v časové oblasti 0.02 h (n)

Aby bylo možné nasimulovat směrové frekvenční charakteristiky elektroakustického měniče, je nutné nejprve proměřit kmitočtovou charakteristiku měniče v několika bodech. Na obrázku 3 je zobrazena modulová kmitočtová charakteristika kytarového komba ENGL měřená v několika bodech vedle sebe. K měření byla použita metoda s využitím signálu SWEEP. Zobrazené kmitočtové charakteristiky se nacházejí ve vzdálenosti 2 cm od sebe od okraje reproduktoru směrem ke středu. Z obrázku je patrné, že naměřené kmitočtové charakteristiky mají na nízkých kmitočtech velmi podobný průběh. U vysokých kmitočtů se kmitočtové charakteristiky liší, ale zároveň lidský sluch má u vyšších kmitočtů horší rozlišovací schopnost. Proto je teoreticky možné rekonstruovat frekvenční charakteristiku i mezi měřenými body. Nejjednodušší a zároveň nejrychlejší metodou interpolace je lineární interpolace daná vztahem

0 −0.02 −0.04 0

0.02

0.04

0.06

0.08

0.1

0.06

0.08

0.1

n [−]

−4

5 chyba

3

−20

x 10

0 −5 −10 0

0.02

0.04 n [−]

y = xi (1 − p) + xi+1 p

p ∈< 0, 1 >,

(3)

kde p je vzdálenost interpolovaného bodu od měřených pozic. Chceme sice interpolovat kmitočtové charakteristiky, k dispozici jsou však naměřená data v podobě impulsních odezev a i filtrace vyžaduje koeficienty v podobě vzorků impulsní odezvy. Vzhledem k linearitě Fourierovy transformace mohou být interpolovány přímo vzorky impulsní odezvy h(n) podle vztahu hint (n) = hi (n)(1 − p) + hi+1 (n)p

n ∈< 0, N >,

(4)

Obrázek 4: Průběh změřené a interpolované impulsní charakteristiky (nahoře) a jejich rozdíl (dole). Z obrázku 4 a 5 je patrné, že rozdíl mezi interpolovanou a změřenou kmitočtovou charakteristikou je velmi malý, obzvláště v okolí nízkých kmitočtů. V oblasti vyšších kmitočtů chyba narůstá, ale zároveň také klesá rozlišovací

25 – 3

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

M(f) [dB]

0 −50

−100

2

10

3

10 f [Hz]

4

10

error [dB]

2 0 −2

2

10

3

10 f [Hz]

4

10

Obrázek 5: Modulová frekvenční charakteristika změřená a interpolovaná (nahoře) a jejich rozdíl (dole).

Obrázek 6: Trojrozměrné pole bodů, ve kterých jsou naměřeny impulsní odezvy a nalezení osmi vrcholů krychle potřebných pro výpočet interpolace.

schopnost lidského sluchu [14]. 3.1

Základní model simulace snímání kytarového reproduktorového boxu

Vzhledem k tomu, že akustické pole před elektroakustickým měničem je trojrozměrné, bude i měření impulsních odezev probíhat ve 3 souřadnicích. Tak získáme trojrozměrné pole bodů s konstantními vzdálenostmi ve tvaru krychle, ve kterých máme naměřené impulsní odezvy. Pro libovolný bod uvnitř tohoto pole tedy můžeme nalézt 8 bodů, které tvoří vrcholy menší krychle, ve které bod leží. Trojrozměrné pole bodů a nalezení vrcholů krychle (zvýrazněny červenou barvou) sousedících s hledaným bodem (zelená barva) je naznačeno na obrázku 6. Pro přehlednost je zobrazeno pole o velikosti pouze 3 × 3 × 3 pozic, v reálném případě je samozřejmě žádoucí mít mnohem více naměřených pozic. Pomocí trilineární interpolace pak z osmi známých impulsních odezev vypočítáme hledanou impulsní odezvu konkrétního bodu. Takto tedy získáme impulsní odezvu pomocí které můžeme přesně nasimulovat pozici „ideálníhoÿ mikrofonu před reproboxem. Ideální mikrofon v tomto případě představuje elektroakustický měnič, u kterého naprosto zanedbáme jeho frekvenční a směrovou charakteristiku, konečné rozměry měniče a také případný vliv vzdálenosti od snímaného zdroje, tzv. proximity efekt [13], který se vyskytuje u měničů prvního řádu. . Rozšíření modelu o simulaci vlivu použitého mikrofonu se přímo nabízí. Pro naše účely by plně postačilo, pokud bychom považovali frekvenční charakteristiku mikrofonu závislou pouze na vzdálenosti od zdroje zvuku. Pak by stačilo změřit frekvenční charakteristiky mikrofonu, potažmo tedy získat impulsní odezvy, v několika pozicích ležících na ose kolmé k ploše reproduktoru. Měřící soustava tak, jak je popsaná v kapitole 2, musí být samozřejmě modifikována ve smyslu použití měrného (kalibrovaného) elektroakustického měniče a měřeného mikrofonu. K výpočtu výsledné impulsní

odezvy popisující vliv mikrofonu by pak stačilo využít lineární interpolace popsané vztahem (4). Zanedbání vlivu směrových charakteristik mikrofonu by však bylo v rozporu s praktickými zkušenostmi se snímáním kytarových boxů. Proto je vhodné doplnit model simulace o parametr určující natočení mikrofonu kolem svislé osy procházející co nejblíže mikrofonní kapsle. Tak můžeme nastavit odklon mikrofonu od osy kolmé k ploše reproduktoru. Měření impulsních odezev pak probíhá tak, že v každé z výše popsaných pozic naměříme impulsní odezvy pro několik úhlů natočení mikrofonu. Počet a úhly natočení musí být samozřejmě pro všechny pozice stejné. Z naměřeného dvojrozměrného pole impulsních odezev pak získáme pro libovolnou vzdálenost a úhel natočení mikrofonu potřebnou impulsní odezvu pomocí bilineární interpolace. Tímto postupem se dostaneme do situace, kdy máme impulsní odezvu hspk (n) popisující vliv pozice mikrofonu vzhledem k reproduktorového boxu a impulsní odezvu hmic (n), která popisuje vlastnosti mikrofonu při dané vzdálenosti od reproduktoru a úhlu natočení. V podstatě se tak jedná o sériové spojení dvou dílčích lineárních časově invariantních systémů a proto můžeme jejich impulsní odezvy sloučit do výsledné h(n) dle vztahu h (n) = hspk (n) ∗ hmic (n) ,

(5)

kde ∗ značí operaci lineární konvoluce. Získanou impulsní odezvu h(n) pak použijeme při filtraci pomocí algoritmu rychlé konvoluce. Popsaný základní model simulace má tedy celkem 6 parametrů: souřadnice x, y, z a úhel α popisující polohu mikrofonu vůči reproduktoru a za poslední dva parametry můžeme považovat konkrétní model reproduktorového boxu a mikrofonu, pro které máme naměřené potřebné impulsní odezvy. Lineární interpolace v časové oblasti byla ověřena při práci se sadou impulsních odezev naměřenou celkem v

25 – 4

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

5×5×5 pozicích před reproduktorovým boxem, frekvenční vlastnosti mikrofonu byly naměřeny pro 5 různých úhlů natočení mikrofonu pro každou z pěti vzdáleností. Nejlepších výsledků tato interpolace dosahuje, pokud nastavená vzdálenost přímo odpovídá některé z „vrstevÿ bodů, ve kterých jsme impulsní odezvy naměřili. V tom případě dochází k interpolaci pouze v horizontálním a vertikálním směru a rozdíly mezi interpolovanými a změřenými charakteristikami jsou srovnatelné s rozdílem znázorněným na obrázku 5. Problém ale nastává, pokud se interpolují impulsní odezvy z různých vzdáleností od reproduktorového boxu. S rostoucí vzdáleností totiž dochází ke zpoždění signálu vlivem konečné rychlosti šíření zvuku. Pokud jsou pak interpolovány dvě impulsní odezvy, které mají podobný tvar, ale jsou vzájemně posunuty, dochází ke vzniku efektu hřebenového filtru, který je nejpatrnější v poloze uprostřed mezi interpolovanými polohami. Pokud jsou změřené polohy vzdáleny od sebe 2 cm, tak dochází ke zpoždění 5, 9 × 10−5 s a první minimum hřebenového filtru má polohu 8,5 kHz. Při změně vzdálenosti snímání reproduktoru v reálném čase jsou pak slyšet rušivé artefakty spojené s výskytem právě hřebenového filtru v mezipolohách. Řešením je kompenzace zpoždění nebo lépe využití interpolace v kmitočtové oblasti.

Simulace směrových charakteristik lineární interpolací v kmitočtové oblasti

0.5 y(n)

Interpolace v časové oblasti je výpočetně velmi jednoduchá, její výsledek je ale vždy závislý na konkrétních vstupních impulsních odezvách. Za předpokladu, že hodláme simulaci pomocí impulsní odezvy použít pouze jako kmitočtový filtr, tj. ne například pro simulaci dozvuku, můžeme se soustředit hlavně na kmitočtové vlastnosti výsledného FIR filtru. Tuto podmínku splníme jednak tím, že měření impulsních odezev provádíme v bezodrazovém prostředí, a také tím, že volíme relativně malou délku impulsní odezvy. Pokud také připustíme, že dle Ohmova akustického zákona lidské ucho nevnímá fázové vztahy mezi kmitočtovými složkami signálu [14], omezíme svůj zájem pouze na modul frekvenční charakteristiky. Ke zpětnému získání impulsní odezvy pak můžeme využít např. metody vzorkování kmitočtové charakteristiky. Prvním krokem interpolace je tedy převedení impulsních odezev do kmitočtové oblasti a zanedbání fáze následujícím postupem Mi (ω) = |FFT (hi (n))| . (6)

1

0

−0.5

0

500

1000 n [−]

1500

2000

1 Hammingovo okno 0.5 h(n)

4

transformaci potřebuje i informaci o fázi. K tomu lze využít fáze z jednoho naměřeného bodu, se kterým interpolujeme. Pomocí zpětné Fourierovy transformace interpolovaného spektra Y (ω) získáme impulsní odezvu y(n), kterou je však nutno dále upravit. Jak je vidět na obrázku 7 nahoře, její délka odpovídá použité velikosti okna Fourierovy transformace, obvykle tedy dvojnásobku délky původních impulsních odezev. Hlavním problémem je ale fakt, že nejmarkantnější impulsy leží (symetricky) na začátku a na konci impulsní odezvy. Při dostatečně dlouhé impulsní odezvě by tedy snadno mohlo dojít ke vzniku nežádoucí ozvěny, echa. Abychom získali použitelné „ jádro filtruÿ [10], musíme provést cyklické posunutí vzorků doprava o čtvrtinu délky okna FFT. Poslední úpravou je vynásobení Hammingovým nebo Blackmanovým oknem, jehož střed leží na polovině původní délky impulsní odezvy, tedy nad nejmarkantnějším vzorkem posunuté odezvy, viz. obrázek 7 dole. Vzorky ležící mimo okno vynulujeme a můžeme je tedy zanedbat. Tak můžeme volbou délky použitého okna měnit délku výsledné hledané impulsní odezvy h(n), ovšem při očekávatelné ztrátě přesnosti, která se projeví zejména na nižších kmitočtech.

0

−0.5

0

500

1000 n [−]

1500

2000

Obrázek 7: Průbeh vypočtené impulsní odezvy před (nahoře) a po cyklickém posunutí a vynásobení oknem (dole).

Tím navíc přejdeme od komplexních čísel zpět k reálným hodnotám a další výpočty se tak zjednoduší a zrychlí. Samotnou interpolaci pak provedeme dle vztahu Y (ω) = (1 − p)Mi (ω) + pMi+1 (ω) p ∈< 0, 1 >,

(7)

kde p je vzdálenost interpolovaného bodu od naměřených pozic. Ke zpětnému získání impulsní odezvy pak můžeme použít metodu vzorkování kmitočtové charakteristiky, resp. její modifikovanou obdobu [10]. Pro zpětnou

Srovnání interpolace v časové a kmitočtové oblasti je zobrazeno na obrázku 8. V průběhu kmitočtové charakteristiky interpolované v časové oblasti je jasně patrný vliv hřebenového filtru, ke kterému dochází v důsledku zpoždění impulsních odezev v různé vzdálenosti. Naopak, interpolace v kmitočtové oblasti má očekávaný průběh.

25 – 5

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

vlastností mikrofonu, bude mít interpolace tvar

15

M(f) [dB]

10 5 0 −5 −10 −15 −20

pozice 1 pozice 2 interpolace v þasové oblasti interpolace ve frekvenþní oblasti 2

3

10

10

4

10

f [Hz]

Obrázek 8: Srovnání interpolace v časové a kmitočtové oblasti.

4.1

Rozšířený model simulace snímání kytarového reproduktorového boxu

Dříve zmíněných 6 parametrů základního modelu obsahuje i model rozšířený. Na základě dané pozice a úhlu natočení mikrofonu tedy obdobným způsobem, avšak pomocí trilineární, resp. bilineární interpolace v kmitočtové oblasti získáme hodnoty modulu kmitočtové charakteristiky Yspk (ω), resp. Ymic (ω). Výsledné hodnoty modulu frekvenční charakteristiky Y (ω), které použijeme pro zpětné získání impulsní odezvy, pak získáme pouhým vynásobením dle vztahu Y (ω) = Ymic (ω)Yspk (ω). (8)

YmicR (ω) = kYmic (ω) + (1 − k)Mflat (ω) k ∈< 0, 1 >, (10) kde k je parametr řídící míru simulace. Dosazením (10) do vztahu (8) a použitím stejného postupu při výpočtu Mkomp (ω) můžeme libovolně ovlivňovat míru simulace i kompenzace každého objektu. Popsané řízení míry simulace by z důvodu linearity Fourierovy transformace bylo teoreticky možno provést i v časové oblasti. Tam by stačilo vypočtenou impulsní odezvu interpolovat k hflat (n) definovanou jako  1 pro n=1 hflat (n) = (11) 0 pro 1 < n ≤ N, kde N je délka impulsní odezvy. Nastává zde ale problém obdobný tomu, když se interpolují impulsní odezvy naměřené v různých vzdálenostech od reproduktorového boxu. Poloha maxima v hflat (n) je totiž zpravidla odlišná od polohy maximální hodnoty v obecné impulsní odezvě a tak se při interpolaci projeví vliv hřebenového filtru. Srovnání řízení míry simulace v časové a kmitočtové oblasti je zobrazeno na obrázku 9. 20 M(f) [dB]

20

0 k=1

−20 1

K získání Mkomp (ω) je samozřejmě možno použít výše popsaný způsob interpolace v kmitočtové oblasti a tak vykompenzovat libovolnou konfiguraci reproduktorového boxu a mikrofonu, pro které máme naměřeny potřebné impulsní odezvy. Druhým zlepšením, které rozšířený model přináší, je možnost řízení míry simulace, resp. kompenzace každého objektu ze simulovaného, resp. kompenzovaného řetězce. Toho dosáhneme lineární interpolací hodnot modulu kmitočtové charakteristiky daného objektu s hodnotami Mflat (ω), které odpovídají zesílení kmitočtových složek rovno 0 dB. Numericky jsou tedy všechny hodnoty modulu kmitočtové charakteristiky Mflat (ω) rovny jedné. Chceme-li tedy např. ovlivnit míru simulace frekvenčních

k=0

3

10

4

10

10

f [Hz]

Fakt, že interpolace probíhá v kmitočtové oblasti, nám přináší další možnosti, o kterých se v časové oblasti dá jen těžko uvažovat. První z nich je možnost kompenzace nechtěných kmitočtových vlastností signálu, které byly zapříčiněny např. mikrofonem použitým pro pořízení nahrávky. Máme-li změřenou jeho impulsní odezvu, dle vztahu 6 získáme modul frekvenční charakteristiky Mkomp (ω). Kompenzaci, tedy dekonvoluci, pak provedeme pouhým podělením výsledku interpolace dle vztahu (9)

k = 0.3

2

10

Ykomp (ω) = Y (ω)/Mkomp (ω).

k = 0.6

M(f) [dB]

20 0 k=1

−20 1

10

k = 0.6

k = 0.3

2

3

10

10

k=0 4

10

f [Hz]

Obrázek 9: Řízení míry simulace v časové (nahoře) a kmitočtové (dole) oblasti.

Další užitečnou možností, kterou nám interpolace v kmitočtové oblasti přináší, je přímá kontrola zesílení kmitočtových složek. Tak můžeme například zabránit příliš velikému zesílení v určitém pásmu kmitočtů, které by mohlo být nebezpečné pro elektroakustické měniče nebo dokonce pro lidský sluch. Taková situace může lehce nastat, pokud kompenzujeme kmitočtovou charakteristiku objektu, ve které se vyskytují výrazná minima. Rozšířený model simulace tedy může mít variabilní počet parametrů, který je závislý pouze na tom, co vše od něj očekáváme. K šestici základních parametrů popisujících konfiguraci simulovaných objektů může přibýt až šest para-

25 – 6

2011/25 – 13. 6. 2011

VOL.13, NO.3, JUNE 2011

metrů popisujících objekty kompenzované. Budeme-li řídit míru simulace a kompenzace každého objektu zvlášť, jsou to další čtyři parametry. Pokud budeme pracovat s kompenzací, je rozumné zavést poslední parametr, kterým omezíme maximální výsledné zesílení dílčích frekvenčních složek. Blokové schéma rozšířeného modelu simulace snímání kytarového reproduktorového boxu je zobrazeno na obrázku 10.

Pro výpočet konvoluce existuje celá řada algoritmů, při zpracování zvukových signálů v reálném čase je však požadována velmi nízké dopravní zpoždění. Zároveň je však třeba udržet nízkou výpočetní náročnost celého algoritmu i pro případ, kdy pracujeme s velice dlouhou impulsní odezvou. Tyto rozporné požadavky dokáže splnit algoritmus konvoluce s dvojitou FDL (frequency delay line) [4]. Tento algoritmus počítá konvoluci v kmitočtové oblasti tím způsobem, že vstupní impulsní odezvu ještě před převodem do kmitočtové oblasti dělí na bloky dvou různých délek. Krátké bloky na začátku umožňují splnit požadavek na nízké dopravní zpoždění, následující delší bloky pak snižují výpočetní náročnost algoritmu. To společně s redukcí počtu potřebných zpětných Fourierových transformací během výpočtu realizací metody přičtení přesahu již v kmitočtové oblasti vede k velmi nízké výpočetní náročnosti tohoto algoritmu. Velikost bloků druhé FDL je možné nastavit dle délky impulsní odezvy. To činí z tohoto algoritmu univerzální prostředek pro výpočet konvoluce v reálném čase pro délky impulsní odezvy od pár stovek vzorků až po hodnoty odpovídající několika desítkám sekund. Podrobnosti o implementaci a analýzu výpočetní náročnosti algoritmu lze nalézt např. v [12].

6

Obrázek 10: Blokové schéma rozšířeného modelu simulace snímání kytarového reproduktorového boxu.

5

Implementace

Nejsložitější částí algoritmu, který bude simulovat reproduktorový box v reálném čase, je interpolace mezi změřenými daty, a to v časové nebo frekvenční oblasti. Pro samotnou interpolaci uvnitř algoritmu je třeba mít k dispozici naměřená data, která mohou být např. uložena v binárním souboru, ze kterého budou čtena. Nad těmito daty se pak provádí interpolace. Při použití interpolace v časové oblasti pracujeme přímo s naměřenými impulsními odezvami a jako výsledek interpolace obdržíme také přímo žádanou impulsní odezvu. V případě interpolace v kmitočtové oblasti jsou z uložených impulsních odezev nejprve pomocí vztahu (6) vypočítány hodnoty modulu kmitočtové charakteristiky. Po té následují výpočty využívající interpolace dle vztahu (7) závislé na tom, co vše má model dělat. Nakonec jsou pomocí výše popsaného postupu vypočteny vzorky výsledné impulsní odezvy. Ty jsou pak použity jako koeficienty FIR filtru a odezva na vstupní signál je vypočítána pomoci konvoluce.

Závěr

Simulace kmitočtových vlastností elektroakustických měničů, a to zvláště v blízkém poli, je poměrně složitý úkol. Avšak lze jej řešit pomocí kombinace měření kmitočtových charakteristik v různých pozicích a jejich vzájemnou interpolací. Tímto způsobem je pak možné interpolovat kmitočtovou charakteristiku v celé oblasti před měřeným elektroakustickým měničem. V tomto článku jsou prezentovány dvě verze interpolace. Interpolace v časové oblasti má obrovskou výhodu ve své jednoduchosti a rychlosti finální implementace. Ovšem objevují se zde problémy, pokud dochází k interpolaci impulsních odezev, které byly naměřeny v různé vzdálenosti od měřeného měniče. Vlivem časového zpoždění dochází při interpolaci k výskytu hřebenového filtru, který výslednou kmitočtovou charakteristiku zkresluje. Tuto nevýhodu odstraňuje interpolace v kmitočtové oblasti, která je však výpočetně náročnější. Na druhou stranu však nabízí mnoho dalších rozšiřujících možností. Výsledkem obou typů interpolace je impulsní odezva, jejíž vzorky se dají použít jako koeficienty FIR filtru. Tímto způsobem je pak možné poměrně jednoduše implementovat algoritmus, který bude simulovat daný elektroakustický měnič v reálném čase. Ideálním řešením pro implementaci výpočtu konvoluce je využití algoritmu dvojité FDL, který má vhodné vlastnosti pro zpracování signálů v reálném čase. Další práce bude směřována k ověření simulace směrových charakteristik pomocí poslechových testů a stanovení velikosti kroku mřížky, ve které je nutné elektroakustický měnič měřit.

25 – 7

2011/25 – 13. 6. 2011

7

VOL.13, NO.3, JUNE 2011

Poděkování

Tento projekt byl realizován za finanční podpory z prostředků státního rozpočtu prostřednictvím Ministerstva průmyslu a obchodu v rámci projektu ev. č. FR-TI1/495.

[12] TICHÝ, V.: Digitální zvukový efekt typu reverb využívající konvoluci signálu s impulsní charakteristikou poslechového prostoru. Diplomová práce, Brno University of Technology, Brno, 2009. [13] VLACHÝ, V.: Praxe zvukové techniky. Praha: Nakladatelství Muzikus, 2000, ISBN ISBN 80-86253-05-8, 257 s.

Literatura [1] BORWICK, J.: Looudspeaker and Headphone Handbook. Oxford: Reed Educational and Professional Publishing Ltd, třetí vydání, 2001, ISBN 0-240-51578-1, 718 s.

[14] ZWICKER, A.; FASTL, H.: Psychoacoustics. New York: Springer, druhé vydání, 1999, ISBN 3-54065063-6, 416 s.

[2] FRIEDLANDER, B.; PORAT, B.: The Modified Yule-Walker Method of ARMA Spectral Estimation. In IEEE Transactions on Aerospace Electronic Systems, ročník AES-20, 1984, s. 158–173.

[15] YEH, T. D.; BANK, B.; KARJALAINEN, M.: Nonlinear Modelling of a Guitar Loudspeaker Cabinet. In Proc. Digital Audio Effects (DAFx-07), Bordeaux, France, Sep. 1-4, 2008, s. 89–96.

[3] FUNG, C. C.; KOK, C.: Mixed-Domain Reduced Order IIR Approximation for FIR Filter. In Proc. of the XII. European Signal Processing Conference EUSIPCO 2004, Vienna, Austria, Sep. 6-10, 2004. [4] GARCIA, G.: Optimal Filter Partition for Efficient Convolution with Short Input/Output Delay. In 113th AES Convention Paper 5660, Los Angeles, US, 2002, s. 1–9. [5] ŠKVOR, Z.: Akustika a elektroakustika. Praha: Academia, 2001, ISBN ISBN 80-200-0461-0. [6] MANOLAKIS, D. G.; PROAKIS, J. G.: Digital Signal Processing: Principles, Algorithms and Applications. Englewood Cliffs, NJ, USA: Prentice Hall, 3rd vydání, 1991. [7] NOVÁK, A.; LOTTON, S.; KADLEC, F.: Modeling of Nonlinear Audio Systems Using Swept-Sine Signals: Application to Audio Effects. In Proceedings of the 12th International Conference on Digital Audio Effects DAFx09, Como, Italy, Sep. 1-4, 2009, s. 1–6. [8] PICINALLI, L.: Techniques for the Extraction of the Impulse Response of a Linear and Time-Invariant System. In Proc. of the DRMN-06: DMRN Doctoral Research Conference 2006., London, Great Britain, 2006. [9] SMETANA, C.: Praktická elektroakustika. Praha: SNTL, 1981. [10] SMITH, S. W.: The Scientist & Engineer’s Guide to Digital Signal Processing. California Technical Pub., první vydání, 1997, ISBN 0966017633, 626 s. [11] Stan, G. B. S.; EMBRECHTS, J. J.; ARCHAMBEAU, D.: Comparison of Different Impulse Response Measurement Techniques. J. Audio Eng. Soc, ročník 50, č. 4, 2002: s. 249–262. URL http://www.aes.org/e-lib/browse.cfm? elib=11083 25 – 8