MODELLO STRUTTURALE DI ORECCHIO ESTERNO PER IL...

Universita di Padova

Dipartimento di Ingegneriadell’Informazione

MODELLO STRUTTURALE

DI ORECCHIO ESTERNO

PER IL RENDERING

AUDIO 3D

Laureando:

Michele Geronazzo

Relatore:

Dott. Federico Avanzini

A. A. 2008 - 2009

ii

Ai miei primi supporter e ammiratori, mamma Rosanna e papa Ermanno,

a nonna Lena,

alla mia Gioia quotidiana Cri :-),

a Diego, amico per sempre nella Luce che ci abbraccia.

Abstract

La Head Related Transfer Function (HRTF) caratterizza i segnali acustici

creati dalla diffusione del suono quando raggiunge il corpo di una persona.

Le HRTF misurate sperimentalmente contengono informazioni su vari feno-

meni: le riflessioni causate dalle parti del corpo (busto, spalle e ginocchia),

la diffrazione della testa e l’effetto di riflessione/difrazzione causato dal pa-

diglione auricolare (termine inglese: pinna). I modelli strutturali (Algazi

et al., 2002; Brown e Duda, 1998) cercano di stabilire una relazione diretta

tra le caratteristiche della HRTF e l’anatomia umana; la semplicita che li

caratterizza permette un’implementazione efficiente per il rendering audio in

applicazioni real-time.

L’obiettivo di questa tesi e analizzare il contributo della pinna all’interno

della HRTF personale, propria del singolo individuo, e modellarne le carat-

teristiche che hanno interesse percettivo per la localizzazione verticale del

suono. L’ applicazione di metodi per l’individuazione delle risonanze e dei

notch presenti nello spettro, assieme ad un partial tracking nello spazio verti-

cale, contribuiscono alla stesura di un modello strutturale per la pinna. Per il

lavoro di analisi sono stati utilizzati i dati del CIPIC HRIR database (Algazi

et al.,2001).

Parole chiave: audio 3D, spatial hearing, head-related transfer function,

pinna, modello strutturale.

iii

iv ABSTRACT

Indice

Abstract iii

Introduzione 1

1 Background e Stato dell’Arte 5

1.1 3D Audio System . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Localizzazione del suono: l’audio binaurale . . . . . . . . . . . 11

1.2.1 Sistema di riferimento . . . . . . . . . . . . . . . . . . 11

1.2.2 HRIR e HRTF . . . . . . . . . . . . . . . . . . . . . . 12

1.2.3 La Testa . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.2.4 L’Orecchio . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.2.5 Busto e Spalle . . . . . . . . . . . . . . . . . . . . . . . 20

2 Modellazione HRIR e HRTF 23

2.1 HRIR e HRTF personali . . . . . . . . . . . . . . . . . . . . . 24

2.1.1 Misurazione acustica . . . . . . . . . . . . . . . . . . . 25

2.1.2 Interpolazione su database . . . . . . . . . . . . . . . . 25

2.1.3 Metodi Numerici . . . . . . . . . . . . . . . . . . . . . 28

2.1.4 Modelli Poli-Zeri . . . . . . . . . . . . . . . . . . . . . 29

2.1.5 Espansioni in serie . . . . . . . . . . . . . . . . . . . . 32

2.2 Modelli Strutturali . . . . . . . . . . . . . . . . . . . . . . . . 37

2.3 PRTF e antropometria . . . . . . . . . . . . . . . . . . . . . . 42

2.3.1 Modelli Strutturali per PRTF . . . . . . . . . . . . . . 44

v

vi INDICE

3 Metodologie e strumenti 55

3.1 Cepstrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.2 Peak Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.3 Modelli Sinusoidali . . . . . . . . . . . . . . . . . . . . . . . . 61

3.3.1 Peak Detection . . . . . . . . . . . . . . . . . . . . . . 63

3.3.2 Partial Tracking . . . . . . . . . . . . . . . . . . . . . . 64

3.4 CIPIC Database . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.4.1 Misurazioni . . . . . . . . . . . . . . . . . . . . . . . . 68

3.4.2 Antropometria . . . . . . . . . . . . . . . . . . . . . . 69

4 Analisi della PRTF 73

4.1 Algoritmo di analisi . . . . . . . . . . . . . . . . . . . . . . . . 74

4.1.1 Inizializzazione . . . . . . . . . . . . . . . . . . . . . . 75

4.1.2 Calcolo del residuo . . . . . . . . . . . . . . . . . . . . 77

4.1.3 Ricerca parametri per il filtro multi-notch . . . . . . . 78

4.1.4 Costruzione del filtro multi-notch . . . . . . . . . . . . 83

4.1.5 Risultato finale . . . . . . . . . . . . . . . . . . . . . . 84

4.2 Analisi dei dati . . . . . . . . . . . . . . . . . . . . . . . . . . 88

4.2.1 Analisi dei notch . . . . . . . . . . . . . . . . . . . . . 89

4.2.2 Analisi delle risonanze . . . . . . . . . . . . . . . . . . 90

5 Un Modello Strutturale per la pinna 95

5.1 Descrizione e motivazione . . . . . . . . . . . . . . . . . . . . 95

5.2 Filtro sintetico . . . . . . . . . . . . . . . . . . . . . . . . . . 96

5.2.1 Filtro sintetico per le risonanze . . . . . . . . . . . . . 97

5.2.2 Filtro sintetico per i notch . . . . . . . . . . . . . . . . 100

5.2.3 Ricostruzione . . . . . . . . . . . . . . . . . . . . . . . 103

Conclusioni e Sviluppi futuri 109

A Algoritmo d’analisi, variazioni parametriche 113

A.1 Parametro Nceps . . . . . . . . . . . . . . . . . . . . . . . . . . 113

A.2 Parametro soglia . . . . . . . . . . . . . . . . . . . . . . . . . 114

INDICE vii

A.3 Parametro rid . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

B Sorgenti MATLAB 119

Bibliografia 149

viii INDICE

Introduzione

Lo scopo di questa tesi e sviluppare algoritmi efficienti ed accurati per la

sintesi di audio binaurale in ambiente real-time, con particolare attenzione

alla localizzazione verticale del suono.

L’idea di audio binaurale trae ispirazione dalla duplex-theory del fisico

Lord Rayleigh, sviluppata piu di cent’anni fa e da allora considerata storica-

mente come punto di partenza sono stati condotti numerosi studi acustici e

percettivi per comprendere i contributi di testa, busto, spalle e orecchio ester-

no, forniti al sistema uditivo binaurale. In questi ultimi decenni la letteratura

si e arricchita di lavori che hanno indagato i segnali di localizzazione spaziale

e successivamente esplorato i fenomeni di riflessione e di diffrazione nel cam-

po vicino (sorgente posta a distanza inferiore di un metro dall’ascoltatore) e

nel campo lontano.

Il continuo aumento della potenza di calcolo ha permesso lo sviluppo di

numerosi sistemi in tempo reale per l’audio (VAD - Virtual Acoustic Display),

rivelando la necessita di incorporare al loro interno dei modelli efficaci ma

allo stesso tempo efficienti. Un modello sferico per la testa, un busto ellissoi-

dale e un orecchio esterno composto da una conca cilindrica sono solo alcu-

ni degli esempi di approssimazioni introdotte per semplificare e soprattutto

controllare gli effetti delle parti anatomiche sull’onda sonora incidente.

Gli studi condotti nel campo dell’audio binaurale hanno alimentano nel-

l’ultimo decennio molteplici applicazioni. Da servizi di teleconferenza a de-

sign di ambienti virtuali immersivi e multimodali; realta virtuale a scopo

intrattenitivo e riabilitativo; infine tecnologie embedded in dispositivi mobili.

1

2 INTRODUZIONE

Il suono emesso da una sorgente giunge ai timpani trasformato profonda-

mente e gli effetti filtranti sono talmente complessi che non tutti si riescono a

sintetizzare in maniera efficace. Alcune problematiche aperte riguardano la

confusione front-back, l’ esternalizzazione del suono e la localizzazione verti-

cale della sorgente sonora. Proprio a quest’ultima e collegato lo studio della

pinna, come strumento principale per questo tipo di discriminazione.

La tesi e concentrata sull’analisi del solo contributo della pinna, come

corpo riflettente e risonante. Sono stati presi degli accorgimenti per limitare

le influenze delle altre componenti anatomiche. I risultati ottenuti sono stati

utilizzati per proporre un primo modello strutturale di pinna, che verra even-

tualmente assemblato in un progetto dallo scopo piu ambizioso: un modello

globale di HRTF (Head Related Transfer Function).

Il lavoro e stato distribuito in tre fasi principali:

1. fase di documentazione e stesura dello stato dell’arte;

2. fase di analisi dei dati;

3. fase di modellazione;

Le informazioni raccolte nella prima fase sono esposte e ordinate nei Capi-

toli 1 e 2. Nel primo capitolo vengono fornite delle nozioni di background

e di carattere generale sugli ambienti audio tridimensionali, assieme ad una

esposizione mirata alla localizzazione del suono in sistemi audio binaura-

li. Il secondo capitolo espone tutte le possibili tecniche, finora sviluppate,

per la modellazione della HRTF personali; vengono approfonditi i modelli

strutturali per la HRTF fino a giungere alla modellazione del solo contribu-

to della pinna (PRTF). Dopo questi capitoli iniziali che contestualizzano il

problema, viene intrapresa la fase di analisi. Il capitolo 3 fa da introduzione,

richiamando risultati teorici e descrivendo gli strumenti utilizzati. Il quar-

to capitolo propone e descrive nel dettaglio un algoritmo di analisi dei dati

in grado di separare i contributi delle risonanze e delle riflessioni. Ognuno

dei due fenomeni verra relazionato a studi autorevoli gia effettuati, e posto

come punto d’inizio per la progettazione di un modello strutturale. Il terzo

INTRODUZIONE 3

capitolo racchiude richiami alle tecniche di DSP utilizzate per la costruzio-

ne dell’algoritmo e del modello (calcolo dell’inviluppo via cepstrum, il peak

tracking, etc.).

Il capitolo finale contiene una prima proposta di filtro sintetico per la

PRTF, dove i due contributi ad opera della pinna (risonanze e riflessioni),

vengono modellati indipendentemente per poi essere combinati assieme.

Le due Appendici forniscono informazioni aggiuntive sul comportamento

dell’algoritmo d’analisi (Appendice A) e mettono a disposizione il codice

MatlabTM sviluppato come supporto all’intero lavoro (Appendice B).

4 INTRODUZIONE

Capitolo 1

Background e Stato dell’Arte

Il suono ricopre molti ruoli nella multimedialita e nel sistema di realta virtua-

le. Viene usato principalmente per la comunicazione (informazioni, linguag-

gio e audio non parlato) o per intrattenimento (musica, ambiente, voce, etc.).

La spazializzazione del suono e il suono 3D hanno una lunga storia di ricerca.

Partendo dalla prima teoria di Lord Rayleigh, chiamata duplex theory, all’ini-

zio del 20esimo secolo (Rayleigh, 1907), il primo concetto di acustica virtuale

viene esplorato nel campo delle teleoperazioni. Negli anni Trenta Spandock

(1934) progetto un modello in scala e cerco di dare spazialita ad una sa-

la da concerto. Durante la Prima Guerra Mondiale, l’individuazione degli

aerei nemici fu migliorata utilizzando dei segnali acustici spazializzati, cosı

da permettere un’immediata localizzazione tridimensionale acustico-visuale,

mediante l’utilizzo di cuffie rudimentali. Queste prime applicazioni hanno lo

stesso scopo dell’acustica virtuale dei giorni d’oggi: il design virtuale di spazi

audiovisivi, un’accurata simulazione e riproduzione tridimensionale di eventi

acustici.

I sistemi per l’audio tridimensionale sono spesso associati con i VAD

(Virtual Acoustic Display). I VAD racchiudono modelli per la sorgente, l’am-

biente e l’ascoltatore in modo da produrre un’esperienza virtuale totalmente

immersiva [39]. I campi delle tele/video-conferenze, dei giochi e della mul-

timedialita sono oggigiorno i principali fruitori delle applicazioni 3D-audio.

5

6 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE

Inoltre con la nascita di standard quali VRML (Virtual Reality Modeling

Language, ISO/IEC, 1997) e MPEG-4 (ISO/IEC, 1999), l’uso dell’acustica

virtuale si e espanso nel contesto delle simulazioni audiovisive. La crescita

delle tecnologie per il posizionamento 3D del suono e stata significativa grazie

all’uso sempre piu frequente di cuffie e casse per una riproduzione stereofoni-

ca ad alta definizione e all’avvento di sistemi multicanale con relativi formati

audio. Questo ha portato allo sviluppo di numerose applicazioni, anche di

uso domestico, quali i sistemi Home Theater.

In questo capito introduttivo tratteremo i concetti basilari dell’audio 3D,

dando una breve panoramica sulle metodologie di rendering per raggiungere

la spazializzazione di una sorgente sonora. Definiremo una terminologia spe-

cialistica che ci aiutera a capirne i concetti fondamentali. Ci soffermeremo

sullo stato dell’arte per i sistemi binaurali, che si appoggiano per la riprodu-

zione sonora alle tradizionali apparecchiature a 2 canali stereofonici (cuffie).

Grazie ai moderni metodi di signal processing e possibile progettare e im-

plementare un ambiente virtuale d’ascolto 3D e molteplice effettistica audio

da poter inserire in prodotti audiovisivi che vengono riprodotti da computer,

televisioni e sistemi hi-fi.

1.1 3D Audio System

Le caratteristiche dell’ascolto umano sono gli elementi che si vogliono model-

lare con la ricerca nel suono 3D e nell’acustica virtuale. La progettazione e

l’implementazione di ambienti virtuali puo essere divisa in tre passi[37]:

a) definizione: La definizione di un ambiente acustico virtuale include la

conoscenza delle caratteristiche del sistema da realizzare, informazioni

sulla sorgente sonora, sulla geometria della stanza e sugli ascoltatori.

Questo studio e portato avanti off-line a priori prima di una possibile

implementazione in simulazioni real-time.

b) modellazione: La parte di modellizzazione e divisa in tre rami che tro-

vano corrispondenza con la fase precedente e rispecchiano il paradigma

1.1. 3D AUDIO SYSTEM 7

Figura 1.1: Concetti per la modellazione di un ambiente virtuale.

source-medium-receiver, ampiamente conosciuto nei sistemi di teleco-

municazioni.

c) riproduzione: Metodi per ricreare l’illusione di un campo sonoro tridi-

mensionali con l’aiuto di trasduttori e tecniche di signal processing.

La parte di modellazione e divisa a sua volta in tre parti:

Modelli della sorgente I modelli della sorgente permettono di produrre un

suono all’interno di una scena audiovisiva. Il metodo piu usato e quello

di impiegare dei suoni pre-registrati, ma in molti sistemi le sorgenti

vengono trattate come puntiformi omnidirezionali. Questa approssi-

mazione pero non sempre e soddisfacente in quanto alcuni strumenti

hanno un pattern di radiazioni dipendente dalla frequenza. Per esem-

pio, una tipica sorgente sonora emana piu energia nell’emisfero frontale;

variando la distanza angolare dall’asse di direzione si avra un filtraggio

passa-basso che ne attenuera le radiazioni.


Modelli del mezzo di trasmissione acustica Questi modelli hanno il com-

pito di descrivere il comportamento della propagazione del suono in uno

spazio acustico (room acoustics). L’obiettivo della simulazione e calco-

lare la curva temporale di energia (ETC) di una stanza e di conseguenza

derivarne gli attributi acustici [38] come il tempo di riverbero (Tr), il

clarity index (Ct), lo strength index (G), etc. I metodi piu usati sono

il ray-based method e l’image-source method, altri metodi computazio-

nalmente piu impegnativi si basano sulle tecniche wave-based come il

metodo agli elementi finiti (FEM), il metodo agli elementi al contor-

no (BEM), e il metodo alle differenze finite nel dominio del tempo

(FDTD). Nelle simulazioni real-time si cerca di trovare un compromes-

so tra risorse computazionali e accuratezza, per questo si modellano

solo il suono diretto e le prime riflessioni, usando l’acustica geometrica,

mentre per le riverberazioni nel campo lontano si usano delle strutture

di filtri digitali ricorsivi.

Modelli del ricevitore-ascoltatore Vengono prese in considerazione le pro-

prieta dell’ascolto umano; si cerca di dare un senso di direzionalita (es.

ILD e ITD v. paragrafo 1.2) e grazie alla head-related transfer

function (HRTF) vengono modellate le riflessioni e filtraggio ad ope-

ra della testa, busto, spalle e orecchio esterno dell’ascoltatore. Grazie

allo sviluppo di tecniche di misurazione precise per la HRTF e proget-

tazione di filtri efficienti e stato possibile utilizzare questo approccio

per il rendering 3D in tempo reale di ambienti sonori virtuali.

Per la parte di riproduzione consideriamo i tre metodi principali:

Riproduzione binaurale Questo metodo utilizza un paio cuffie e la HRTF

misurata viene approssimata dalla progettazione di filtri che catturano

le caratteristiche acustico-percettive piu importanti. In Fig. 1.2a un se-

gnale monofonico nel dominio del tempo xm(n) e filtrato con due HRTF

approssimate Hl(z) e Hr(z) per creare un’immagine stereofonica della

singola sorgente virtuale. Con questo metodo si ha il vantaggio di evi-

1.1. 3D AUDIO SYSTEM 9

Figura 1.2: a) Processo binaurale e b) processo binaurale con cancellazione

dell’effetto crosstalk.

tare l’influenze di disturbo nella simulazione dovute alla stanza in cui

e immerso l’ascoltatore e alla sua posizione. Tuttavia sono necessarie

delle tecniche compensative che tengono conto del movimento della te-

sta nello spazio simulato (head-tracking); inoltre per avere un risultato

efficace nella localizzazione e nella esternalizzazione di un suono virtua-

le la HRTF deve essere personale e le tecnologie di misurazione sono

troppo costose per essere impiegate a livello commerciale. Questo con-

cetto sara piu chiaro con le spiegazioni nei prossimi paragrafi in quanto

tratteremo ampiamente la tematica riguardante l’audio binaurale, e

ci renderemo conto di come la ricostruzione di HRTF personalizzabili

soddisfacenti sia tutt’oggi un problema aperto.

Riproduzione binaurale con crosstalk canceled Questo metodo utiliz-

za degli altoparlanti per la riproduzione binaurale dello spazio acustico

virtuale ed e molto diverso dal metodo precedente. La direzionalita

caratteristica dei suoni binaurali e messa a rischio dal fenomeno di

crosstalk per il quale un suono riprodotto da un altoparlante posizio-

nato per esempio a destra viene sentito da entrambe le orecchie. Il

problema di cancellare il crosstalk e stato formalizzato gia negli anni

’60 e sono state proposte numerose risoluzioni. Nella sintesi binaurale


per una coppia di altoparlanti (v. Fig. 1.2b) i due segnali privi della

componente di crosstalk xl(n) e xr(n) sono riprodotti in modo da gene-

rare gli output desiderati yl(n) e yr(n). Le loudspeaker-to-ear transfer

function dipendenti dalla direzione Hi(z) e Hc(z) (in questo caso con-

sideriamo i contributi destro e sinistro simmetrici) devono ricreare lo

stesso effetto dell’ascolto in cuffia. Possiamo vedere la progettazione

della HRTF separata dal filtro di cancellazione dell’effetto crosstalk,

oppure combinare le due componenti in strutture shuffler. Le due limi-

tazioni principali di questo metodo sono: 1) la criticita della posizione

dell’ascoltatore 2) la criticita delle condizioni della stanza d’ascolto.

Lo spazio sonoro dove l’ascolto risulta migliore(sweet spot) viene cosı

molto limitato, anche se in quella posizione l’ascoltatore puo avere una

percezione molto convincente della simulazione.

Riproduzione multicanale Una soluzione naturale per produrre uno spa-

zio d’ascolto tridimensionale e quella di utilizzare piu altoparlanti po-

sizionandoli in un ambiente e dirigere il segnale che si vuole riprodurre

agendo sul panning tra le casse. Sistemi digitali audio multicanale per

l’home theater e per l’intrattenimento offrono un suono 3D che vie-

ne decodificato da un materiale stereofonico (come Dolby ProLogic) o

usano un decodificatore multicanale come nel caso di Dolby Digital e

del Digital Theater System (DTS). Lo standard ISO/MPEG-2 AAC

offre la possibilita di definire 5.1 canali al rate di 320 kb/s. Il proble-

ma principale che coinvolge i sistemi multicanale con piu di 2 canali e

l’utilizzo e il posizionamento di N altoparlanti. Nonostante lo spazio

d’ascolto sia notevolmente ampliato rispetto ad un ascolto binaurale

con cancellazione del crosstalk, l’ascoltatore resta vincolato alla stan-

za in cui il sistema e stato testato e incontra maggiori complicazioni

nel caso di spostamento delle apparecchiature. Infine un importante

vantaggio di questi sistemi sta nel non utilizzo di alcuna modellazione

dell’ascoltatore e quindi di nessuna HRTF.

1.2. LOCALIZZAZIONE DEL SUONO: L’AUDIO BINAURALE 11

Dopo una panoramica generale sulle tecniche e tecnologie nel campo

dell’audio 3D, focalizziamo la nostra attenzione sui modelli dell’ascoltatore

(HRTF-modeling) per approfondire l’argomento della localizzazione di una

sorgente sonora in favore di una riproduzione binaurale tramite l’uso di cuffie.

1.2 Localizzazione del suono: l’audio binau-

rale

La percezione di una sorgente sonora e trasformata dall’orecchio esterno e

da alcune parti del corpo individuate come la testa, le spalle e il busto. Il

segnale sonoro, cosı trasformato, percorre il canale uditivo fino a percuotere

il timpano. Questi cambiamenti ad opera delle componenti corporee conten-

gono le informazioni spaziali che vengono codificate in attributi temporali e

spettrali. Molte ricerche hanno dimostrato come l’uomo usi alcuni segnali

di localizzazione (localization cues); tre principali sono: la Interaural Time

Difference (ITD), la Interaural Level Difference (ILD) e pinna cues. Varia-

no in funzione dell’azimuth, dell’elevazione e della frequenza e forniscono le

informazioni all’ascoltatore per localizzare un suono nello spazio. In questo

paragrafo tratteremo i segnali di localizzazione, mettendoli in relazione con

le parti del corpo che li generano. Completeremo una panoramica genera-

le sulle caratteristiche percettive fondamentali contenute in un modello di

ascoltatore.

1.2.1 Sistema di riferimento

Definiamo in questo paragrafo la terminologia per descrivere il sistema di

riferimento utilizzato. Il sistema di coordinate e del tipo interaurale-polare1,

quindi le dimensioni che descrivono lo spazio tridimensionale sono:

1per interaurale si intende semplicemente tra le orecchie. Una definizione piu rigorosa:

tutto cio che riguarda il confronto tra la percezione di un segnale audio rilevata ad un

orecchio e la percezione dello stesso segnale rilevata all’altro orecchio.


Figura 1.3: Sistema di coordinate.

θ: azimuth;

φ: elevazione;

r: raggio.

In Fig. 1.3 vengono evidenziati la direzione frontale (θ = 0o e φ = 0o) e

quella posteriore (θ = 0o e φ = 180o), il piano frontale (θ = ±90o), il piano

mediano (θ = 0o) e il piano orizzontale (φ = 0o).

1.2.2 HRIR e HRTF

La propagazione delle onde sonore attraverso l’aria e un processo lineare e

permette una semplice caratterizzazione matematica dei segnali sonori di lo-

calizzazione. Per un suono lontano che si propaga nello spazio dalla sorgente

all’orecchio di un ascoltatore, la funzione di trasferimento che descrive come

l’onda sonora incidente viene modificata dal corpo dell’ascoltatore stesso e

definita Head Related Transfer Function, HRTF. La HRTF puo essere

descritta come un filtro lineare che contiene informazioni sull’interazione

dell’onda sonora con il busto, la testa e l’orecchio esterno. Va sottolineato che

non vi e alcuna informazione sulla riverberazione e sulle caratteristiche della

stanza (spazio di propagazione dell’onda sonora), pur considerati dei segnali

di localizzazione importanti per l’esternalizzazione di un suono. I percorsi di

un onda sonora per giungere alle due orecchie sono diversi, quindi definiamo

xL(t) e xR(t) come i segnali ricevuti rispettivamente dall’orecchio sinistro e

destro ottenuti dalla seguenti convoluzioni:


xL(t) =

∫hL(τ)x(t− τ)dτ (1.1)

xR(t) =

∫hR(τ)x(t− τ)dτ (1.2)

dove hL(t) e hR(t) sono le head-related impulse response (HRIR)

per l’orecchio sinistro e destro e x(t) e il segnale sonoro emesso dalla sorgen-

te. Possiamo cosı esprimere la seguente relazione: HRTF come trasformata

di Fourier di HRIR. La HRTF (HRIR) racchiude informazioni riguardanti i

fenomeni di riflessione e diffrazione e a seconda della posizione della sorgen-

te2, l’onda sonora percorrera dei percorsi diversi che si ripercuoteranno sulla

HRTF, provocando un determinato andamento spettrale. Queste modifiche

dipendono direttamente dalla struttura anatomica dell’individuo, quindi pos-

siamo affermare che, dall’unicita di ogni essere umano, deriva l’unicita della

HRTF, personale per ogni ascoltatore. Un’operazione naturale e quella di

provare a scomporre la HRTF nelle componenti che corrispondono ai contri-

buti delle diverse parti del corpo. Presentiamo in Fig. 1.4 la semplificazione

di un modello presentato da Algazi et al [4].

Figura 1.4: Scomposizione strutturale della HRTF.

Il diagramma a blocchi suggerisce una cascata degli effetti prodotti sepa-

ratamente dal busto, dalla testa ed infine dell’orecchio esterno. Possiamo cosı

ricavare i contributi delle singole componenti sottraendo opportunamente gli

effetti delle varie parti del corpo.

2Un certo azimuth, una certa elevazione e una certa distanza rispetto al sistema di

riferimento interaurale-polare centrato nella testa dell’ascoltatore.


Figura 1.5: Rappresentazione schematica di ITD e ILD (o IID).

1.2.3 La Testa

La testa puo essere considerata acusticamente come un corpo rigido a cui

sono ancorate circa alla stessa altezza e ai due lati opposti le orecchie. Si

interpone da ostacolo alla libera propagazione del suono e gia dal 1907 il fisico

Lord Rayleigh ne studio gli effetti cercando una correlazione con la posizione

della sorgente sonora. La sua teoria, chiamata duplex theory, individuo due

parametri interaurali fondamentali per una prima localizzazione del suono:

interaural time difference (ITD)

interaural level difference (ILD)

La ITD e la differenza temporale tra gli istanti di arrivo di un suono alle due

orecchie, e si misura in secondi. La ILD (in db), chiamata anche Interaural

Intensity Difference (IID), e il rapporto tra le ampiezze dei suoni in ingresso

ai due canali uditivi, dovuto all’effetto schermante ad opera della testa.

Consideriamo un modello sferico per la testa, e una sorgente sufficiente-

mente lontana. Questo scenario ci permette di considerare le onde sonore che

giungono sul piano orizzontale, come planari (v. Fig. 1.6); la distanza extra

∆x necessaria ad un raggio sonoro per raggiungere l’orecchio piu distante e

descritta dalla formula di Woodworth:


Figura 1.6: Stima della ITD nel caso di testa sferica e sorgente sonora lontana

(onda planare).

ITD =a(sin θ + θ)

c(1.3)

dove c e la velocita del suono3, a e il raggio della testa sferica e θ rappresenta

l’angolo di azimut4 . Quindi la ITD e uguale a zero quando la sorgente e

direttamente di fronte all’ascoltatore (θ = 0) ed ha il suo massimo (a/c(π/2+

1)) quando la sorgente proviene dalla parte di una delle due orecchie (θ =

π/2). Una ITD del valore di 0.6 ms si ottiene per un valore realistico del

raggio della testa, a=8.5 cm.

Queste approssimazioni rendono la ITD indipendente dalla frequenza,

mentre la ILD e fortemente dipendente dalla frequenza. Per frequenze bas-

se, che hanno lunghezza d’onda maggiore del diametro della testa, la ITD

costituisce una differenza di fase tra i segnali che giungono alle due orecchie;

questa e un’informazione molto potente per discriminare la direzione del suo-

no nel piano orizzontale. La ILD non aggiunge alcuna caratteristica in quanto

3La velocita del suono varia a seconda del mezzo e delle sue proprieta. Nell’aria a

temperatura ambiente, essa e approssimativamente 343 m/s.4Definisce la direzione dell’onda sonora nel piano orizzontale, in relazione al punto

d’ossevazione coincidente con il naso.


le basse frequenze attraversano la testa senza essere particolarmente attenua-

te sul lato opposto alla sorgente. Per alte frequenze, che hanno lunghezza

d’onda minore del diametro della testa (dell’ordine dei centimetri), essendo

l’orecchio umano sensibile soltanto alla fase e non alla differenza temporale

assoluta, la ITD risulta poco utile5. L’effetto di schermatura effettuato dalla

testa sulle onde sonore rende la ILD molto importante per le alte frequenze.

Consideriamo la seguente situazione semplificata: testa sferica con raggio a

e una sorgente sonora puntiforme posizionata alla distanza r > a dal centro

della sfera. La diffrazione di un onda acustica sulle sfera vista da un punto

scelto della sfera stessa e espressa dalla funzione di trasferimento:

Hsfera(ρ, θinc, µ) = −ρµe−iµρ

∞∑m=0

(2m+ 1)Pm(cos θinc)hm(µρ)

h′m(µ)(1.4)

dove µ e la frequenza normalizzata e ρ e la distanza normalizzata:

µ = f2πa

c(1.5)

ρ =r

a(1.6)

e Pm e hm sono rispettivamente il polinomio di Legendre di ordine m e la

funzione di Hankel sferica e θinc e l’angolo di incidenza6.

La Fig. 1.7 viene riportata la relativa risposta in ampiezza in funzione

della frequenza normalizzata per vari angoli di incidenza; possiamo esprimere

le seguenti considerazioni:

a basse frequenze la funzione di trasferimento non e dipendente dalla

direzione del suono e l’ampiezza e in prossimita dell’unita per ogni

angolo di incidenza;

5Un segnale ad alta frequenza provoca una ITD maggiore di un periodo del segnale. Se

si considerano gli onset del suono si puo migliorare l’individuazione temporale.6L’angolo tra il raggio che congiunge il centro della sfera con la sorgente puntiforme e il

raggio che congiunge il centro della sfera con il punto di misurazione. Considerando la figu-

ra 1.6, le funzioni di trasferimento in corrispondenza delle orecchie avranno rispettivamente

θ(r)inc = π/2− θ e θ(l)inc = π/2 + θ


Figura 1.7: Risposta in frequenza di una sfera rigida ideale.

ad alte frequenze (µ > 1) si fa sentire la dipendenza direzionale; nel

caso di incidenza normale (θ = 0 ) si osserva un guadagno di 6 dB, in

quanto questa situazione e paragonabile ad un’onda piana che incide

perpendicolarmente con una superficie piana rigida.

all’aumentare dell’angolo di incidenza il guadagno diminuisce; spostan-

do la sorgente dalla parte controlaterale7 della testa si ha attenuazione

della pressione sonora, con la presenza di varie oscillazioni dovute alla

propagazione delle onde in diverse direzioni attorno la sfera;

la risposta minima non corrisponde ad una sorgente situata esattamente

dalla parte opposta (θinc = π); in un particolare punto si nota l’effetto

bright spot dove tutte le onde, che si propagano attorno alla sfera,

arrivano in fase.

7Si dice che la sorgente e posizionata nella parte ipsilaterale della sfera (della testa) se

l’onda sonora incide normalmente la superficie in un punto appartenente all’emisfero che

ha il polo nel punto di osservazione; controlaterale se l’onda incide la sfera nell’emisfero

opposto.


Figura 1.8: Anatomia della pinna.

Questo modello sferico e l’unico trattabile analiticamente, ma e anche

il meno preciso. Per questo motivo vi sono molti lavori in letteratura che

cercano di migliorare l’accuratezza del modello, considerando, per esempio,

considerando la testa come un ellissoide[30] e usando tecniche di ray-tracing

per studiarne il comportamento.

1.2.4 L’Orecchio

L’orecchio esterno e formato dalla pinna e dal canale uditivo che si congiunge

al timpano. Ad una prima approssimazione il canale uditivo si comporta co-

me un risonatore monodimensionale, mentre gli effetti dovuti alla pinna sono

piu articolati e saranno oggetto di studio in questa tesi8. Le cavita risonanti

della pinna amplificano alcune frequenze mentre la sua geometria genera at-

tenuazioni e cancellazioni di altre frequenze. La risposta in frequenza della

pinna e altamente dipendente dalla posizione della sorgente sonora relativa

all’orecchio.

Shaw[23] espresse l’orecchio esterno come un sistema acustico leggermente

smorzato e descrisse i modi normali di risonanza.

8Da ora in poi con il temine orecchio esterno ci riferiremo principalmente alla pinna

con le sue caratteristiche antropometriche, trascurando il contributo del canale uditivo.


Modo 1 presente indipendentemente dalla direzione e appare mediamente

attorno alla frequenza di 4.2 kHz;

Modo 2 appare alla frequenza 7.1 kHz, dove si osserva una separazione tra

due zone di pressione vicino alla cymba(v. Fig 1.8); questo modo e

meglio stimolato per un angolo di elevazione di 68 nel piano mediano.

Modo 3 appare alla frequenza 9.6 kHz e ha due transizioni tra zone di pres-

sione. Una zona positiva attorno alla fossa che passa ad una negativa

all’interno e attorno alla cymba o alla crus helias ; piu giu si nota un’al-

tra zona di pressione positiva nella concha. Questo modo e meglio

stimolato per un angolo di elevazione di 73 nel piano mediano.

Modo 4 ha 3 transizioni di pressione e appare alla frequenza 12.1 kHz. Vi

e una zona negativa attorno alla fossa che diventa positiva muovendosi

verso la cymba. Questo modo e meglio stimolato per un angolo di

elevazione di -60 nel piano mediano.

Modo 5 e 6 appaiono rispettivamente alle frequenze 14.4 kHz e 16.7 kHz.

Sono simili per il patter delle transizioni di pressione e sono meglio

stimolate per un angolo di elevazione medio di 7 nel piano mediano.

Il primo modo e quello dominante e dipende principalmente dalla pro-

fondita della conca. I modi 2 e 3 sono risonanze trasversali e si possono

descrivere come una coppia di risonanze verticali che viene eccitata ad alti

angoli di elevazione. In maniera analoga i modi 4, 5 e 6 si possono descrivere

come risonanze orizzontali perche meglio eccitate ad elevazioni prossime allo

zero. L’onda sonora per arrivare al timpano percorre piu di un cammino,

uno e quello diretto e gli altri, piu lunghi, dovuti alla riflessione della pinna.

A basse frequenze la lunghezza d’onda e piu grande della dimensione della

pinna quindi il segnale diretto e quello riflesso arrivano in fase, mentre ad alte

frequenze il segnale riflesso e sfasato rispetto a quello diretto e genera un’in-

terferenza distruttiva. Nello spettro si formano cosı dei notch in relazione

alla posizione della sorgente e all’antropometria dell’orecchio. La pinna viene


Figura 1.9: I sei modi normali di risonanza della conca umana.

utilizzata anche per discriminare la front-back confusion, in quanto maschera

e attenua i suoni che provengono da dietro l’ascoltatore in quanto sporgente

dalla testa. Un concetto fondamentale che verra ampiamente trattato nel

corso della tesi riguarda l’individualita della forma dell’orecchio che genera

caratteristiche del tutto personali nella HRTF. La trattazione approfondita

sulla pinna continua nella descrizione dei modelli strutturali e nei seguenti

capitoli di analisi.

1.2.5 Busto e Spalle

Gli ultimi elementi che, assieme alla testa e all’orecchio esterno, contribui-

scono al filtraggio del segnale sonoro che arriva al timpano sono il busto e

le spalle. Provocano riflessioni che vanno a sommarsi con il suono diretto e

schermano i raggi sonori che provengono da dietro l’ascoltatore.

Tutti i lavori storici a cui si fa riferimento nel seguito sono citati in [31],

che ne fornisce una rassegna esaustiva. Il ruolo di disturbo a basse frequenze

dell’onda incidente e stato riconosciuto gia dal 1944 ad opera di Hanson e


Figura 1.10: Cono d’ombra del busto.

successivamente meglio studiato da Kuhn e Guernsey nel 1983. Tuttavia,

l’effetto del busto e relativamente debole e l’importanza percettiva alle basse

frequenze non e del tutto compresa. Per esempio Theile e Spikofski (1982)

hanno concluso dai loro esperimenti che il busto non fornisce un segnale signi-

ficativo per la discriminazione front/back, informazione che potrebbe essere

contenuta solo nella parte alta dello spettro. Mentre Asano et al. (1990) ri-

conoscono anche nella parte bassa dello spettro il contributo importante del

busto per la tale discriminazione. Gli effetti del busto nella localizzazione

verticale nel piano mediano sono stati studiati sistematicamente da Gardner

(1973). Mise alla prova il senso di localizzazione di alcuni soggetti; ne osservo

la diminuzione nel caso vengano rimosse alla sorgente sonora le alte frequen-

ze e pote scoprire come alcuni individui fossero in grado di localizzare suoni

provenienti da altoparlanti posti nel piano mediano frontale anche se privi

di contenuto spettrale sopra i 4 kHz. Confrontando i risultati ottenuti per

elevazioni tra i 18 e i -18 concluse che l’orecchio esterno non ha influenza

sotto i 3.5 kHz e che il busto introduce dei segnale spaziali per frequenze tra

0.7 e 3.5 kHz.

La geometria del busto e delle spalle e particolarmente complicata e le


approssimazioni piu semplici ed efficaci sono quelle di considerare un busto

ellisoidale[31] o sferico[32][33]. L’approssimazione della testa per mezzo di

una sfera o di un ellisoide aggiunta alle approssimazioni per il busto danno

origine ai modelli snowman.

In generale si ha che il tempo di ritardo tra il raggio sonoro diretto e quello

riflesso dal busto e maggiore quando la sorgente e sopra alla testa. Man mano

che si riduce l’elevazione il ritardo temporale si riduce fino a diventare zero

se il raggio dalla sorgente all’orecchio e tangente al busto. L’insieme dei

raggi tangenti al busto formano quello che viene chiamato cono d’ombra del

busto (v. Fig. 1.9). Una sorgente di suono fuori dal cono d’ombra induce

una riflessione ad opera del busto; la riflessione scompare quando la sorgente

entra all’interno del cono d’ombra, dove viene schermata dal busto stesso. Il

ritardo tra il raggio diretto e quello riflesso dal busto non varia di molto se la

sorgente si muove sul piano orizzontale lungo una circonferenza , soprattutto

se il raggio della testa e molto piu piccolo rispetto al raggio della circonferenza

che sta percorrendo la sorgente.

Nel dominio della frequenza le riflessioni del busto si manifestano come

un comb filter che introduce notch periodici nello spettro. Le frequenze dei

notch sono strettamente legate al tempo di ritardo ed hanno un pattern

caratteristico al variare dell’elevazione della sorgente sonora. Vi sono anche

effetti di diffrazione e diffusione sonora che producono forti attenuazioni per

le alte frequenze9.

Tuttavia gli effetti acustici del busto e delle spalle non sono cosı forti

come quelli causati dall’orecchio esterno, anche se la loro importanza e indi-

scutibile per le basse frequenze, dove il segnale sonoro ha la maggior parte

della sua energia e dove la risposta in frequenza dell’orecchio esterno e piat-

ta. Possiamo concludere dicendo che gli effetti spettrali di busto e orecchio

esterno si completano a vicenda.

9Per lunghezze d’onda minori o uguali alla taglia del busto.

Capitolo 2

Modellazione HRIR e HRTF

Si possono generare forti effetti di localizzazione spaziale di una sorgente so-

nora mediante la convoluzione di un segnale monoaurale con le HRIR delle

due orecchie e il risultato verra opportunamente compensato da delle cuf-

fie e o degli altoparlanti cross-talk-canceled. Per accrescerli ulteriormente e

sensato utilizzare delle HRIR spazio-varianti in relazione ai movimenti della

testa. Purtroppo questo approccio, pur essendo concettualmente semplice

e computazionalmente troppo dispendioso. Un sistema tipico implementa

un movimento della testa come un’interpolazione real-time tra un insieme

piuttosto grande di coefficienti FIR ricavati da un processo di identificazione

delle HRIR misurate, indicizzato per azimuth ed elevazione.

Per produrre effetti convincenti (soprattutto d’elevazione), le HRIR devo-

no essere misurate separatamente per ogni ascoltatore; questa e una soluzione

molto dispendiosa in termini di tempo e risorse e per semplificare i sistemi

audio sono state proposte numerose ricerche che mirano a sostituire le HRIR

misurate con modelli computazionali. Gli effetti di azimuth possono venir

prodotti semplicemente introducendo le appropiate ITD e ILD. Gli effet-

ti di elevazione possono venir prodotti introducendo opportuni notch della

HRTF, a seconda dele differenze tra persone che rendono il meccanismo di

localizzazione verticale difficile da simulare. Inoltre questi sistemi soffrono di

altri problemi quali la percezione del suono troppo vicino e una confusione

23

24 CAPITOLO 2. MODELLAZIONE HRIR E HRTF

front/back.

In questo capitolo tratteremo le principali tecniche, proposte in lettera-

tura, per ottenere delle HRIR/HRTF personali. Ci soffermeremo sul metodo

dei modelli strutturali che utilizziamo in questa tesi.

2.1 HRIR e HRTF personali

La HRIR e la HRTF sono funzioni in quattro variabili: le tre coordinate

spaziali e il tempo o la frequenza. Entrambe queste funzioni sono piuttosto

complicate e variano significativamente da persona a persona. Da molto

tempo si cerca di sostituire le HRIR/HRTF con funzioni approssimate, ma

la progettazione di tali sistemi di identificazione viene complicata da alcuni

problemi importanti:

1. si incontrano difficolta nell’approssimazione di effetti di propagazio-

ne e diffrazione dell’onda sonora mediante semplici modelli di ordine

contenuto che riescano ad essere accurati;

2. non si riesce a calcolare la HRIR/HRTF nelle quattro variabili perche

queste interagiscono strettamente tra loro. Se i parametri del modello

variassero in maniera complessa al variare della posizione spaziale si

rischierebbe di memorizzare questi cambiamenti in una grande mole

di dati, vanificando lo sforzo fatto per eliminare la tabelle di HRIR

misurate.

3. non esiste un criterio quantitativo per misurare la capacita di un’ap-

prossimazione nel catturare le informazioni direzionali che sono percet-

tivamente rilevanti. Esistono solo test psico-acustici per valutare un

modello, e molto spesso non sono provvisti di uno standard.

4. Un’approssimazione che funziona bene per un individuo potrebbe non

funzionare per un’altro. Questo genera la necessita di creare in maniera

diversificata l’identificazione per ogni ascoltatore.

2.1. HRIR E HRTF PERSONALI 25

Vengono ora trattate le seguenti tecniche per produrre HRIR/HRTF1

personali: 1) misurazione acustica; 2) interpolazione su database; 3) modelli

poli-zeri; 4) espansione in serie; 5) metodi numerici; 6) modelli struttura-

li. Per ogni tecnica verranno elencate piu possibilita e ci soffermeremo su

alcune varianti in modo da orientarne la comprensione globale. Ai modelli

strutturali verra dedicata la seconda parte del capitolo.

2.1.1 Misurazione acustica

Una mole considerevole di dati riguardanti la HRTF e stata ottenuta da mi-

surazioni acustiche in ambienti anecoici o quasi-anecoici e parecchi database

sono stati resi pubblici[3] per uso didattico e ricerca. Esistono anche alcu-

ni prodotti commerciali per misurare le HRTF personali, tuttavia restano

sistemi troppo costosi e sofisticati per un uso domestico. Per avere una mi-

surazione accurata e ad alta risoluzione il processo di acquisizione dei dati

necessita di apparecchiature speciali e diventa molto dispendioso in termini

di tempo.

2.1.2 Interpolazione su database

Consideriamo in questo paragrafo alcune metodologie per la riduzione delle

misurazioni acustiche. Utilizziamo l’interpolazione per ricondurre un ascol-

tatore ad un soggetto del database (evitando nuove acquisizioni di dati) e

per limitare le misurazioni acustiche per ogni soggetto.

Interpolazione di individuo Avendo a disposizione un database di HRTF,

risulta un alternativa attraente cercare dei modi non acustici per misu-

rare la similarita di un ascoltatore con i soggetti gia presenti, ricercando

la HRTF che meglio gli assomiglia. Una possibilita consiste nel ridurre il

numero di prototipi possibili tra i quali l’ascoltatore puo iterativamente

1Da questo momento potremmo utilizzare i due termini interscambiandoli, tenendo

presente che la HRTF e la trasformata di Fourier della HRIR e viceversa.


scegliere quello con il quale sente meglio. Un’altro approccio conside-

ra le correlazioni tra le caratteristiche della HRTF e l’antropometria.

Se le misurazioni antropometriche possono venir fatte in tempi brevi,

usando tecniche di computer vision, il miglior matching nel database

corrispondera al soggetto che piu assomiglia all’ascoltatore. Purtrop-

po le caratteristiche della HRTF sono molto sensibili a cambiamenti

relativamente piccoli, soprattutto per la forma dell’orecchio esterno, e

trovare il matching migliore nel database rimane ancora un problema

aperto.

Interpolazione spaziale Le misurazioni acustiche per la creazione di un

database di HRTF vengono fatte considerando un insieme finito di po-

sizioni e quando vi e la necessita di effettuare il rendering per posizioni

intermedie si ricorre all’interpolazione tra i dati presenti. Se la HRTF

intermedia venisse ricondotta alla HRTF piu vicina verrebbero genera-

ti degli artefatti acustici. Per effettuare l’interpolazione direttamente

sui campioni delle HRIR si puo utilizzare il metodo bilineare[39]: vie-

ne calcolata la risposta di un punto, individuato dalla coppia di valori

(θ, φ), come media pesata delle quattro risposte piu vicine, presenti nel

database. Se l’insieme di HRIR e misurato su una griglia sferica con

passo di campionamento spaziale θgrid e φgrid, la risposta all’impulso

stimata, h, nel punto (θ, φ), puo essere calcolata come:

h[n] = (1− cθ)(1− cφ)ha[n] + cθ(1− cφ)hb[n] + cθcφhc[n] + (1− cθ)cφhd[n]

(2.1)

dove hα[n](α = a, b, c, d) sono le HRIR associate ai quattro punti vicini

al punto desiderato. I parametri cθ e cφ sono calcolati come:

cθ =θ mod θgrid

θgrid, cφ =

φ mod φgridφgrid

(2.2)


Figura 2.1: a) interpolazione bilineare delle HRIR; b) interpolazione di zeri

per modelli poli-zeri sintetici di HRTF

Si possono applicare numerosi miglioramenti a questa tecnica[38], pre-

stando attenzione all’introduzione da parte dell’interpolazione di arte-

fatti udibili dovuti ai cambiamenti di fase, al sampling rate, etc.

L’interpolazione puo essere adoperata anche con HRTF sintetiche mo-

dellate tramite filtri poli-zeri(v. paragrafo successivo). Nel semplice

caso di due filtri FIR Hα(z)(α = a, b) con q zeri della forma:

Hα(z) = 1 +

q∑k=1

bα,kz−k =

q∏k=0

(1− cα,kz−1), bα,0 = 1 (2.3)

gli zeri di entrambi i filtri sono ordinati secondo la fase. Il filtro risul-

tante dall’interpolazione, H(z) =∏q

k=0(1− ckz−1), puo essere ottenuto

associando gli zeri per prossimita angolare e calcolando lo zero come

ck = (1 − ρ)ca,k + ρcb,k con k = 1,...,q. Essendo Hα a fase minima,

anche il filtro interpolato sara a fase minima. Nel caso di filtri IIR,

l’interpolazione dei poli e degli zeri diventa piu complicata e necessita

di algoritmi specifici per garantire la stabilita del risultato.


2.1.3 Metodi Numerici

In teoria le HRTF potrebbero venir calcolate risolvendo l’equazione delle

onde, soggette alle condizioni al contorno imposte dalla superficie del corpo

su cui si infrangono. In pratica una soluzione analitica risulta impossibile e

si utilizzano metodi numerici come:

teoria della diffrazione di Kirchhoff nel dominio della frequenza;

boundary-element methods nel dominio della frequenza;

boundary-element methods incrementato dall’espansione multipolo;

metodo alle differenze finite nel dominio del tempo;

sintesi differenziale della pressione;

metodi ispirati alla computer graphics basati su ray tracing.

Le due principali problematiche che limitano l’utilizzo di queste tecni-

che per ottenere HRTF personali sono: la necessita di avere delle superfi-

ci mesh molto accurate soprattutto per modellare la pinna, e l’alto costo

computazionale per simulazioni su banda larga.

Ci soffermiamo sul metodo alle differenze finite nel dominio di tempo

(FDTD) adoperato in [22]. Yee propose, nel 1966, questo algoritmo nu-

merico per risolvere efficientemente il problema al contorno che coinvolge le

equazioni di Maxwell (riportate di seguito).

∇× E =−∂H∂t

1

µ(2.4)

∇ · E = ρ (2.5)

∇×H =−∂D∂t

+ J (2.6)


∇ ·B = 0 (2.7)

Dalla relazione temporale e spaziale tra il campo magnetico e il campo elet-

trico, Yee penso di esprimere il campo attuale utilizzando gli stati passati.

Si puo porre quest’idea in analogia con una macchina a stati finiti in cui

lo stato attuale dipende solamente da quello precedente ed influenzera so-

lamente lo stato successivo. Nel il periodo storico in cui e stato proposto

il metodo, il costo computazionale era insostenibile, ma negli ultimi dieci

anni e stato riscoperto grazie alla sua semplicita e soprattutto grazie all’au-

mento della potenza di calcolo della macchine odierne. L’idea alla base del

metodo FDTD puo essere facilmente applicata alla propagazione delle onde

sonore nello spazio. Per esempio, nella simulazione di campi aperti, dove si

vogliono eliminare le riflessioni, e necessario definire la spazio di calcolo per

avere delle griglie di calcolo limitate. Vengono introdotte delle condizioni di

assorbimento al contorno (Absorbing Boundary Condition, ABC ) all’inter-

no dello spazio, che simulano un materiale che genera riflessioni ridotte. I

parametri che generalmente si possono decidere e controllare sono: la taglia

della griglia spaziale, il passo di campionamento spaziale della griglia, il pas-

so di campionamento temporale (la lunghezza degli incrementi temporali tra

uno stato e l’altro), il punto di osservazione, il numero di sorgenti sonore, la

banda spettrale da analizzare e il tipo di ABC. La scelta di questi parametri

influenza notevolmente la performance e l’accuratezza della simulazione.

2.1.4 Modelli Poli-Zeri

La HRTF per una certa direzione (θ, φ) puo essere approssimata con un

modello poli-zeri o modello ARMA (AutoregRessive Moving Average) con

funzione di trasferimento nella seguente forma:

H(z) =b0 +

∑qk=1 bkz

−k

1−∑p

k=1 akz−k =

B(z)

A(z)(2.8)

dove i coefficienti bk e ak dipendono da θ, φ e possono venire individuati da

numerosi sistemi di identificazione fino ad oggi utilizzati.


Nel caso di p = 0, la H(z) e un modello FIR tutti zeri. Una tecnica

immediata per ottenerlo consiste nella finestratura della risposta all’impulso

h(n) corrispondente a H(z).

h(n) =1

N

N−1∑k=0

H(ejwk)ejwkn (2.9)

il segnale nel tempo sara h(n) = h(n)w(n). La lunghezza di w(n) e inferiore

inferiore ad h(n) a seconda dell’approssimazione desiderata e, in relazione al

tipo di finestra utilizzata, si avranno diversi effetti sulla HRTF. Questo meto-

do puo essere modificato introducendo dei pesi nel dominio della frequenza,

in modo da modellare la risoluzione spettrale non uniforme dell’orecchio;

in modo analogo le frequenze si possono distorcere con metodi di warping,

come bilinear warping, che sostituisce un’unita di ritardo con una sezione

passa-tutto del primo ordine:

z−1 ← D1(z) =z−1 − λ1− λz−1

(2.10)

con λ denominato coefficiente di warping. Il procedimento e reversibile, viene

chiamato unwarping, e si realizza sostituendo −λ a λ. In Fig. 2.2 si puo

osservare l’andamento della curva di warping al variare di λ.

Nel caso di q = 0, la H(z) e un modello IIR tutti poli. Un metodo

per calcolare i coefficienti ak che meglio approssimano H(z) e la Linear Pre-

diction. In generale, il design di un filtro approssimatore si riconduce alla

minimizzazione di una certa funzione d’errore nella forma:

E(ejw) = H(ejw)− H(ejw) (2.11)

da cui possiamo citare alcuni criteri d’errore comunemente utilizzati:

Least-Squares Error: EH − H = 1L

∑L−1k=0 (H(ωk) − H(ωk))

2 con ωk =

2kπ/LFs;


Figura 2.2: Warping delle frequenze per una sezione passa-tutto D1(z) del

primo ordine al variare di λ. Le frequenze sono normalizzate rispetto alla

Nyquist rate.

Log-Magnitude responce Error: ElogH − H = 1L

∑L−1k=0 (ln|H(ωk)| −

ln|H(ωk)|)2 criterio di errore guidato da una motivazione percettiva;

Weighted Least Square: EH− H = 1L

∑L−1k=0 (H(ωk)−W (ejw)H(ωk)|)2

con W (ejw) funzione peso per modellare la risoluzione spettrale dell’a-

scoltatore.

In Fig. 2.3 e riportata la tabella presentata da Houpaniemi in [37] del-

la tassonomia negli anni ’90 per design di filtri che approssimano la HRTF.

Come si puo osservare i risultati sono molto variegati. Questo e dovuto al-

l’utilizzo, per alcuni di questi studi, di errori spettrali definiti teoricamente

o per ispezione visiva, per altri di test d’ascolto.2 Un’altra motivazione che

rende cosı diversi questi progetti risiede nella varieta di tipologie di dati usati

2Purtroppo non esistono test d’ascolto definiti formalmente e che danno un risultato

statistico probabilmente uniforme.


Figura 2.3: Tassomania anni ’90 dei design di filtri per HRTF binuarali. I

filtri con [Study = Empirical? ] esibiscono risultati dedotti da test informali

nei singoli studi. Le HRTF possono venir equalizzate rispetto alle condizioni

di campo vicino oppure rispetto ad un particolare tipo di cuffie e i dati posso-

no riferirsi a un manichino (KEMAR) oppure essere personali o impersonali.

Infine puo essere stata applicata o meno una ricostruzione a fase-minima.

I ricercatori hanno utilizzato varie tecniche di filter design, sistemi di iden-

tificazione e reti neurali per cercare una corrispondenza tra i modelli multi-

parametro e dati sperimentali. Sfortunatamente, molti dei filtri progettati

hanno coefficienti che si legano con funzioni molto complicate dell’azimu-

th e dell’elevazione: modelli che hanno un numero sufficiente di parametri

per catturare i segnali di localizzazione personali non forniscono vantaggi

computazionali significativi.

2.1.5 Espansioni in serie

Sebbene la HRTF possa parere complicata, si puo ipotizzare, su basi fisiche, la

sua completa determinazione tramite un insieme relativamente piccolo di pa-


Figura 2.4: Un insieme di dati bidimensionale e due vettori di base (o assi

principale) PC1 e PC2 estratti utilizzando la PCA.

rametri fisici: il raggio medio della testa, l’eccentricita della testa, il diametro

massimo della pinna, la profondita della conca, ecc. La dimensionalita della

HRTF risulterebbe piccola, nonostante la sua complessita potrebbe portare

ad una visione errata dei parametri. Uno dei metodi che e stato utilizzato

per fare chiarezza e la principal component analysis, PCA (anche conosciuta

come trasformazione di Karhunue-Loeve) applicata alla HRTF complessa o

al logaritmo dell’ampiezza della HRTF. Si produce un insieme direzionalmen-

te indipendente di funzioni base e un insieme direzionalmente dipendente di

pesi per combinare le funzioni di base che rappresentano lo spazio dei pa-

rametri. Sfortunatamente, nella sintesi real-time, viene richiesto uno sforzo

computazionale considerevole; quando vi e un movimento della testa dell’a-

scoltatore o uno spostamento della sorgente sonora, i pesi si relazionano con

funzioni complicate a seconda dell’azimuth e dell’elevazione, che devono venir

catalogate in tempi brevi.

Per dare un’idea piu chiara di come operano questi metodi, spieghere-

mo brevemente la PCA, procedura statistica molto popolare che permette

di ridurre la dimensionalita di una grande mole di dati e mantenere il piu

possibile coerenti le variazioni tra gli stessi.


Si vuole rappresentare un insieme di M vettori a valori reali x1,...,xM ,

di dimensione N, attraverso la proiezione dei dati su una linea che passi

attraverso la media, di equazione

xk = m + ake, k = 1, ...,M (2.12)

dove e e un vettore unitario nella direzione della linea, e gli scalari ak corri-

spondono alla distanza di ogni punto dalla media m = 1M

∑Mk=1 xk. Ottenia-

mo l’insieme ottimo di coefficienti ak minimizzando l’errore quadratico nella

funzione

E(a1, ..., aM , e) =M∑k=1

‖(m + ake)− xk‖2 =M∑k=1

‖ake− (xk −m)‖2

=M∑k=1

a2k‖e‖2 − 2

M∑k=1

akeT (xk −m) +

M∑k=1

‖xk −m‖2.

(2.13)

Sapendo che ‖e‖ = 1, si pone la derivata parziale rispetto ad ak pari a zero:

∂E(a1, ..., aM , e)

∂ak= 2ak − 2eT (xk −m), (2.14)

i coefficienti ottimi si trovano proiettiamo il vettore xk sulla linea e che passa

attraverso la media, ottenendo:

ak = eT (xk −m). (2.15)

Per trovare la miglior direzione per e definiamo, come prima cosa, la matrice

di covarianza dell’insieme di vettori,

S =1

M

M∑k=1

(xk −m)(xk −m)T . (2.16)


e cerchiamo di minimizzare l’errore, trovando il massimo della funzione f(e) =

eTSe, con il vincolo ‖e‖ = 1. Illustriamo brevemente i passaggi algebrici che

ci portano a questo risultato:

E(a1, ..., ak, e) =M∑k=1

a2k − 2

M∑k=1

a2k +

M∑k=1

‖xk −m‖2

= −M∑k=1

[eT (xk −m)]2 +M∑k=1

‖xk −m‖2

= −M∑k=1

eT (xk −m)(xk −m)Te +M∑k=1

‖xk −m‖2

= −MetSe +M∑k=1

‖xk −m‖2;

(2.17)

Ora si puo utilizzare il metodo dei moltiplicatori di Lagrange3

ed ottenere la seguente funzione di Lagrange e il suo gradiente

L(e, λ) = eTSe− λ(ete− 1) (2.19)

∇eL(e, λ) =∂L(e, λ)

∂e= 2Se− 2λe (2.20)

Concludendo, i punti di e che massimizzano f(e) sono quelli che soddisfano:

Se = λe, (2.21)

3Procedimento matematico per trovare gli estremi di una funzione f(x) soggetta ad

uno o piu vincoli nella forma g(x) = 0. Si costruisce la funzione Lagrangiana

L(x, λ) = f(x) + λg(x), (2.18)

dove λ e uno scalare chiamato moltiplicatore incognito di Lagrange e si cercano gli zeri del

gradiente ∇xL(x, λ)


e deve essere un autovettore della matrice di covarianza S, con rispettivo au-

tovalore λ. Il miglior autovettore possibile secondo il criterio della differenza

quadratica minima e quello a cui corrisponde l’autovalore della matrice di

covarianza maggiore che garantisce la massimizzazione di eTSe = λeTe = λ.

Il risultato puo essere esteso a proiezioni ad uno scenario multidimensio-

nale. Per una proiezione q-dimensionale attraverso la media dei campioni si

otterra:

x = m +

q∑i=1

ak,iei; (2.22)

dove i coefficienti ak,i si dicono componenti principali. Si procede quindi in

maniera analoga al caso unidimensionale proiettando i dati su i q autovettori

di S (assi principali) corrispondenti ai q autovalori piu grandi. Se si usano

tutti gli autovettori (q = M) si otterranno i dati originali; il numero di

assi principali q, necessari ad una adeguata rappresentazione dei dati, e una

funzione dell’ammontare di ridondanza4 o correlazione presente nel dataset.

Il primo ricercatore che adopero la PCA sui dati di HRTF fu Martens

(1987); descriviamo brevemente un possibile utilizzo della PCA in questo

campo. Siano H(θk, φk, ωj) le HRTF su M direzioni θk, φk(k = 1...M) e

su N frequenze ωj(j = 1...N). Costruiamo M vettori N-dimensionali xk =

log|H(θk, φk, ωj)|, in modo da sfruttare percettivamente l’approssimazione

logaritmica dell’ampiezza. Il risultato e composto da un insieme di q basi ei,

con ei,j = ei(ωj), in modo da approssimare la direzione (θk, φk) nella seguente

maniera:

log|H(θk, φk, ωj)| ∼q∑i=1

ai(θk, φk)ei(ωj) (2.23)

Gli studi effettuati su questo metodo hanno individuato le prime cinque

funzioni base (q=5) come rappresentative dell’insieme delle ampiezze delle

4Piu ridondanza e presente, piu si potra diminuire q.

2.2. MODELLI STRUTTURALI 37

HRTF; il livello di correlazione tra le HRTF sintetizzate e quelle misurate e

buono. Tuttavia evidenziamo che le dipendenze spaziali e di frequenza sono

state disassociate semplificando molto la computazione e che le grandezze

antropometriche collegate alle funzioni base sono di difficile misurazione.

2.2 Modelli Strutturali

L’ultimo approccio che presentiamo e quello che ha guidato le ricerche di

questo lavoro. La HRTF puo venir modellata sulla base di un’analisi sem-

plificata della fisica della propagazione e della diffrazione delle onde sonore.

Possiamo vedere il modello sferico di Lord Rayleigh (v. par. 1.2.3) come un

primo passo in questa direzione.

Uno dei lavori piu ambiziosi e stato fatto nel 1983 da Genuit, dove rilevo

34 grandezze antropometriche che caratterizzano la forma delle spalle, della

testa e della pinna e le uso per approssimare la HRTF come composizione

delle corrispondenti funzioni di trasferimento. Da questi studi prese forma

il concetto di modello strutturale. Purtroppo la documentazione sul mo-

dello di Genuit e scarsa, non sono chiare le relazioni tra alcuni parametri

del modello in relazioni alle corrispondenti grandezze antropometriche e non

ci sono informazioni per comprendere se le 34 grandezze sono necessarie o

sufficienti per la determinazione della HRTF. Ma non vogliamo soffermarci

sui lati in ombra di questo primo modello, bensı vogliamo sottolineare le

caratteristiche importanti che aprono gli orizzonti dello studio dei modelli

strutturali.

ogni componente considerata nel modello e ritenuta responsabile di

fenomeni fisici significativi e identificati;

il modello e adatto ad applicazioni real-time perche computazionalmen-

te economico;


Figura 2.5: Un modello strutturale proposto da Brown e Duda nel 1998

per approssimare la HRTF. Moduli separati descrivono la schermatura della

testa, le riflessioni delle spalle e gli effetti della pinna.

si possono relazionare i parametri del modello con le grandezze cor-

poree, dando la possibilita di determinare HRTF personali ricavate

dall’acustica e dell’antropometria;

non e un modello a cascata o multi-cammino, ma un modello naturale

che si adatta al problema.

La HRTF viene modellata come una combinazione di filtri collegati in

maniera congruente al meccanismo di ascolto umano: il contributo della te-

sta, del busto, delle spalle e le riverberazioni dovute all’acustica della stanza

vengono processate indipendentemente in parallelo per poi essere incanalate

verso il filtro che racchiude le caratteristiche della pinna. Separare gli ef-

fetti delle varie componenti anatomiche in strutture indipendenti di filtri e

un’approssimazione euristica che non considera le interazioni tra le onde che

si infrangono da una all’altra, ma che porta comunque a buoni risultati. Un

esempio di modello strutturale sviluppato da Brown e Duda [24] e riportato

in Fig. 2.5.

Il modello denominato snowman model (v. par. 1.2.5) puo essere visto

come un esempio di estensione strutturale del modello sferico di Lord Raylei-


Figura 2.6: Il modello strutturale snowman, chiamato anche head-and-torso

(HAT) model.

gh [32]. Un busto sferico e posizionato sotto una testa sferica che considera

anche le riflessioni delle spalle. I parametri coinvolti sono solamente quattro

(a,h,b, e δ, v. Fig. 2.6) e forniscono una buona approssimazione del com-

portamento a basse frequenze della HRTF, in quanto non viene considerato

il contributo della pinna. I parametri a raggio della testa, b raggio del buto

e h altezza del collo sono facilmente misurabili.

Il filtro di Brown e Duda che rappresenta l’effetto schermante della testa

(da cui il nome shadowing filter) e composto da un polo e da uno zero che

approssimano la soluzione analitica (formula 1.4) riducendo la forza del bright

spot (v. par. 1.2.3). La funzione di trasferimento e data da:

H(s, θ, a) =ατs+ 1

τs+ 1(2.24)

dove il parametro α dipende da θ e il costante di tempo dipende da a, con

τ =2a

c(2.25)


Siccome H(0, θ, a) = 1, il DC gain e 0 db per qualsiasi valore di θ o a. Il

parametro α(θ) e il guadagno asintotico per le alte frequenze con formula:

α(θ) =

(1 +

αmin2

)+

(1− αmin

2

)cos

(θ

θminπ

)(2.26)

con il valore αmin = 0.1 e θmin = 5π/6rad o 150 si raggiunge una buo-

na corrispondenza con i risultati teorici e un bright spot attenuato. Una

proprieta fondamentale del shadowing filter e la dipendenza della posizio-

ne dei poli solo dal raggio della testa dell’ascoltatore. Per una testa di

grandezza media, raggio di a = 8.75 cm, si ha una frequenza per il polo

di fc = 1/2πτ = c/4πa ≈ 312Hz. La posizione dello zero varia a seconda

dell’angolo di osservazione. Ad un angolo θflat, dove α = 1, lo zero e il polo

si cancellano e la risposta in frequenza e piatta. Dall’equazione (2.26) si puo

ricavare l’angolo critico:

θflat = θmin

(1

2+

1

πsin−1 αmin

2− αmin

)(2.27)

Se αmin = 0.1 e θmin = 150, θflat ≈ 77.5. Se θ > θflat, il filtro taglia

le alte frequenze, mentre se θ < θflat le aumenta. La Fig. 2.7.b mostra la

famiglia completa delle curve della risposta in frequenza per il filtro sintetico.

Il taglio massimo alle alte frequenze e attorno ai 20db (=20log10αmin) e

avviene a θ = θmin(150). Se θ = 0 le alte frequenze sono aumentate di

6db. Quando θ = π la risposta trova qualche difficolta nell’approssimare la

soluzione analitica5.

I modi per migliorare questo semplice modello sono numerosi, ne citiamo

alcuni:

inserire un fractional delay filter FITD(θ, z) che considera l’ITD;

trattare la testa come un ellissoide;

5Per avere un fit migliore con la soluzione teorica bisognerebbe utilizzare una equazione

per α piu complicata che renderebbe vano lo sforzo di semplificazione del modello.


Figura 2.7: Risposta in frequenza di un semplice modello di testa sferica del

primo ordine al variare dell’angolo di osservazione θ. (a) soluzione teorica

esatta, (b) approssimazione con filtro sintetico

posizionare i punti d’osservazione (le orecchie) in maniera realistica,

non per forza lungo il diametro della testa.

Per quanto riguarda il busto si puo semplicemente rappresentare la rifles-

sione principale, per mezzo di un ritardo che dipende dall’azimuth e dall’e-

levazione τ (t)(θ, φ). Il busto si puo modellare utilizzando un unico filtro a

ritardo frazionario:

Htorso = g(t)Fτ (t)(θ, φ, z) (2.28)

dove g(t) e il coefficiente di riflessione. Facciamo notare come questo modello

sia valido solo per φ > 0; se la sorgente scende in elevazione ci sara un punto

in cui subentrera l’effetto schermante del busto.

Gli effetti della pinna sono molto piu complessi da modellare perche e

difficile estrarre i suoi parametri antropometrici dai dati misurati. Nel pros-


Figura 2.8: Modello strutturale per la pinna proposto da Batteau e

formalizzato da Watkins

simo paragrafo, concluderemo questo capitolo esponendo i vari modelli per

pinna studiati e approfonditi nella letteratura.

2.3 PRTF e antropometria

Per PRTF si intende Pinna Related Transfer Function, cioe la risposta

in frequenza della pinna isolata dal resto del corpo e dalla testa6.

Negli anni ’70 Batteau mostro sperimentalmente che la pinna forniva dei

segnali molto importanti per la discriminazione dell’elevazione,oltre ad es-

serlo per l’azimuth. Propose un semplice modello dove la parete della conca

e l’elice (helix o rim) agiscono da riflettori. Le onde riflesse interferiscono

con il cammino diretto dell’onda sonora producendo dei profondi notch nello

spettro. Ad ogni prima riflessione e associato un ritardo T rispetto al suo-

no diretto, che produce dei notch spettrali con frequenza 1/T; siccome T

varia con l’elevazione, anche la frequenza dei notch varia. Possiamo quindi

ipotizzare che la forma spettrale e le frequenze dei notch siano i principali

segnali per la stima dell’elevazione di una sorgente sonora. Nel 1978 Watkins

formalizzo il modello proposto da Batteau in un single-delay-and-add system.

In Fig. 2.8 possiamo notare come al suono diretto si sommino due

riflessioni:

riflessione della conca: caratterizzata da un coefficiente di riflessione fisso

ρA e un tempo di ritardo fisso τA; quindi una riflessione sempre presente

6Si puo pensare la pinna montata su un piano infinito.

2.3. PRTF E ANTROPOMETRIA 43

e indipendente;

riflessione dell’elice: caratterizzata da un coefficiente di riflessione fisso ρV

e da un ritardo dipendente dall’elevazione τV (φ) = d/c, con d la diffe-

renza di cammino rispetto all’onda diretta e c la velocita del suono. La

formulazione single-delay-and-add fornisce una collezione di frequen-

ze di notch date dall’espressione: fn = n(1/2τ(φ)) con n=1,3,5... (v.

paragrafo successivo) .

Watkins dimostro che questo modello ha la capacita di approssimare bene

i dati reali per elevazioni all’interno dell’intervallo 0 − 45. La grande po-

tenzialita di questo modello e la sua semplicita, anche se i notch spettrali

prodotti si discostano abbastanza dai notch dei dati sperimentali e il coef-

ficiente di riflessione dovrebbe essere dipendente dalla frequenza7. E’ stato

largamente utilizzato in cascata all’uscita dell’HAT model[4] in modo da otte-

nere un modello strutturale di HRTF completo, semplice e allo stesso tempo

piuttosto efficace.

Lopez-Poveda e Meddis [26] estesero il modello di Batteau/Watkins in-

cludendo l’effetto di diffrazione per spiegare la rilevazione delle riflessioni da

parte della parete posteriore della conca per sorgenti nel piano mediano, sia

laterale che verticale. Utilizzarono misurazioni acustiche su un pezzo di la-

miera e calcoli numerici basati sulla teoria della diffrazione di Kirchhoff su

una forma cilindrica di conca. Attraverso varie approssimazioni studiarono

la pressione di diffrazione associato ad ogni punto della conca, considerando

come relazione di base per la funzione di trasferimento della conca:

pTp0

= 1 +pRp0

(2.29)

dove pT e la pressione totale all’entrata dell’apertura della conca, pR e la pres-

sione di riflessione della parete della conca misurata all’entrata dell’apertura

della conca, e p0 e la pressione del suono incidente. Al variare dell’elevazio-

ne e dell’azimuth Lopez-Poveda analizzarono il comportamento della HRTF

7Soprattutto quando la superficie riflettente e piccola rispetto alla lunghezza d’onda.


sperimentale in relazione al modello sviluppato. Trovarono tre notch princi-

pali: un notch principale associato alla profondita della conca, un altro meno

intenso associato alla crus helias e un ultimo (non individuato dal modello)

probabilmente associato al tragus. Tutti e tre i notch hanno un’evoluzione

della frequenza centrale in base alle variazioni di θ e φ.

2.3.1 Modelli Strutturali per PRTF

In questo paragrafo conclusivo trattiamo tre lavori che hanno maggiormente

guidato le idee su cui si basa questa tesi.

Barreto et al. [34][35][36] approfondisce la struttura risonanza piu

ritardo e cerca di relazionarne i risultati all’antropometria della pinna;

Raykar et al. [25][14] sviluppa un procedimento per l’individuazione dei

notch piu significativi e cerca di fornirne una mappatura sulla superficie

della pinna;

Satarzadeh et al. [21][22] focalizza la sua attenzione sul caso specifico

θ = 0 e φ = 0, fornendo un modello risonanza piu notch.

Barreto et al. ha proposto un modello per la pinna nel quale le trasforma-

zioni del suono che viaggia verso il timpano sono create dalla sovrapposizione

di un numero di riflessioni che vengono condizionate dal risonatore conca. Co-

me nel primo modello di Watkins/Batteau i cammini paralleli rappresentano

i multipli rimbalzi del suono sulla struttura geometrica della pinna.

Diverse traiettorie hanno differenti lunghezze, modellate dai delay τi. La

perdita di energia per ogni riflessione del suono e modellata dal fattore di

ampiezza ρi presente in ogni cammino del diagramma a blocchi. La risonan-

za puo venir rappresentata con una funzione di trasferimento del II ordine,

quindi da due parametri, l’angolo e il raggio dei poli nel dominio z. Questo

modello richiede in tutto la definizione di 9 parametri da relazionare all’an-

tropometria dell’orecchio esterno. Barreto procede scomponendo la HRIR


Figura 2.9: Diagramma a blocchi del modello di pinna risonanza piu ritardo.

nella sovrapposizione di 4 sinusoidi scalate e smorzate 8, in accordo al mo-

dello in Fig. 2.9. Ogni sinusoide sara scalata di un fattore ρi e ritardata

con latenza τi (τ1 = 0). Nel dominio z una singola componente sinusoidale

smorzata si puo modellare come un dipolo:

X(z) =kz2

(z − p1)(z − p2)(2.30)

dove k e uno scalare e p1 e p2 sono poli complessi coniugati. Il corrispondente

segnale nel dominio del tempo ha la forma:

xi = edin sin(ωdπn) (2.31)

con n = 1...N, N e lunghezza del segnale, di e il fattore negativo di smorza-

mento e ωd e la frequenza della sinusoide.

Vengono proposti due processi di decomposizione dell HRIR in compo-

nenti sinusoidali. Il primo [34] e eseguito da un algoritmo sequenziale che

racchiude una modellazione del secondo ordine tramite il metodo Steiglitz-

McBride (STMCB) su finestre consecutive definite a partire dalla HRIR mi-

surata. Lo scopo del processo sequenziale e quello di restringere l’analisi a

piu finestre parziali di dati di lunghezze diverse, in modo da trovare quella

dove solo una componente sinusoidale (o meglio la sua approssimazione del

8La risposta di un risonatore ad un impulso fornisce un segnale sinusoidale smorzato.


secondo ordine) e presente. L’altro algoritmo[35], procede in maniera diversa

e con un costo computazionale inferiore, ma risultati leggermente peggiori.

Inizia approssimando la HRIR con un modello 2*N poli utilizzando STMCB,

con N numero di componenti sinusoidali da ricercare; vengono calcolate le

sinusoidi smorzate candidate considerando ogni coppia di poli coniugati e

si sottraggono alla HRIR originale ottenendo N residui. Il processo si ite-

ra fino ad avere un modello del secondo ordine che fornira N! possibilita,

corrispondenti alle foglie dell’albero delle iterazioni. Si scegliera il cammi-

no che dalla radice (il modello di ordine 6) alle foglie approssima meglio la

HRIR9. La relazione con i parametri antropometrici e stata studiata median-

te regressione lineare su 8 parametri al variare dell’elevazione. Purtroppo la

documentazione riguardante questa fase del lavoro e poco chiara e le grandez-

ze antropometriche considerate sono state individuate tramite scanner 3D.

Le informazioni necessarie sono di difficile misurazione e la semplicita del

modello viene vanificata.

Raykar et al presenta un lavoro per estrarre automaticamente le frequenze

corrispondenti ai notch spettrali. La difficolta maggiore riguarda la separa-

zione di molti effetti combinati nella HRIR misurata, quali la diffrazione della

testa, la spalle, il busto e le riflessioni artificiali10 delle ginocchia. Vengono

sviluppate tecniche robuste di signal processing per individuare solamente in

notch dovuti alla pinna. Si utilizzano le seguenti tecniche: residuo di un mo-

dello LP (linear prediction), funzione di autocorrelazione finestrata, funzione

group-delay, modellazione tutti poli, guidate dalla conoscenza a priori della

9L’algoritmo prevede di non generare l’albero delle iterazioni completo; viene introdotta

una soglia sull’ampiezza del picco iniziale della sinusoide per la quale, al mancato supera-

mento, non viene calcolata la soluzione, interrompendo di fatto un ramo per le iterazioni

successive10Le misurazioni per ricavare i dati del CIPIC database sono state effettuate su soggetti

in posizione seduta. Si e potuto notare come nella HRIR apparisse una riflessione molto

ritardata (dopo 2 ms dal suono diretto), dovuta alle ginocchia. Questa osservazione non

si nota per suoni con elevazione φ > 90, dietro all’ascoltatore[25].


fisica del problema. Le sequenza di operazioni puo essere riassunta in questi

passi:

1. determinare l’onset iniziale della HRIR e usare la HRIR da quell’istante;

2. derivare il residuo trovato come differenza tra la HRIR del passo (1) e

il modello LP di ordine 10-12;

3. effettuare la finestratura del residuo usando meta di una finestra Hann

di 1.0 ms circa;

4. effettuare l’autocorrelazione del segnale ottenuto al passo (3);

5. effettuare la finestratura dell’autocorrelazione utilizzando una finestra

Hann di 1.0 ms circa;

6. calare la funzione di group-delay del segnale ottenuto al passo prece-

dente;

7. individuare i minimi locali utilizzando una soglia sulla funzione di

group-delay;

La LP analisi rimuove la parte ridondante del segnale, togliendo la parte pre-

vedibile. Le operazioni di finestratura sono eseguite per eliminare l’effetto

della spalle e del busto che si presenta dopo un tempo superiore al 1 ms. La

funzione di autocorrelazione di un segnale produce ampiezze decrescenti al

di fuori dei picchi (notch), questo comportamento favorisce il calcolo della

funzione di group-delay e nello stesso tempo preserva molte delle caratteristi-

che dell’inviluppo spettrale. La soglia utilizzata nella group delay e di -1, in

modo da poter eliminare i nulli spurii creati dall’operazione di finestratura.

Se x(t) e un’onda incidente, il segnale misurato all’ingresso del canale

uditivo (y(t)) sara la somma dell’onda diretta e di quella riflessa, come visto

in precedenza. Da cui otteniamo:

y(t) = x(t) + αx(t− τ(φ)) (2.32)


dove α e il coefficiente di riflessione e τ(φ) e il tempo di ritardo dato da:

τ(φ) =2d(φ)

c(2.33)

con 2d(φ) la distanza corrispondente al ritardo e c la velocita del suono. La

distanza d(φ) dipende dall’angolo φ e dalla forma della pinna. Il ritardo τ(φ)

causa nello spettro dei notch periodici della frequenza data da:

fn(φ) =(2n+ 1)

2τ(φ)=c(2n+ 1)

4d(φ), n = 0, 1... (2.34)

La frequenza del primo notch nello spettro si ha per n = 0 :

f0(φ) =c

4d(φ)(2.35)

In pratica ci sono una moltitudine di riflessioni all’interno della pinna,

ogni riflessione genera una serie di notch periodici. Al variare di dell’angolo

φ la frequenza del notch cambia in accordo con la forma della pinna: avendo a

disposizione la frequenza f0(φ) si puo ricavare la distanza d(φ). Hanno notato

che il primo notch appare a causa della riflessione da parte della conca e che

al variare dell’elevazione ne viene tracciata la forma. Il terzo notch potrebbe

essere causato dal crus helias che divide la conca in due parti.

Questo lavoro fornisce un’idea molto attraente: prendere un’immagine

dell’orecchio esterno e con tecniche di computer graphics estrarne le caratte-

ristiche antropometriche che genereranno i corrispondenti notch.

Il lavoro di Satarzadeh e centrato sulla HRTF per sorgenti posizionate di-

rettamente di fronte all’ascoltatore (θ = 0 e φ = 0) e cerca di trovare una

soluzione piu semplice possibile che possa racchiudere le caratteristiche piu

importanti dal punto di vista percettivo. Il punto di partenza e rappresen-

tato dagli studi di Shaw (v. par. 1.2.4) in cui viene proposto un modello di

conca rettangolare con alcune modifiche (chiamato model E ), Fig 2.10:

per ridurre l’eccitamento verticale della frequenza risonante e stata in-

trodotta una barriera parziale, che rappresenta la crus helais, che e

determinata dal valore F;


Figura 2.10: Modello finale di Shaw, denominato E model.

il terzo modo risonante viene rappresentato introducendo un canale

simile ad una fossa, che e determinato dai parametri E e C.

I parametri A,B e D che rappresentano l’altezza totale della conca, F rap-

presenta la larghezza e C e E rappresentano la lunghezza e il diametro della

fossa. Infine il parametro G rappresenta la profondita della conca.

Lo studio del modello e stato condotto utilizzando della simulazioni con

il metodo FDTD (v. par. 2.1.3) per elevazioni che variano da −30 a 30 ad

intervalli di 5. I risultati ottenuti sono stati usati in relazione ai parametri

seguenti:

F1 frequenza del primo picco;

G1 guadagno del primo picco;

F3dbl frequenza del primo punto con threshold di 3db;

Fn frequenza del primo notch;

Gn differenza tra il guadagno del picco e del notch;

F2 frequenza del secondo picco;


G2 guadagno del secondo picco;

Per ogni insieme di possibili input (A-G,φ) sono stati estratti gli output

(F1, G1, F3dbl, Fn, Gn, F2 e G2) e si sono analizzate le relazioni tra i parame-

tri di input e le loro risposte. La piu grande deviazione standard per ogni

output si ottiene per i parametri di input F e G. Questo indica che la lar-

ghezza e la profondita della conca dominano l’effetto degli altri parametri

antropometrici.

Dopo questa fase preliminare di analisi, Sataradeh propone dei sistemi

automatici per l’identificazione delle caratteristiche spettrali della PRTF:

Group Delay Analysis del cepstrum per l’individuazione dei notch

piu significativi; viene descritto come metodo robusto rispetto all’ope-

razione di finestratura che produce degli spostamenti degli zeri nelle

PRTF;

Scale-Space Filtering per individuare le risonanze della PRTF;

Il primo viene utilizzato per l’identificazione degli zeri dell PRTF, mentre

il secondo per i poli.

Riportiamo brevemente le due definizioni che stanno alla base dei due

metodi. La funzione di group delay e definita come la derivata negativa della

fase di una f.d.t rispetto alla frequenza angolare:

τd =d

dwargH(ejw) (2.36)

La rappresentazione scalata nello spazio, L, di un qualsiasi segnale ψ si ottiene

da:

L(µ, σ) = ψ(x) ∗ 1

2πσ2e−(x−µ)2

2σ2 (2.37)

La convoluzione con una Gaussiana con deviazione standard σ e media µ sop-

primera le caratteristiche di ψ(x) con una lunghezza minore di σ 11. Viene

11Le caratteristiche per un segnale monodimensionale sono gli estremi locali.


Figura 2.11: Modello strutturale proposto da Satarzadeh.

cosı proposto un modello strutturale in cui i due modi risonanti principa-

li, associati alla profondita e alle larghezza della conca, sono in parallelo e

connessi in serie ad un comb filter che cerca di tenere conto dei notch piu

importanti:

Hmodel(s) = (Hdepth(s) +Hwidth(s))Hcomb(s) (2.38)

H(s) =

(K1(s+ ω/3)

(s+ σ1 − jω1)(s+ σ1 + jω1)+

K3(s+ ω/4)

(s+ σ2 − jω2)(s+ σ2 + jω2)

)(1 + Γe−jwτ )

(2.39)

dove Γ e τ rappresentano il coefficiente di riflessione e il tempo di ritardo del

comb filter. Le risonanze hanno la forma:

H(s) =K1(s+K2)

s2 + ωbs+ (ωb/2)2 + ω2c

(2.40)

dove ωc determina la frequenza di risonanza, ωb la banda passante e K1 e K2

sono rispettivamente il guadagno della risonanza e la risposta del filtro alla

DC.

L’ultima parte del lavoro di Satarzadeh et al. riguarda la possibilita

di associare i pochi parametri del modello proposto con l’antropometria

dell’orecchio esterno. I parametri trattati sono:

Time Delay τ : si vuole utilizzare il tempo di ritardo per identificare la

superficie riflettente, che nel caso di questo studio e un solo punto che


causa dei notch periodici (comb filter) nello spettro della PRTF. Sia

w la distanza tra il tragus e la superficie riflettente. L’onda che passa

il tragus si riflette sulla superficie riflettente e ritorna all’entrata del

canale uditivo percorrendo una distanza di 2w da cui:

w =cτ

2(2.41)

Nel caso di superficie riflettente si avra un coefficiente di riflessione

ρ > 0 e se 2w e meta della lunghezza d’onda (o un multiplo dispari di

meta lunghezza d’onda) verra prodotto un notch spettrale:

2w = cτ =λ

2=

c

2f1

(2.42)

τ =1

2f1

(2.43)

con f1 la frequenza del primo notch.

Nel caso di ρ < 0 l’onda passa oltre il rim della pinna e incontra la bassa

impedenza dello spazio dietro la conca. In questo caso l’interferenza

distruttiva appare se 2w = λ da cui:

τ =1

f1

(2.44)

Traducendo queste formule sull’antropometria dell’orecchio esterno pos-

siamo osservare come la distanza calcolata per il punto di riflessione sia

o la parte posteriore della conca o il contorno dell’elice. Questo proba-

bilmente dipende dal flare angle, cioe l’angolo tra il piano che contiene

la pinna e il piano tangente alla testa nel punto di attacco dell’orecchio.

Depth Resonance ω1: viene utilizzato un modello cilindrico di conca. La

caratteristica risonante di un cilindro e conosciuta e la risonanza fon-

damentale associata alla profondita e data da:

λmax4

= d+ 0.822r (2.45)


dove d e la profondita del cilindro e r il raggio. Il termine 0.822r

rappresenta la correzione per un cilindro vuoto dove la risonanza di

profondita e influenzata dalla larghezza del cilindro. Da cui la frequenza

massima di risonanza:

fmax =c

4(d+ 0.822r)(2.46)

Width Resonance ω2: la risonanza dovuta alla larghezza della conca vie-

ne associata, per motivi fisici, al tempo di ritardo τ e in un secondo

momento all’antropometria. Da osservazioni empiriche la relazione tra

ω2 e τ dipende da quale, tra la conca o l’elice sia il principale riflettore.

A seconda dei due casi viene preso in considerazione il III o il IV notch

del comb filter:

ωn =1 + 2nπ

2τ, n = 1, 2, 3... (2.47)

Da cui la frequenza di risonanza fwidth:

fwidth =

72τ

n = 3, conca92τ

n = 4, elice(2.48)

Da questi tre lavori analizzati emergono ancora piu chiare le due componenti

fondamentali che stanno alla base di una modellazione di PRTF:

le risonanze: l’orecchio esterno ha delle proprie caratteristica come

corpo risonatore;

i notch: l’orecchio esterno utilizza la propria personale geometria per

fornire dei segnali spaziali;

Capitolo 3

Metodologie e strumenti

Questo capitolo si pone come agile aiuto per la descrizione delle metodologie e

degli strumenti utilizzati in questa tesi per la fase di analisi e di modellazione.

Per ogni strumento verranno presentate le caratteristiche piu importanti in

relazione allo scopo per cui e stato impiegato.

Lo poniamo ad introduzione dei capitoli centrali della tesi, come supporto

per la buona comprensione e a completamento dell’esposizione. In certi casi,

presentando argomenti consolidati, ci limiteremo a semplici richiami teorici,

rintracciabili facilmente in letteratura.

3.1 Cepstrum

Il metodo del cepstrum (inverso di spec-trum) permette di stimare l’inviluppo

spettrale di un segnale finestrato x(n), a partire dalla sua trasformata discreta

di Fourier X(k) (con k indice delle frequenze discrete).

X(k) =N−1∑n=0

x(n)W knN = |X(k)|ejϕx(k), k = 0, 1, ..., N − 1 (3.1)

Da X(k) viene preso il logaritmo

X(k) = logX(k) = log |X(k)|+ jϕx(k) (3.2)

55

56 CAPITOLO 3. METODOLOGIE E STRUMENTI

ed effettuata la IFFT di X(k), che genera il cepstrum complesso:

x(n) =1

N

N−1∑k=0

X(k)W−knN (3.3)

Ricavata la parte reale di X(k),

XR(k) = log |X(k)|; (3.4)

ed eseguendone la IFFT si ottiene il real cepstrum:

c(n) =1

N

N−1∑k=0

XR(k)W−knN (3.5)

Quindi il real cepstrum c(n) e la IFFT della log magnitude della FFT di

un segnale finestrato x(n). Per generare un cepstrum pesato, la finestra

con caratteristiche low-pass si presenta nella seguente forma:

wLP (n) =

1 n = 0, N1

2 1 ≤ n < N1

20 N1 < n ≤ N − 1

, N1 ≤ N/2 (3.6)

La FFT di un cepstrum finestrato cLP (n) e CLP = FFT [cLP (n)] e rappresen-

ta la versione smussata dello spettro X(k) in db. N1 e il numero di coefficienti

cepstrali che contribuiscono alla determinazione dell’inviluppo. L’avvicinarsi

di N1 a N produce un risultato sempre piu simile a X(K).

Il cepstrum e anche utilizzato per estrarre lo spettro stazionario. Nei

modelli source-filter permette la separazione di un segnale y(n) = x(n) ∗h(n) nel segnale sorgente x(n) e nella risposta impulsiva h(n) del filtro. La

corrispondente relazione nel dominio della trasformata di Fourier diventa

Y (ejΩ) = X(ejΩ) · H(ejΩ). Riscrivendo H(ejΩ) = |H(ejΩ)|e−6 H(ejΩ), si puo

ipotizzare fermamente che la risposta in frequenza del filtro sia a valori reali

e che la fase sia assegnata al segnale di sorgente.

3.1. CEPSTRUM 57

Il cepstrum reale fornisce una stima dell’inviluppo spettrale in relazione

all’ampiezza del segnale, utilizzando la proprieta matematica del logaritmo

(log(a · b) = log(a) + log(b)):

|Y (ejΩ)| = |X(ejΩ)| · |H(ejΩ)|log|Y (ejΩ)| = log|X(ejΩ)|+ log|H(ejΩ)|

(3.7)

Separare log|H(ejΩ)| da log|X(ejΩ)| significa separare le variazioni lente del

filtro, da quelle rapide della sorgente. Se consideriamo il diagramma di

log|Y (ejΩ)| come un segnale nel dominio del tempo, possiamo distinguere

due componenti: un’oscillazione veloce, dovuta alla struttura armonica, e un

andamento piu lento in relazione alle risonanze del filtro (inviluppo spettrale).

La separazione della sorgente dal filtro puo essere attuata pesando di

applicare al cepstrum (c(n)) due finestre: una lowpass window (wLP (n)) e

una highpass window (wHP (n)):

c(n) = cx(n) + ch(n) (3.8)

cx(n) = c(n) · wLP (n)

ch(n) = c(n) · wHP (n)

c(n) = c(n) · wLP (n) + c(n) · wHP (n)

(3.9)

I valori di tempo bassi (basse quefrencies) ricavati dal filtraggio passa basso

restituiscono l’inviluppo spettrale in db, log|H(ejΩ)|, e i valori di tmpo alti

(alte quefrencies) la stima della sorgente log|X(ejΩ)|. Per separare le due

componenti si puo usare un metodo basato sulla IFFT:

IFFT (log|Y (ejΩ)|) = IFFT (log|X(ejΩ)|) + IFFT (log|H(ejΩ)|) (3.10)

La soglia tra le basse e le alte quefrencies e fondamentale per un giusto

equilibrio tra sorgente e filtro. Abbassandola si otterra un inviluppo spettrale

piu smussato, mentre alzandola si includeranno alcune armoniche o parziali

in piu nell’inviluppo.


Figura 3.1: Diagramma riassuntivo per il calcolo del cepstrum.

3.2 Peak Filter

La funzione di trasferimento di un filtro shelving del primo ordine a basse

frequenze puo essere decomposta come [41]:

H(s) =s+ V0ωcs+ ωc

= 1 +H0ωc

s+ ωc

= 1 +H0

2

[1− s− ωc

s+ ωc

] (3.11)

3.2. PEAK FILTER 59

con

V0 = H(s = 0), (3.12)

H0 = V0 − 1, (3.13)

V0 = 10G/20(G in dB), (3.14)

La funzione di trasferimento e composta da un cammino diretto e da un filtro

passa basso. Il filtro passa basso del primo ordine e implementato utilizzando

una decomposizione passa tutto. Applicando la trasformata bilineare a (3.11)

si ottiene:

H(z) = 1 +H0

2[1− A(z)] (3.15)

con

A(z) = − z−1 + aB1 + aBz−1

(3.16)

Ripercorriamo lo stesso procedimento per un filtro shelving che taglia le basse

frequenze, considerando la seguente decomposizione:

H(s) =s+ ωc

s+ ωc/V0

= 1 + (V0 − 1)︸︷︷︸H0

ωc/V0

s+ ωc/V0

= 1 +H0

2

[1− s− ωc/V0

s+ ωc/V0

] (3.17)

La trasformazione bilineare applicata a (3.17) ritorna la (3.15). La struttura

del filtro e identica sia nel caso di boost sia di cut delle basse frequenze. I

parametri aB per il boost e aC per il cut hanno le seguenti forme:

aB =tan(ωcT/2)− 1

tan(ωcT/2) + 1(3.18)

aC =tan(ωcT/2)− V0

tan(ωcT/2) + V0

(3.19)

con ωc pulsazione di taglio. Assemblando assieme (3.15) e (3.16) con la

distinzione di aB e aC , appena definita, si ottiene la funzione di trasferimento


Figura 3.2: Filtro shelving per le basse frequenze e filtro passa basso,

entrambi del primo ordine.

di un filtro shelving per le basse frequenze del primo ordine nel dominio

discreto:

H(z) =1 + (1 + aBC)H0

2+ (aBC + (1 + aBC)H0

2)z−1

1 + aBCz−1(3.20)

Con A1(z) = −A(z) il diagramma in Fig. 3.2 descrive un filtro passa basso

del primo ordine e un filtro shelving per la basse frequenze del primo ordine.

Un peak filter del secondo ordine puo essere ricavato dalla trasforma-

zione di un filtro shelving del primo ordine da passa basso a passa banda.

L’aggiunta di un filtro passa banda del secondo ordine oltre al cammino co-

stante diretto (come nel caso precedente) portano come risultato un peck

filter.

Utilizziamo una implementazione passa tutto per creare un filtro passa

banda:

H(z) =1

2[1− A2(z)] (3.21)

A2(z) =−aB + (d− daB)z−1 + z−2

1 + (d− daB)z−1 − aBz−2(3.22)

3.3. MODELLI SINUSOIDALI 61

da cui un peak filter del secondo ordine puo essere espresso similmente a

(3.15)

H(z) = 1 +H0

2[1− A2(z)] (3.23)

I parametri per la banda passante sono aB o aC a seconda che si tratti di un

picco (boost) o di un notch1 (cut).

aB =tan(ωbT/2)− 1

tan(ωbT/2) + 1(3.24)

aC =tan(ωbT/2)− V0

tan(ωbT/2) + V0

(3.25)

con ωb pulsazione di banda passante ai 3 dB. Il parametro che controlla la

frequenza centrale e d, e H0 controlla il guadagno.

d = − cos(ωc) (3.26)

V0 = H(ejωc), (3.27)

H0 = V0 − 1, (3.28)

La funzione di trasferimento di un peak(notch) filter a guadagno controllato

del secondo ordine ha equazione 3.29 e diagramma rappresentato in Fig. 3.3.

H(z) =1 + (1 + aBC)H0

2+ d(1− aBC)z−1 + (−aBC − (1 + aBC)H0

2)z−2

1 + d(1− aBC)z−1 − aBCz−2

(3.29)

3.3 Modelli Sinusoidali

La modellazione sinusoidale trae origine dalla tecnica di sintesi additiva. L’i-

dea guida deriva dal teorema di Fourier cha afferma la possibilita di modellare

una qualsiasi funzione periodica, per mezzo di una somma di sinusoidi di varia

ampiezza e dalle frequenze armoniche. Per un suono stazionario, le ampiezze

e le frequenze evolvono lentamente con il tempo in maniera continua, dando

1Si intende un notch a guadagno controllato


Figura 3.3: Peak (Notch) filter e filtro passa banda, entrambi del secondo

ordine.

origine ad un insieme di oscillatori pseudo-sinusoidali chiamati comunemente

parziali.

La modellazione sinusoidale descrive un segnale audio x(t) come somma

di P parziali.

x(t) =P∑p=1

Ap(t) sin(φp(t)) (3.30)

φp(t) = φp(0) + 2π

∫ t

0

fp(u)du (3.31)

dove Ap(t) e l’ampiezza modulata e φp(t) e la fase modulata della parziale

p. L’equazione (3.30) viene comunemente espressa nel dominio discreto in

questo modo:

x(n) =P∑p=1

Ap(n) sin(φP (n)) (3.32)

Ap(n) ≈ Ap (3.33)

φp(n) ≈ Ωpn+ φp(0) (3.34)

con Ap e Ωp considerati costanti all’interno di ogni frame di analisi, preso

sufficientemente piccolo.


Figura 3.4: Porzione di un sistema di analisi sinusoidale.

Lo scopo principale dell’analisi sinusoidale e la stima degli Ap e Ωp lungo

il segnale diviso in frame consecutivi. Una porzione di questo procedimento

e illustrata in Fig. 3.4.

La fase di Time/Frequency Decomposition calcola la STFT (Short-Time

Fourier Transform) del segnale audio x(n):

X(m, k) = STFT (x(n,m)) =1

N

N−1∑n=0

w(n)x(n+mH)e−jk2πNn (3.35)

dove w(n) e la funzione finestra di lunghezza N (es. una Hamming window),

k l’indice discreto delle frequenze, m l’indice di frame, H il salto in campioni

tra un frame e l’altro. A seguire la fase di peak detection (par. 3.3.1) seleziona

i picchi che corrispondono a sinusoidi stazionarie presenti nel frame m. Infine

il partial tracking (par. 3.3.2) raggruppa i picchi lungo i frame consecutivi,

creando le tracce per le parziali.

3.3.1 Peak Detection

Il modello sinusoidale assume che ogni spettro della rappresentazione STFT

possa essere espresso da un serie di sinusoidi. Per una risoluzione frequen-

ziale che usi abbastanza punti, una sinusoide puo essere descritta dalla sua

forma. Teoricamente una sinusoide che e stabile in ampiezza e frequenza ha

una rappresentazione ben definita. In pratica raramente le parziali sono per-

fettamente periodiche, ben spaziate e i picchi sono chiaramente definiti nel

dominio della frequenza. Le interazioni tra le diverse componenti rendono

complicata l’individuazione della forma dei picchi.


Una soluzione pratica e quella di trovare piu picchi possibili, vincolando

la ricerca con pochi fattori, e ritardando la decisione di quali siano le vere

parziali al prossimo passo d’analisi (partial tracking).

Definiamo un picco come un massimo locale nella risposta in ampiezza

di un spettro. Gli unici vincoli imposti sono il range di frequenza e la soglia

di ampiezza. A causa della natura campionata dello spettro, ogni picco avra

l’accuratezza di mezzo campione 2. La tecnica di zero-padding nel dominio

del tempo aumenta il numero di campioni spettrali per Hz e quindi aumenta

l’accuratezza di un semplici algoritmo di peak detection. Puo venir utiliz-

zato uno schema di interpolazione piu efficace (quadratico, spline, etc.), in

modo da investigare i campioni attorno al punto di massimo e raggiungere

l’accuratezza dello 0.1 %.

Di uguale interesse per questa tesi e lo scenario speculare che coinvolge i

minimi locali (notch) per cui possono venire utilizzate le medesime tecniche.

3.3.2 Partial Tracking

Il tracking delle parziali e stato proposto da McAulay e Quatieri per l’analisi

del parlato nella rappresentazione sinusoidale [15] nel 1986 ed ha avuto una

immediata applicazione per molti altri tipi di segnali. Vi presentiamo l’algo-

ritmo originale e alcuni miglioramenti successivamente proposti in letteratura

per aumentarne l’efficacia.

L’analisi di un segnale tempo variante viene effettuata per frame con-

secutivi, in altre parole una finestra scorre il segnale determinando il focus

dell’algoritmo. Se il numero dei picchi fosse costante da frame a frame,

l’associazione di un picco con il suo successivo si ridurrebbe ad un banale

ordinamento in frequenza dei picchi. In pratica i picchi compaiono e scom-

paiono a causa dei complicati fenomeni acustici presenti nel segnale. Nel

caso originale del parlato la posizione dei picchi cambia con il pitch e con la

2Uno spettro campionato viene rappresentato ad intervalli di fs/N Hz, dove fs e la

frequenza di campionamento e N e la taglia della FFT.


Figura 3.5: McAulay - Quartieri algoritmo di partial tracking.

rapida variazione delle parti vocalizzate o meno. Si puo definire il concetto

di nascita e morte di componenti sinusoidali.

Supponiamo che nel frame k-esimo siano selezionati N picchi con frequen-

ze ωk0 , ωk1 , ..., ω

kN−1 e nel frame k+1 M picchi con frequenze ωk+1

0 , ωk+11 , ...ωk+1

M−1

(in generale N 6= M) ordinate in maniera crescente all’interno del proprio

frame. Il processo che permette di associare ogni frequenze del frame k, ωkn

(0 ≤ n < N), con qualche frequenza del frame k+1, ωkm (0 ≤ m < M), e

descritta da tre passi:

1. si supponga che sia stato trovato un match per le frequenze ωk0 , ωk1 , ..., ω

kn−1

e si cerca di trovarlo per ωkn. La Fig. 3.5.a rappresenta il caso in

cui tutte le frequenze del frame k + 1 stanno fuori l’intervallo di


corrispondenza ∆ di ωkn

|ωkn − ωk+1m | ≥ ∆, per ogni m (3.36)

In questo caso la track associata alla frequenza ωkn viene dichiarata

morta: le viene assegnata ampiezza nulla nel frame k + 1 e non viene

piu considerata nella prosecuzione dell’algoritmo. Questo primo passo

viene ripetuto per le tutte le frequenze, continuando da ωkn+1.

Se esiste ωk+1m nel frame k+ 1 che sta all’interno di ∆ ed e la frequenza

piu vicina a ωkn:

|ωkn − ωk+1m | < |ωkn − ωk+1

i | < ∆ (3.37)

con tutti gli i 6= m. Il match si potra dire definitivo solamente se

risultera il migliore alla fine dell’esecuzione dei tre passi.

2. un tentativo ha associato ωkn con ωk+1m ; se ωk+1

m non ottiene candidati

migliori quello che gli e stato associato viene considerato definitivo.

Questo caso e illustrato in Fig. 3.5.c in cui

|ωk+1m − ωkn| < |ωk+1

m − ωki+1|per i > n (3.38)

l’analisi per ωkn e ωk+1m si conclude per proseguire con ωkn+1. Se questa

condizione non e verificata ωk+1m ha un match migliore con ωkn+1 e si

creano due sottocasi illustrati in Fig 3.5.d e 3.5.e . Nel primo sottocaso

(d) la frequenza minore adiacente a ωk+1m , ωk+1

m−1, non rientra nell’inter-

vallo di corrispondenza di ωkn. In questo caso la traccia associata a ωkn

viene dichiarata morta all’entrata del frame k + 1. Nel secondo sotto-

caso (e) la frequenza ωk+1m−1 giace all’interno di ∆ di ωkn e viene definito

il match. In entrambi i casi viene ripetuto il passo 1 usando ωkn+1 come

frequenza da analizzare.

3. quando tutti i picchi del frame k sono stati assegnati a tracce vive

o morte, nel frame k + 1 potrebbero rimanere delle frequenze libere.

Ipotizziamo che ωk+1m sia una di queste, possiamo dire che e nata una

nuova possibile parziale nel frame k di ampiezza zero.

3.4. CIPIC DATABASE 67

Un semplice miglioramento per questo algoritmo introduce due parametri

in aggiunta a ∆: E,S.

E: se una parziale appena nata non trova una continuazione dopo E frame,

viene scartata; questo permette di rimuovere le tracce corte, che po-

trebbero esser state erroneamente identificate in una fase precedente

all’applicazione dell’algoritmo.

S: se una parziale non trova una continuazione nel frame k + 1 ha a disposi-

zione S frame per trovarla, dopodiche viene scartata; questo permette

di diminuire le discontinuita tra picchi erroneamente non individuati.

In letteratura sono proposte moltissime altre varianti di questo famoso algo-

ritmo [16][17][18] applicate a diverse tipologie di segnale. In questo lavoro di

tesi ne viene implementata la versione base, considerando i frame non tem-

poralmente consecutivi, ma spazialmente adiacenti. Il procedimento resta il

medesimo.

3.4 CIPIC Database

Il CIPIC Interface Laboratory di U.C. Davis ha misurato le HRTF con un

alta risoluzione spaziale per 90 soggetti. Sono stati resi di dominio pubblico

(all’indirizzo http://interface.cipic.ucdavis.edu) i dati riguardanti 45 sogget-

ti: 43 umani piu due manichini KEMAR (Knowles Electronics Manikin for

Acoustic Research3)con pinna piccola o grande. In aggiunta alle 1250 rispo-

ste all’impulso per le varie direzioni, il database contiene le informazioni di

alcune grandezze antropometriche per ogni soggetto. Vengono inoltre fornite

delle utility in MatlabTM per navigare, in modo ordinato, nel database.

Questa tesi utilizza le HRTF per sorgenti che provengono direttamente

dalla direzione frontale rispetto all’ascoltatore. Se la HRIR e finestrata per

rimuovere le riflessioni delle spalle che arrivano con ritardo TS, per frequenze

al di sopra di 1/TS il risultato e essenzialmente lo stesso di una pinna isolata.

3Rappresentazione Antropomorfa dell’adulto medio, sotto forma di manichino.

http://interface.cipic.ucdavis.edu


Figura 3.6: Posizione dei punti di ripresa, veduta frontale (a) e laterale (b).

Questo ci permette di utilizzare i dati del CIPIC HRTF database per studiare

gli effetti della pinna soltanto (PRTF).

3.4.1 Misurazioni

Tutte le misurazioni sono state effettuate su soggetti seduti al centro di un

cerchio di raggio 1 m allineato con l’asse interaurale. La posizione della testa e

lasciata svincolata, ma il soggetto deve mantenerla approssimativamente fissa

per non creare discontinuita spettrali. Lungo il cerchio sono stati montati

degli altoparlanti Bose AcoustimassTM dal cono di 5.8 cm di diametro per

l’emissione di segnali Golay-code. I due canali uditivi del soggetto sono

stato bloccati e sono stati posti all’ingresso dei microfoni Etymotic Research

ER-7C per raccogliere il segnale filtrato. L’output dei microfoni e stato

digitalizzato a 44.1 kHz con una risoluzione di 16 bit. Una finestra Hamming

e stata applicata alla HRIR grezza per rimuovere le riflessioni della stanza e

il risultato e stato compensato con una ripresa free-field4 per correggere le

caratteristiche spettrali dei trasduttori. La lunghezza di ogni HRIR e di 200

campioni, corrispondenti alla durate di 4.5 ms. La posizione della sorgente

sonora e descritta dalle coordinate polari interaurali, angolo di azimuth e di

elevazione (θ, φ). L’elevazione e stata uniformemente campionata con passo

4La risposta free-field e stata misurata nella posizione corrispondente al centro della

testa, in assenza del soggetto.


360/64=5.625 da -45 a +230.625. Per ottenere approssimativamente la

densita uniforme della sfera, l’azimuth e stato campionato a -80, -65, -55

e da -45 a 45 con passo di 5, 55, 65 e 80.

3.4.2 Antropometria

Il CIPIC HRTF database include un insieme di misurazioni antropometriche

che possono venir utilizzate per vari scopi: nel nostro caso aiutare a trovare

una corrispondenza tra l’HRIR e la forma della pinna.

L’insieme e le specifiche di un insieme di misurazioni antropometriche

definite e rilevanti e ancora un problema irrisolto. Soprattutto per quanto ri-

guarda la pinna, dove piccole variazioni possono produrre grossi cambiamenti

all’andamento della HRTF.

Le grandezze misurate traggono spunto dal lavoro di Genuit (v. par. 2.2).

Sono definite 24 misurazioni antropometriche: 17 per la testa e per ii busto

(Fig. 3.8) e 10 per la pinna (Fig. 3.9). Altri criteri che hanno guidato la

scelta di questi parametri sono:

la caratteristica antropometrica sembra avere un’influenza significativa

sulla HRTF;

la caratteristica antropometrica e facilmente misurabile;

Descriviamo le modalita di misurazione dei vari parametri:

la maggior parte dei parametri e ricavabile da foto digitali ad alta

definizione;

i parametri di altezza e circonferenza x14, x15, x16, x17 sono stati misu-

rati con un metro a nastro.

la profondita della conca, d8, e l’angolo di attacco della pinna5, θ2, sono

stati misurati con Polhemus 3D stylus digitizer ;

5L’angolo d’attacco della pinna e definito in riferimento al piano tangente alla testa,

determinato dall’acquisizione 3D.


Figura 3.7: Parametri antropometrici.


per ogni soggetto vengono anche riportati anche sesso, eta e peso cor-

poreo.

In [3] e stato condotto uno studio sull’intero database per analizzare la pre-

senza di relazioni tra i vari parametri. In generale sono stati rilevati risultati

statistici interessanti, ma con correlazioni deboli tra le coppie di misurazioni.

Figura 3.8: Misurazioni testa, busto e spalla.


Figura 3.9: Misurazioni orecchio esterno.

Capitolo 4

Analisi della PRTF

Questo capitolo e dedicato alla descrizione e alla discussione di una prima ver-

sione di algoritmo per la separazione del contributo dei notch dal contributo

delle risonanze in una PRTF personale. Queste informazioni sono fondamen-

tali in fase di analisi per poter comprendere a fondo il giusto contributo dei

fenomeni fisici che coinvolgono l’orecchio esterno.

Per iniziare il nostro lavoro di analisi su alcuni individui del CIPIC da-

tabase [3], abbiamo prestato molta attenzione allo studio approfondito dei

notch e delle risonanze effettuati da Raykar [25] e Satarzadeh [22] (esposti

nel par. 2.3.1). In particolare l’interesse per i modelli strutturali, come so-

luzione costruita ad hoc per il problema della localizzazione del suono, ci ha

indirizzato verso una scomposizione netta tra i due fenomeni. In questa breve

introduzione vorremmo soffermarci sui pro e i contro dei due lavori proposti,

in modo da giustificare alcune scelte guida nella progettazione dell’algoritmo.

Il lavoro di Raykar ha evidenziato l’importanza fondamentale dei notch

e la loro relazione con la forma della pinna, confermando il modello single-

delay-and-add gia proposto da Watkins. Tuttavia, tenendo presente cosı

tanti notch, il modello risulatante viene appesantito dalla considerazione di

un numero elevato di riflessioni che producono a loro volta una serie comb di

notch.

Il lavoro di Satarzadeh e molto prezioso per la ricerca di un modello sem-

73

74 CAPITOLO 4. ANALISI DELLA PRTF

plice ed e progettato tenendo ben presente la conoscenza della percezione acu-

stica dell’uomo. Soprattutto e molto completo nella descrizione riguardante

la scelta delle risonanze dominanti e della loro relazione con l’antropometria

dell’orecchio, avvalorando gli studi di Lopez-Poveda/Meddis[26]. Purtroppo,

oltre a limitarsi alla considerazione della sola direzione frontale del suono,

cerca di ricondurre i notch spettrali ad un’unica riflessione principale, non

trovando dei risultati uniformi.

La nostra analisi ha come punto di partenza e di arrivo la struttura di

un filtro multi-notch essenziale, composto dalla somma dei notch della

PRTF percettivamente importanti. Il processo di calcolo inizia con una stima

rozza di notch filter multiplo per arrivare ad un filtro-risultato realistico in

frequenze, guadagni e bande passanti. Dal risultato finale si puo ricavare il

contributo delle risonanze come sottrazione dalla PRTF. La modellazione dei

risultati verra trattata nel capitolo conclusivo.

4.1 Algoritmo di analisi

L’algoritmo proposto per l’analisi dei dati delle HRTF e di tipo iterativo e

utilizza tecniche DSP (Digital Signal Processing). La HRTF presa in consi-

derazione e ottenuta dalla trasformata di Fourier della HRIR finestrata con

una Hann window della lunghezza di 1 ms, in modo da escludere gli effetti

delle riflessioni ad opera del busto e delle spalle. Inoltre prendiamo in con-

siderazione il piano mediano (θ = 0), variando l’elevazione da −45 a 90;

quindi, in questa situazione, gli effetti dovuti alla testa sono trascurabili. La

condizioni di partenza sono le seguenti:

HRTF finestrata, che chiamiamo PRTF (per le considerazioni sopra

esposte).

Filtro multi-notch con spettro piatto.

L’idea di questo algoritmo e che ad ogni iterazione la PRTF verra compensata

con un filtro multi-notch approssimato che andra a sommare il suo contributo

4.1. ALGORITMO DI ANALISI 75

ai filtri creati nelle iterazioni precedenti; quando si arrivera alla convergenza,

la PRTF conterra solamente le risonanze, mentre la somma dei filtri dara il

multi-notch filter completo.

Cerchiamo di descrivere formalmente il funzionamento dell’algoritmo, espri-

mendolo in pseudo-codice, per poi soffermarci sulle scelte di progetto effet-

tuate. Possiamo vedere l’algoritmo come suddiviso nelle seguenti fasi:

Algorithm 4.1 Algoritmo di analisi, Fasi

1 - Inizializzazione

while ci sono notch sotto una certa soglia do

2 - Calcolo del residuo

3 - Ricerca parametri per filtro multi-notch

4 - Costruzione filtro multi-notch

5 - Aggiornamento dello stato

end while

Per comprendere e analizzare il comportamento dell’algoritmo, descrivia-

mo ogni fase facendo riferimento allo pseudocodice di Alg 4.2 e riportando,

quando necessario, le righe di codice MatlabTM sviluppate (in Appendice B

si trovano tutti i sorgenti delle routine implementate).

4.1.1 Inizializzazione

Nella fase di inizializzazione vengono istanziate le condizioni iniziali dell’algo-

ritmo e definiti parametri che influenzano in maniera determinante il risultato

finale.

Come esposto ad inizio capitolo le condizioni iniziali all’ingresso del ciclo

iterativo sono:

lo spettro da analizzare e quello della HRTF finestrata per sorgente

sonora sul piano mediano;

lo spettro del filtro multi-notch e piatto;


Algorithm 4.2 Algoritmo di analisi, Descrizione formale

1: PRTF0 ← HRTFwindowed; —Condizioni iniziali—2: Hmulti−notch ← flat spectrum;

3: Nceps,maxit,soglia,rid; —Parametri dell’algoritmo—4: it← 0;

5: nn;

6: while nn > 0 and it ≤ maxit do

7: CEPSi ← Re[CEPSTRUM] di PRTFi con Nceps coefficienti cepstrali;

8: |RESi| ← |PRTFi| − CEPSi; in dB9: ERESi ← |RESi|2;

10: notch← PickNotch(|PRTFi|);11: fc, Amp, fb; —Ricerca parametri per filtro multi-notch—12: A, B; —Costruzione filtro multi-notch—13: nn← 0;

14: for j = 1 to length(notch) do

15: Amp← |RESi(notch(j))|;16: if Amp ≥ soglia then

17: fc← notch(j);

18: fb←ricerca su ERESi ;

19: (b, a)← Peack(Amp, fc, fb/rid);

20: A← A ∗ a;

21: B ← B ∗ b;22: nn← nn+ 1;

23: end if

24: end for

25: Hi(z)← B(z)A(z)

;

26: PRTFi ← PRTFi−1/Hi; —Aggiornamento—27: Hmulti−notch ← Hmulti−notch ∗Hi;

28: i← i+ 1;

29: end while


nn ← il numero dei notch rilevati algebricamente dallo spettro da

analizzare (PRTF0);

i← 0, numero iterazioni iniziali;

I parametri che contribuiscono a tarare il sistema di analisi sono:

Nceps: numero di coefficienti cepstrali considerati per la stima dell’invi-

luppo della PRTFi;

rid: fattore di riduzione della banda passante dei notch che compongono il

filtro multi-notch in ogni iterazione;

maxit: limite massimo delle iterazioni effettuate dall’algoritmo;

soglia: soglia minima per l’ampiezza dei notch compensativi da considerare.

Ogni parametro ha un effetto ben preciso sulla resa dell’algoritmo e verra

esposto nei prossimi paragrafi in relazione alla fase in cui e contenuto.

4.1.2 Calcolo del residuo

Il calcolo del residuo viene posto in relazione al Cepstrum ricavato dalla

f.d.t (funzione di trasferimento) presa in considerazione all’i-esima iterazione

(PRTFi). Il codice MatlabTM che implementa questo passaggio e il seguente:

flog = A HRTF;

N1 = 4;

cep = ifft(flog);

cep cut = [cep(1);2*cep(2:N1);cep(N1+1);

zeros(length(cep)−N1−1,1)];flog cut = real(fft(cep cut));

Dove flog rappresenta la log magnitude di PRTFi e N1 e il numero di

coefficienti cepstrali. Il procedimento per ricavare il cepstrum e illustrato

nel par. 3.1. A seconda della scelta del numero di coefficienti cepstrali per


CEPSi (individuato nel codice MatlabTM da flog-cut) si avra un inviluppo

spettrale piu dettagliato (parte reale del cepstrum).

Il numero di coefficienti cepstrali coincide approssimativamente con il

numero di picchi con cui viene approssimata la f.d.t (funzione di trasferimen-

to) presa in esame. Da questa osservazione possiamo scegliere il numero di

coefficienti in base alle risonanze che cerchiamo alla luce degli studi acustici

compiuti negli anni. Tenendo in considerazioni il lavoro di Shaw possiamo si-

curamente individuare una risonanza sempre presente (Modo 1) e due gruppi

di risonanze, quelle trasversali (Modi 2 e 3) e quelle orizzontali (Modi 4, 5

e 6) per angoli di elevazioni diverse (v. par. 1.2.4). Quindi per elevazioni

prossime allo 0, potremmo trovare il contributo di 3-4 risonanze (Modo 1, 4,

5 e 6) mentre per altri angoli solamente 2.

Questo ci puo guidare nel scelta di 4 coefficienti cepstrali. Il comporta-

mento dell’andamento del cepstrum nell’evoluzione dell’algoritmo verra espo-

sto in Appendice A. Il caso in cui la presenza di un numero di risonanze sia

minore di 4, congiunto all’utilizzo di 4 coefficienti cepstrali (caso di sovra-

stima del numero di risonanze), sara inserito in un’ottica piu generale, dopo

aver spiegato in maniera completa tutto l’algoritmo.

4.1.3 Ricerca parametri per il filtro multi-notch

Questa fase e molto delicata e determinante per la bonta dei risultati. Si

cercano i valori dei parametri che serviranno alla costruzione del filtro multi-

notch: frequenza centrale (fc), ampiezza (Amp) e banda passante (fb) per

ogni notch considerato nell’i-esima iterazione.

Le frequenze centrali dei notch vengono ricavate sui dati della PRTFi

tramite un algoritmo di notch-picking (v. par 3.3.1), con un numero massimo

di 20 notch individuabili per ogni iterazione. Le ampiezze dei notch vengono

ricavare valutando la risposta in ampiezza del residuo RESi nelle frequenze

centrali dei notch:

Amp(j) = |RES(fc(j)|, j = 1, ...N (4.1)


Figura 4.1: Caso 1) E(n) < 0, E(pl) > 0, E(pr) > 0, Caso 2) E(n) <

0, E(pl) < 0, E(pr) < 0, Caso 3) E(n) < 0, E(pl) < 0, E(pr) > 0, Caso 4)

E(n) < 0, E(pl) > 0, E(pr) < 0.

con N numero di notch individuati nella PRTFi. Per stimare la banda

passante dei noth si analizza l’energia del residuo:

ERESi = |RESi|2 (4.2)

Def. Per Banda Passante (fb) si intende 4f = f2− f1 dove f2 e f1 sono

frequenze ai 3db, cioe che possiedono, rispetto al picco/notch di frequenza

centrale fc, un’energia dimezzata (se si tratta di picco) o raddoppiata (se si

tratta di notch).

La ricerca delle frequenze ai 3db dipende dal target energetico (etarget) in

modo da ottenere:

ERESi(f2) ≈ ERESi(f1) ≈ etarget (4.3)

Il modo di scegliere etarget e fondamentale e viene calcolato a seconda di

come si presenta lo spettro del residuo, iterazione dopo iterazione; possiamo

riconoscere 4 casi principali:

1. notch sotto i 0 db di energia e picco destro e sinistro sopra i 0 db;

(E(n) < 0, E(pl) > 0, E(pr) > 0)


2. notch sotto i 0 db di energia e picco destro e sinistro sotto i 0 db;

(E(n) < 0, E(pl) < 0, E(pr) < 0)

3. notch sotto i 0 db di energia, picco destro sopra i 0 db e picco sinistro

sotto i 0 db; (E(n) < 0, E(pl) < 0, E(pr) > 0)

4. notch sotto i 0 db di energia, picco destro sotto i 0 db e picco sinistro

sopra i 0 db; (E(n) < 0, E(pl) > 0, E(pr) < 0)

In Fig. 4.1 riassumiamo graficamente i 4 casi dando un nome alle variabili

coinvolte per agevolare la spiegazione. Il seguente frammento di codice viene

applicato ai casi 1,3,4.

target = 0;

nEnergia = E res2(nftloc(no)); %−no db

nEnergia db = E res(nftloc(no));

if (amp ≤ −5) %−− Caso 1

target = nEnergia db + 3;

else

target = (1 − nEnergia)/2 + nEnergia; %−no db

target = 10*log10(target);

end

Nel caso il notch abbia ampiezza minore di -5 db si andranno a cercare le

frequenze f2 e f1 con energia doppia rispetto al notch di riferimento. Per

notch di ampiezza maggiore di -5 db si potrebbero ottenere delle energie

target in db positive non considerate valide, in quanto l’incremento deve

mantenersi sotto i 0 db1. Abbiamo imposto il limite di 0 db di energia

come banda-passante di riferimento calcolando l’energia target non come

1Ogni compensazione effettuata da un notch non deve portare la PRTFi al superamento

dell’inviluppo cepstrale nella banda che coinvolge. Ovviamente ci sono parti di spettro

della PRTFi che sono sopra l’inviluppo, ma al crescere delle iterazioni si assottigliano per

giungere alla convergenza.


raddoppio, ma come incremento relativo tra notch e soglia a 0 db.

etarget = 10log10

((1− Eno−db

RESi(fc))

2+ Eno−db

RESi(fc)

)(4.4)

I casi 3 e 4 vengono ricondotti al primo in quanto consideriamo, come pos-

sibile andamento del notch compensativo, la parte che attraversa gli 0 db

per giungere al picco positivo. La differenza rispetto al caso 1 sta nel calcolo

della larghezza di banda:

Caso 3⇒ fb = 2(f2 − fc) (4.5)

Caso 4⇒ fb = 2(fc − f1) (4.6)

Come si puo notare dallo spezzone di codice, nel caso 3 non si effettuera la

ricerca di f1, mentre nel caso 4 quella di f2, rispettivamente a sinistra e a

destra di fc.

e max peack db = max([e peak sx e peak dx]);

[e min peack db,in] = min([e peak sx e peak dx]);

no left = 1; no right = 1; %−− Caso 3−4if (e min peack db < 0) && (e max peack db > 0)

if (in == 1)

no left = 0;

elseif (in == 2)

no right = 0;

end

end

Il caso 2 e quello piu particolare, perche viene effettuato un cambio di ri-

ferimento: viene spostato in corrispondenza del picco piu vicino alla soglia

degli 0 db e piu vicino in frequenza a fc; viene poi individuato etarget come

nel caso 1, con unico accorgimento il cambio di riferimento, ed infine si effet-

tua la ricerca di f1 e f2. Giustifichiamo questo trattamento speciale perche

consideriamo il caso 2 come una situazione di passaggio. Vogliamo ottenere


Figura 4.2: Subject 010 left pinna φ = 51 - Contributo dei notch ottenuto

alla fine dell’algoritmo di analisi.

un filtro multi-notch con banda di riferimento a 0 db, approssimativamente

piatta (a 0 db) dove non ci sono interferenze da parte dei notch nella PRTF

(v. Fig. 4.2).

if (e max peack db < 0) %−− Caso 2

e max peack = 10ˆ(e max peack db/10);

amp max peack db = 20*log10(sqrt(e max peack));

amp rel = amp − amp max peack db;

if (abs(amp rel) ≥ soglia) %−− se false,0db come riferimento

amp = amp rel;

if (amp ≤ −5)target = amp + 3;

else

target = (e max peack − nEnergia)/2 + nEnergia;


end

end

end


4.1.4 Costruzione del filtro multi-notch

In quest’ultima fase i parametri stimati vengono elaborati per la costruzione

del filtro multi-notch compensativo, relativo all’iterazione corrente. Per pri-

ma cosa, viene controllata l’ampiezza di ogni notch. Se l’ampiezza relativa

di un notch e superiore al parametro soglia, viene creato un notch filter (sin-

golo) a guadagno controllato con i tre suoi parametri corrispondenti: Amp,

fc e fb. L’algoritmo utilizzato implementa un peak filter 2 (v. par. 3.2) del II

ordine:

(b, a)← PeackF ilter(Amp, fc, fb/rid) (4.7)

dove b e a sono i vettori dei coefficienti del filtro creato. Ogni notch (nk), che

soddisfa la condizione di soglia, possiede la coppia di vettori (bk,ak) che lo

rappresenta. Il filtro multi-notch compensativo (Hi) viene composto dal pro-

dotto delle f.d.t dei notch filter creati e i suoi coefficienti sono la convoluzione

dei b(i)k e a

(i)k :

Hi =B(i)

A(i)(4.8)

B(i) = b(i)1 ∗ ... ∗ b

(i)N (4.9)

A(i) = a(i)1 ∗ ... ∗ a

(i)N (4.10)

con N numero di notch considerati nell’i-esema iterazione.

Il parametro soglia permette di regolare la precisione dell’algoritmo nella

stima dell’inviluppo e il grado di convergenza. Piu la soglia e bassa piu

si possono considerare notch piccoli. Aumentare la soglia ha due effetti

principali:

1. per arrivare alla convergenza l’algoritmo impiega piu iterazioni, in quan-

to prende in considerazione un numero di notch piu elevato;

2. il filtro multi-notch finale risulta piu affinato;

2Un peak filter con guadagno negativo viene considerato come un filtro notch a

guadagno controllato.


Il parametro rid rappresenta il fattore di riduzione per la banda passante

dei notch, considerata nella costruzione dei notch filter. In pratica, questo

parametro permette di regolare il numero di notch filter impiegati per la

compensazione di un singolo notch, presente nella PRTF. L’introduzione di

questo parametro e stata resa necessaria analizzando la compensazione ef-

fettuata con le bande passanti calcolate alla fase precedente3. Si e osservato

che per notch molto profondi nella PRTFi, il notch compensativo, creato

dall’algoritmo peak filter, sovrastima la banda passante sopra i 3 db; si ottie-

ne quindi un’interferenza del notch compensativo sulle frequenze adiacenti,

aumentando il guadagno dove non e richiesto e creando un innalzamento di

CEPSi. Questo comportamento provoca un effetto a catena nelle iterazioni

successive che cercano di raggiungere la stima dell’inviluppo che si e alzata

in maniera innaturale. La relazione tra banda passante e peak filter verra

ripreso e approfondito nel capitolo 5.

Alzare il fattore di riduzione ha due effetti principali:

1. per compensare un notch della PRTF l’algoritmo impiega piu iterazioni,

in quanto vengono creati un numero di notch piu elevato;

2. il filtro multi-notch finale risulta piu affinato nell’individuazione dei

notch, ma con un andamento meno delineato;

L’andamento seghettato e evidente nella Fig. 4.3 per rid = 4 i i = 4.

I risultati al variare di questi parametri verranno esposti e commentati in

dettaglio in Appendice A.

4.1.5 Risultato finale

Esponiamo a titolo d’esempio la seguente situazione, rappresentata dai grafici

in Fig. 4.4 e 4.5, dove l’algoritmo e tarato con parametri sperimentalmente

buoni.

Subject 010, left pinna, φ = −45;

3Si potrebbe pensare di porre rid = 1, in modo da rendere neutrale il suo contributo.


(a) Situazione iniziale, i = 0.

(b) 4 iterazioni i = 1, ..., 4.

Figura 4.3: Subject 010 left pinna φ = −45 particolare 6000-8000 Hz, notch

di profondita -18 db - Prime 4 iterazioni di compensazione di PRTFi al

variare del parametro rid.


Nceps = 4;

soglia = 10−3;

rid = 2;

Convergenza dell’algoritmo dopo 59 iterazioni;

In Fig. 4.4 soffermiamo la nostra attenzione sul grafico del residuo e della

sua energia (i due grafici a destra). Possiamo notare come abbiano un anda-

mento approssimativamente piatto nella banda di frequenze 4 - 16 kHz, che

corrisponde alla banda di interesse dove agisce la pinna (v. par 1.2.4). Pos-

siamo quindi supporre che il risultato ottenuto dall’algoritmo di analisi sia

attendibile per quel range di frequenze, obiettivo della nostra analisi. Se si

desidera ampliare il range di frequenze, e necessario aumentare il numero di

coefficienti cepstrali per la rappresentazione dell’inviluppo, provocando pero

dei risultati contrastanti (come vedremo in Appendice A).

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Re(CEPSi)

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Inviluppo spettrale

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

|RESi| = |PRTF

i| − |CEPS

i|

|PRTFi|

Re(CEPSi)

Subject: 010 − Elevation: −45

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

Enr

gy v

alue

− d

b

Energia del residuo

Figura 4.4: Subject 010 left pinna φ = −45 - Risultati al termine

dell’algoritmo iterativo di analisi.

La Fig. 4.5 contiene i seguenti grafici (da sinistra a destra):


la PRTF iniziale dove sono contrassegnati i notch e le risonanze signi-

ficative;

il contributo spettrale delle risonanze ricavato dalle compensazioni;

lo stato iniziale e finale dell’inviluppo cepstrale;

il contributo spettrale dei notch ricavato dalla somma dei filtri multi-

notch generati per ogni iterazione;

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20Left HRTF − windowed signal

f (Hz)

dB

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Evoluzione Inviluppo spettrale

BEGINEND

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

dbContributo spettraledelle risonanze

Subject: 010 − Elevation: −45Subject: 010 − Elevation: −45

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Contributo spettrale dei notch

Figura 4.5: Subject 010 left pinna φ = −45 - Risultati al termine

dell’algoritmo iterativo di analisi.

Come prima cosa possiamo notare nuovamente come l’andamento dell’in-

viluppo condizioni la banda di lavoro dell’algoritmo. Un’osservazione mol-

to importante va fatta in relazione all’ultimo grafico: la composizione di

tanti filtri multi-notch di compensazione da come risultato un filtro multi-

notch4; abbiamo avuto una conferma incoraggiante dell’idea che ha guidato

la progettazione della modalita di compensazione della PRTF.

4Riproponiamo brevemente la definizione di filtro multi-notch: tanti notch rappresen-

tati in un unico filtro; l’interferenza tra le bande passanti dei notch e minima e lo spettro

si presenta approssimativamente piatto al di fuori di esse.


4.2 Analisi dei dati

In questo paragrafo verranno esposti e commentati i risultati ottenuti con

l’algoritmo d’analisi che abbiamo sviluppato.

Consideriamo come configurazione di riferimento l’algoritmo che possiede

i parametri con i seguenti valori:

Nceps = 4;

soglia = 10−3;

rid = 2;

Per un primo test sul CIPIC database abbiamo scelto 4 soggetti, 3 gia citati

nel paper di Raykar et el [25] (Subject 010, Subject 134 e Subject 165 ) e

uno utilizzato da Satarzadeh [22] (Subject 048 ), in modo da poter avere

dei punti di riferimento. Per ogni individuo sono state analizzate le HRTF

dell’orecchio sinistro per elevazioni tra i φ = −45 e φ = 90 e θ = 0. A titolo

di esempio, in Fig. 4.6, sono riportati i grafici tridimensionali che contengono

separatamente l’evoluzione spaziale del contributo delle riflessioni (a) e delle

risonanze (b). Ci focalizzeremo prevalentemente sulla banda d’interesse dei

−50

0

50

100

00.5

11.5

22.5 x 10

4

−40

−30

−20

−10

0

10

f (hz)

Notches LEFT−hrir − Subject: 010

elev (gradi)

mag

nitu

de

−40

−35

−30

−25

−20

−15

−10

−5

0

(a) Grafico 3D, Contributo dei notch

in relazione all’elevazione

−50

0

50

100

00.5

11.5

22.5x 10

4

−60

−40

−20

0

20

elev (gradi)

Risonanze LEFT−hrir − Subject: 010

f (kHz)

mag

nitu

de

−60

−50

−40

−30

−20

−10

0

10

20

30

(b) Grafico 3D, Contributo delle riso-

nanze in relazione all’elevazione

Figura 4.6: Subject 010 left pinna - Grafici 3D.

3 - 18 kHz. Cercheremo di interpretare i risultati in modo da dimostrare la

valenza del procedimento che e stato assemblato.

4.2. ANALISI DEI DATI 89

4.2.1 Analisi dei notch

Per l’analisi dei risultati ci serviamo innanzitutto di un grafico tridimensio-

nale (Fig. 4.6.a) in cui sull’asse delle ascisse c’e la frequenza, sull’asse delle

ordinate l’elevazione e sull’asse z (secondo la colorbar) l’ampiezze in db. Il

colore rosso scuro corrisponde ad ampiezze prossime ai 0 db mentre man

mano che il colore diventa freddo aumenta la profondita del notch.

Essendo il nostro studio interessato ai segnali di localizzazione verticale

abbiamo cercato un modo per monitorare l’evoluzione spaziale dei notch.

A tale scopo abbiamo utilizzato la tecnica di partial tracking, che solitamente

viene utilizzata in relazione al tempo (v. par. 3.3.2), nel dominio dello spazio

verticale. L’idea che ha guidato la scelta di questo strumento e stata la

possibilita di collegare i notch tra loro tracciandoli lungo traiettorie comuni

(Fig. 4.6.b).

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0

(a) Contributo notch

−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(b) Notch tracking (c) Contributo notch + Not-

ch tracking

Figura 4.7: Subject 010 left pinna - Contributo dei notch in relazione

all’elevazione.

Iniziamo considerando il Subject 010. Per questo soggetto possiamo fare

almeno 3 osservazioni:

il notch che si sviluppa tra i 1 - 1.2 kHz per elevazioni dai −45 ai

0 e piuttosto profondo (-15 db ca) e, secondo l’algoritmo d’analisi,

la sua banda passante e di 1 kHz. E’ molto facile accorgersi, dalla

sovrapposizione dei due grafici (Fig. 4.6.c), che per elevazioni dai −45

ai−20 la banda passante e molto piu ampia. Raykar infatti individua 2


notch ravvicinati in quell’intervallo di frequenze. Nel capitolo dedicato

alla modellizzazione vedremo come si comporta l’approssimazione per

questo caso particolare.

A elevazione 45 c’e un notch a profondita -46 db che sembra apparire

come somma di due tracce di notch alle elevazioni precedenti. Purtrop-

po lo spazio di ripresa dei dati del database e campionato con passo

troppo grande per comprendere a fondo questa situazione.

Tra elevazioni 35 e 90 c’e una traccia di notch dell’ampiezza di -1 db.

Potrebbe essere un effetto molto attenuato dell’elice vista la relativa

bassa frequenza; la sua importanza percettiva e da provare.

A fine capitolo sono esposti i risultati per gli altri tre soggetti (il soggetto

165 e un KEMAR). Le considerazioni sono analoghe e si puo vedere una

corrispondenza incoraggiante con i risultati di Raykar. Il tracking ci puo

dare le informazioni per collegare i tracciati rilevati con la forma della pinna.

4.2.2 Analisi delle risonanze

Per l’analisi delle risonanze ci siamo serviti, come per i notch, di un grafico

tridimensionale (Fig. 4.7.a) in cui sull’asse delle ascisse c’e la frequenza, sul-

l’asse delle ordinate l’elevazione e sull’asse z (secondo la colorbar) l’ampiezze

in db. Gli estremi della scale di colore sono [-60 db, 30 db], piu il colore e

caldo piu la risposta in ampiezza della risonanza e elevata. Secondo la scale

adottata, il colore giallo rappresenta gli 0 db.

Un oggetto risonatore (nel nostro caso l’orecchio) possiede struttural-

mente delle risonanze. Ogni risonanza ha delle frequenze caratteristiche

che vengono stimolate a seconda dell’elevazione. Nella Fig. 4.7.a possiamo

distinguere, mediante un’ispezione visuale, almeno 2 risonanze principali:

1. una risonanza onnipresente a ∼ 4kHz di intensita tra i 10 e 14 db;

2. una risonanza che va scomparendo per φ grandi (φ > 50 con frequenza

centrale ∼ 1.7kHz;


−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (kHz)


−60

−50

−40

−30

−20

−10

0

10

20

30

(a) contributo risonanze

−40 −20 0 20 40 60 80 1000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2x 10

4

Elevation

Fre

quen

cy (

Hz)

Peak tracking

0

2

4

6

8

10

12

14

16

18

20

(b) peak tracking

Figura 4.8: Subject 010 left pinna - Contributo delle risonanze in relazione

all’elevazione.

Per calcolare le frequenze centrali abbiamo usato un sistema di identi-

ficazione che utilizza un modello ARMA del VI ordine, proposto per una

eventuale modellazione del contributo delle risonanze e che verra trattato nel

cap. 5. E’ stato effettuato il tracking nello spazio (Fig. 4.7.b) per le frequen-

ze ricavate dai poli complessi coniugati del modello e associato ad ognuna di

esse un pallino colorato secondo il valore del modulo.

Abbiamo messo in relazione le due risonanze ottenute con il risultato teo-

rico ottenuto da Lopez-Poveda e Meddis [26] e osserviamo le corrispondenze

con la semplificazione operata da Satarzadeh [22]. Entrambi considerano 2

risonanze principali5 una orizzontale, associata alla profondita della conca,

e una verticale associata alla larghezza della conca. Notiamo comunque un

risultato coerente nell’individuazione dei modi di Shaw. Per elevazioni pros-

sime ai 70, la banda della prima risonanza e molto ampia, facendo pensare

alla presenza di altre risonanze (Modo 2 e 3) a frequenze vicine che contri-

buiscono alla sua espansione. Un discorso analogo puo essere fatto per la

seconda risonanza ad elevazioni prossime allo 0 (Modi 4,5 e 6).

A fine capitolo sono esposti i risultati ottenuti per gli altri tre soggetti. Vo-

5Lopez-Poveda e Meddis considerano principalmente i notch di un modello di conca

privo delle altre componenti di orecchio esterno. Quindi si puo relazionare ogni notch ad

una risonanza corrispondente.


gliamo far notare come per il soggetto 165 (testa KEMAR) siano accentuati

i tratti che riguardano la sovrapposizione degli effetti di risonanze vicine.


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 80 1000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8x 10

4

Elevation

Fre

quen

cy (

Hz)

Peak tracking

0

2

4

6

8

10

12

14

16

18

20

(b) peak tracking


all’elevazione.

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 80 1000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8x 10

4

Elevation

Fre

quen

cy (

Hz)

Peak tracking

0

2

4

6

8

10

12

14

16

18

20

(b) peak tracking


all’elevazione.

−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (kHz)


−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 80 1000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8x 10

4

Elevation

Fre

quen

cy (

Hz)

Peak tracking

0

2

4

6

8

10

12

14

16

18

20

(b) peak tracking


all’elevazione.


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking


ch tracking


all’elevazione.

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking


ch tracking


all’elevazione.

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking


ch tracking


all’elevazione.

Capitolo 5

Un Modello Strutturale per la

pinna

In questo capitolo conclusivo introduciamo un primo modello per sintetizzare

la PRTF. Ci poniamo in un’ottica di ricerca futura, in quanto alcune delle

considerazioni che riportiamo sono delle ipotesi da verificare rigorosamente.

Utilizzeremo i risultati ottenuti dalla fase di analisi dei dati, grazie alla

quale possiamo trattare distintamente il contributo dovuto alle risonanze e

quello dovuto ai notch. Facciamo rientrare la nostra proposta all’interno del-

la tipologia modelli strutturali per due motivi: (1) colleghiamo assieme dei

fenomeni acustici trattati in maniera indipendente; (2) la progettazione di un

modello strutturale permette di individuare un insieme di parametri control-

labili da grandezze antropometriche direttamente misurabili sull’ascoltatore,

caratteristica che garantisce la personalizzazione della PRTF sintetizzata.

5.1 Descrizione e motivazione

Le caratteristiche principali di un modello strutturale (una PRTF sintetica)

che ne guidano la progettazione sono le seguenti:

possibilita di individuare un mapping semplice tra parametri del mo-

dello e grandezze antropometriche;

95

96 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA PINNA

grandezze antropometriche facilmente reperibili;

complessita ridotta per applicazioni real-time;

efficacia dal punto di vista acustico-percettivo;

Un primo tentativo di modello proposto, deriva dall’approssimazione struttu-

rale di una PRTF misurata e viene ricavato dalla catena di operazioni descrit-

ta nello schema a blocchi di Fig. 5.1. L’ingresso HRIR viene finestrato, per

Figura 5.1: Schema a blocchi del sistema di costruzione di una PRTF

sintetica.

poi essere ricavata la PRTF tramite la Trasformata di Fourier. L’algoritmo

di analisi viene applicato alla PRTF e restituisce come risultati i contribu-

ti dei notch e delle risonanze e il loro andamento spaziale, disassociando le

interferenze tra gli effetti di riflessione, risonanza e diffrazione.

Viene spontaneo cercare di modellare le risonanze e i notch con due fil-

tri indipendenti e totalmente controllabili nei loro parametri. Nel prossi-

mo paragrafo vediamo approfonditamente una loro descrizione e modalita di

collegamento.

5.2 Filtro sintetico

Per filtro sintetico intendiamo una combinazione di componenti manipolabili

a piacimento che possano approssimare in maniera convincente una PRTF.

5.2. FILTRO SINTETICO 97

La struttura che proponiamo e rappresentata in Fig 5.2. Il filtro sintetico

Figura 5.2: Filtro sintetico.

e composto da un filtro ARMA del VI ordine collegato in serie ad un filtro

multi-notch, il tutto viene posto in cascata con un filtro passa banda (BP, v.

par. 5.2.3). Cerchiamo di comprenderne a fondo il funzionamento.

5.2.1 Filtro sintetico per le risonanze

Il filtro sintetico per la rappresentazione delle risonanze e un del tipo AR-

MA(p,q) con p=6 (poli) e q=6 (zeri):

Hrisonanze(z) =Bq(z)

Ap(z)=

∑qk=1 bq(k)z−k∑pk=1 ap(k)z−k

(5.1)

I coefficienti bq e ap vengono calcolati da un metodo di identificazione che

viene applicato alla risposta all’impulso, ottenuta antitrasformando la f.d.t

contenente il contributo delle risonanze. Le tecniche di identificazione che

sono state applicate sono STMB (Steiglitz-McBride iteration method) [9][34]

e Linear Prediction.

Illustriamo due situazioni particolari in cui le risonanze da approssimare

sono ravvicinate e analizziamo il comportamento del modello ARMA. In

Fig. 5.3, nell’intervallo di frequenze 12 - 16 kHz, possiamo riconoscere i

contributi dei Modi 4,5 e 6 di Shaw (elevazione prossima agli ∼ 0); in Fig.

5.4, nell’intervallo di frequenze 4 - 10 kHz, possiamo riconoscere i contributi

dei Modi 1,2 e 3 (elevazione prossima ai ∼ 70). Possiamo definire queste

due situazioni come due zone critiche1 per il corpo risonatore orecchio, in

1Si parla di zone critiche in quanto la risonanza ha un’intensita massima per una certe

elevazione, al di fuori della quale il contributo e minore ma non nullo.


0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − Contributo risonanza

f (Hz)

dB

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − Risonanze − AR model

−1 −0.5 0 0.5 1−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Real Part

Imag

inar

y P

art

Left HRTF − Risonanze − Z−plane

Figura 5.3: Subject 134 left pinna φ = 6, Filtro risonanze.

quanto piu risonanze contribuiscono allo spettro finale, interagendo tra loro.

Il filtro risonanze (filtro sintetico per le risonanze) approssima media-

mente bene l’andamento ottenuto dall’algoritmo di analisi e possiamo rico-

noscere dal grafico poli-zeri un accoppiamento tra coppie di poli e coppie di

zeri. Il modello e cosı scomponibile in tre picchi ognuno dei quali della forma:

Hpeak(z) =B2(z)

A2(z)=

∑2k=1 b2(k)z−k∑2k=1 a2(k)z−k

(5.2)

Dal tracking spaziale dei poli, che richiamiamo dal capitolo 4, possiamo

studiare l’andamento dei picchi di risonanza.

l’intensita delle risonanze mediamente diminuisce all’aumentare dell’e-

levazioni;

la risonanza di Modo 1 e sempre presente e modellata da un picco di

frequenza ∼4.2 kHz;

i restanti modi vengono gestiti dagli altri due picchi; in particolare

un picco di frequenza ∼17 kHz che aggiusta, con guadagno inversa-


0.5 1 1.5 2

x 104

−30

−20

−10

0

10


f (Hz)

dB

0.5 1 1.5 2

x 104

−30

−20

−10

0

10


−1 −0.5 0 0.5 1−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Real Part

Imag

inar

y P

art


Figura 5.4: Subject 165 left pinna φ = 68, Filtro risonanze.

mente proporzionale all’elevazione, l’andamento del modello nelle alte

frequenze, e un picco variabile tra 8 e 14 kHz nelle medie frequenze.

Mentre il primo e il terzo picco hanno un comportamento stabile nelle PRTF

modellate, il secondo ha un andamento che varia da persona e persona.

Il secondo picco, per elevazioni basse (φ = −45), ha una frequenza di

risonanza mediamente di ∼ 12 kHz, mentre, per elevazioni alte (φ = 90)

prossime ai ∼ 9 kHz. Il passaggio tra queste due frequenze, all’aumentare

dell’elevazioni, varia molto tra i soggetti:

Subject 010 tra φ = 0 e φ = 34;




Per ora, non sappiamo dire esattamente i fattori che determinano l’andamen-

to del secondo picco, ma per avere un controllo totale sui parametri sara uno

dei primi studi da effettuare.


−40 −20 0 20 40 60 80 1000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8x 10

4

Elevation

Fre

quen

cy (

Hz)

Peak tracking

0

2

4

6

8

10

12

14

16

18

20

(a) Subject 134 left pinna.

−40 −20 0 20 40 60 80 1000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8x 10

4

Elevation

Fre

quen

cy (

Hz)

Peak tracking

0

2

4

6

8

10

12

14

16

18

20

(b) Subject 165 left pinna.

Figura 5.5: Peak Tracking nello spazio.

Abbiamo scelto questo tipo di filtro per due motivi principali: 1. l’e-

ventuale possibilita di modellare le risonanze come combinazioni di peak

filter, approccio speculare a quello elaborato per i notch; 2. un modello

di ordine contenuto che approssima in maniera soddisfacente i dati ricavati

dall’algoritmo di analisi.

5.2.2 Filtro sintetico per i notch

Il filtro che modella i notch e ricavato dalla re-sintesi della f.d.t ottenuta dal

prodotto dei filtri multi-notch compensativi ottenuti durante le iterazioni

dell’algoritmo di analisi. Il procedimento di sintesi e molto semplice: le

regole sone le medesime di quelle utilizzate per la costruzione dei filtri multi-

notch di compensazione. Vi sono solamente alcune osservazioni da fare che

distinguono le due fasi:

1. per il calcolo della banda passante dei notch viene considerato solamen-

te il caso 2; non vi sono picchi al di sopra degli 0 db, quindi la ricerca

viene condotta con banda passante di riferimento a 0 db o inferiore.

2. non viene impiegato il fattore di riduzione; nella costruzione del notch

filter, viene impiegata esattamente la banda passante calcolata;

3. non vi e alcuna soglia in ampiezza;


Algorithm 5.1 Costruzione del filtro sintetico per notch, Fasi

1 - Calcolo dei notch piu importanti mediante l’algoritmo notch pick ;

for j = 1 to length(notch) do

if f jc ≤ 18kHz then

(aj, bj)←notch filter con i parametri ricavati come Alg. 4.2, caso 2;

A← A ∗ aj;B ← B ∗ bj;

end if

end for

3 - Costruzione filtro sintetico multi-notch (B,A);

4. viene introdotta una soglia in frequenza (fc ≥18 kHz), per elimina-

re l’influenza dei notch ad alta frequenza, dove sono percettivamente

indistinguibili.

Il filtro che otteniamo e di ordine 2 ∗N , con N numero di notch considerati

validi.

In Fig. 5.6 e 5.7 sono illustrati i risultati per due soggetti, 048 e 165,

a due elevazioni φ = −34 e φ = 79 rispettivamente. Il grafico in alto a

sinistra rappresenta l’output dell’algoritmo di analisi, mentre quello in basso

a sinistra la ri-sintesi. I due casi sono stati scelti per confrontare due livelli

di intensita riflessiva: nel primo caso (φ = −34) sono presenti molti notch

di profondita elevata, nel secondo (φ = 79) pochi notch e limitati.

Analizziamo il comportamento dell’operazione di ricostruzione nelle basse

e nelle alte frequenze. Come sottolineato in Alg. 5.1, la soglia di 18 kHz fa

in modo che non vengono considerati notch come quello a frequenza 19.5

kHz (Fig. 5.6). A 1.6 kHz e presente un picco di +2.3 db che non viene

considerato in quanto sopra la soglia degli 0db. La causa della sua presenza e

l’inaccuratezza dell’algoritmo di analisi fuori dalla banda di lavoro 4 - 16 kHz;

parte del notch a 2.2 kHz, presente nella PRTF iniziale (Fig. 5.11.b), viene

compensato con un notch a 2.3 kHz, mentre, erroneamente, la parte restante

con un picco a 1.6 kHz. Considerando la parte compensata adeguatamente


0.5 1 1.5 2

x 104

−40

−30

−20

−10

0

Left HRTF − Contributo notch

f (Hz)

dB

0.5 1 1.5 2

x 104

−40

−30

−20

−10

0

Left HRTF − Filtro Notch sintetico

−1 −0.5 0 0.5 1−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Real Part

Imag

inar

y P

art


Figura 5.6: Subject 048 left pinna φ = −34, Filtro notch.

si ha comunque una buona stima, anche se la valenza acustico-percettiva e

da provare.

La ricostruzione per notch poco profondi (Fig. 5.7) avviene con risultati

visivamente soddisfacenti. Gia da questa semplice situazione, possiamo pero

notare una caratteristica che diventa molto evidente in Fig. 5.6: la banda

passante per ogni notch, in prossimita degli 0db, e molto piu ampia rispetto

al dovuto. L’effetto e molto piu rilevante, man mano che si considerano notch

con profondita elevate. Il notch a frequenza 14.9 kHz di intensita -33.5 dB

viene ri-sintetizzato influenzando l’andamento di tutto lo spettro del filtro

globale.

Le forme dei notch incontrati in fase di analisi sono tra le piu svariate, e

abbiamo cercato di approssimarle con un notch filter a guadagno controllato

del secondo ordine con specifica di banda passante ai 3db. La molteplicita

di forme, fa in modo che la specifica di banda passante ai 3db sia troppo

restrittiva e non adatta a quelli con profondita elevate. Proponiamo due

possibili soluzioni:


0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − Contributo notch

f (Hz)

dB

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − Filtro Notch sintetico

−1 −0.5 0 0.5 1−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Real Part

Imag

inar

y P

art


Figura 5.7: Subject 165 left pinna φ = 79, Filtro notch.

1. per ogni notch calcolare una banda passante ad hoc, in modo da ap-

prossimarne al meglio l’andamento;

2. utilizzare un peak filter di ordine piu elevato per avere piu punti di

controllo nell’approssimazione;

5.2.3 Ricostruzione

In questo paragrafo assembliamo i due filtri sintetici in modo da ottenere

un’approssimazione delle PRTF. Descriviamo l’intero procedimento attraver-

so i commenti di un esempio completo riguardante Subject 134 per φ = 0.

In Fig. 5.8 sono riportati i risultati dell’algoritmo di analisi. Notiamo come

l’evoluzione dell’inviluppo, durante le iterazioni dell’algoritmo, sia evidente

sopra i ∼ 12 kHz e la presenza di 5 notch: tre molto pronunciati a 6.9 kHz,

10.9 kHz e 15.7 kHz e 2 meno profondi a 9.1 kHz e 19.3 kHz.

In Fig. 5.9 sono presentati i filtri sintetici per ognuno dei due contributi.

Il filtro per le risonanze (a) posiziona i 3 picchi a 3.7 kHz, 12.5 kHz e 17 kHz


0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − windowed signal

f (Hz)

dB

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Evoluzione Inviluppo spettrale

BEGINEND

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Contributo spettraledelle risonanze

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20

f (hz)

db

Contributo spettrale dei notch

Figura 5.8: Subject 134 left pinna φ = 0, Risultato ottenuti al termine

dell’algoritmo di analisi.

0.5 1 1.5 2

x 104

−30

−20

−10

0

10


f (Hz)

dB

0.5 1 1.5 2

x 104

−30

−20

−10

0

10


−1 −0.5 0 0.5 1−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Real Part

Imag

inar

y P

art


(a) Contributo Risonanze.

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − Contributo notch

f (Hz)

dB

0.5 1 1.5 2

x 104

−30

−20

−10

0

10

20Left HRTF − Filtro Notch sintetico

−1 −0.5 0 0.5 1−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Real Part

Imag

inar

y P

art


(b) Contributo Notch.

Figura 5.9: Subject 134 left pinna φ = 0, filtri sintetici


bilanciando il suo contributo nella parte medio-alta dello spettro2. Il filtro

per i notch (b) non considera il notch a frequenza 19.3 kHz e soffre della

sovrastima del notch a 15.7 kHz di profondita elevata, -22 db.

I due filtri vengono collegati in serie dando origine alla PRTF sintetica

riportata in Fig. 5.10. Ci limitiamo ad alcune osservazioni che restano da

provare attraverso studi futuri piu approfonditi.

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−25

−20

−15

−10

−5

0

5

10

15

20Left HRTF − Filtro Notch sintetico Subject: 134 − Elevation: 0

PRTF originalePRTF sintetica

Figura 5.10: Subject 134 left pinna φ = 0, Filtro sintetico per PRTF.

I notch vengono ben approssimati in frequenza e in ampiezza, mentre

l’andamento dei picchi viene contaminato dall’interferenza delle bande pas-

santi troppo ampie. L’andamento degli estremi della PRTF sintetica e stato

trascurato perche esterno alla banda di interesse. Post-porre un filtro passa

banda puo permettere l’inseguimento degli andamenti nel fondo scala; questa

operazione e desiderabile se il risultato globale deve venire assemblato con i

contributi di testa, busto e spalle.

A titolo di esempio, presentiamo il risultato finale di filtro sintetico per

tre situazioni gia affrontate nei paragrafi precedenti.

In Fig. 5.11.a altro esempio, in cui e presente la necessita di attenuare

le alte frequenze per un’eventuale inserimento all’interno di un modello

2Come ci aspettavamo dall’influenza dei Modi 4,5 e 6.


completo3.

In Fig. 5.11.b notiamo come l’interferenze di piu notch ad elevata

profondita generi un’approssimazione poco efficace.

In Fig. 5.11.c e riportato un esempio molto promettente: l’andamento

delle risonanze e ben rispettato e i notch, di media profondita, vengono

approssimati in maniera adeguata.

3Ricordiamo che l’energia audio sopra i 18 kHz e poca e percettivamente ha poca

importanza.


0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−25

−20

−15

−10

−5

0

5

10

15



(a)

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−25

−20

−15

−10

−5

0

5

10

15

20Left HRTF − Filtro Notch sintetico Subject: 048 − Elevation: −34


(b)

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 104

−30

−25

−20

−15

−10

−5

0

5

10

15



(c)

Figura 5.11: Alcuni esempi di filtri sintetici ottenuti nella fase di

modellazione, messi a confronto con le corrispettive PRTF reali.

Conclusioni e Sviluppi futuri

Il lavoro di tesi si e concentrato sulla ricerca di un nuovo approccio che

potesse contribuire al problema aperto riguardante la localizzazione verticale

del suono nella sintesi audio binaurale. Molte energie sono state impiegate

nello studio delle strade gia percorse in letteratura rivolte alla comprensione

dei fenomeni che coinvolgono la pinna.

E’ stato presentato un algoritmo di analisi il cui scopo e quello di stimare i

contributi distinti di risonanze e riflessioni dovuti all’orecchio esterno. I risul-

tati sono stati relazionati a conoscenze consolidate in letteratura; l’algoritmo

e stato progettato in maniera mirata in modo da ottenere degli input disas-

sociati da fornire ad un modello strutturale, nel quale ogni effetto acustico

viene considerato indipendentemente e successivamente collegato al sistema

globale.

Il modello strutturale proposto fornisce visivamente una buona appros-

simazione, anche se riconosciamo molti punti che necessitano di ulteriori

indagini. Per alcuni sono state fornite delle idee risolutive, per altri e nata

la necessita di uno studio approfondito. La problematica aperta piu im-

portante e complementare a questo lavoro riguarda la descrizione dell’an-

damento spaziale del filtro sintetico; variando questo considerevolmente da

persona a persona, risulta ancora difficile poterlo controllare con variazioni

parametriche.

Possiamo definire i risultati di questa tesi come promettenti. La loro diffi-

cile interpretazione non ci permette di stabilirne la loro reale bonta, tantopiu

che tuttora permane l’assenza di un criterio d’errore oggettivo che possa so-

109

110 CONCLUSIONI E SVILUPPI FUTURI

stituirsi alla validazione fornita dai test di ascolto. Il giudizio potra essere

rigoroso nel momento in cui si giungera ad una simulazione coinvolgente e

convincente per un ascoltatore o verra definito un modello di ascoltatore

adatto a tale scopo.

Sviluppi futuri

Gli studi da compiere per continuare questo lavoro sono molteplici: dal mi-

glioramento dell’algoritmo d’analisi al test del modello strutturale. Le strade

da percorrere sono numerose:

avendo a disposizione degli ascoltatori e le loro relative HRTF si do-

vranno condurre dei listening test per comprendere la bonta del modello

proposto;

avendo a disposizione solo i dati antropometrici e le HRTF corrispon-

denti, gli studi dovranno proseguire cercando la relazione tra i parame-

tri del modello in relazione alle grandezze della pinna;

sara necessario comprendere rigorosamente se e possibile disassociare

completamente i contributi delle risonanze da quelle dei notch e cercare

cosı nuove tecniche di compensazione piu idonee;

Scendendo piu nel particolare di questa tesi, le problematiche da risolvere

riguardano:

lo studio dell’ordine del modello per le risonanze; testare un modello

di ordine superiore per osservare l’andamento dei picchi e cercare una

spiegazione piu rigorosa per il loro andamento nello spazio verticale.

l’allargamento della banda di lavoro dell’algoritmo, in modo da avere

un’analisi tutto spettro;

la modellazione fedele del singolo notch, riducendo le interferenze in-

trodotte da un’approssimazione inappropriata;

CONCLUSIONI E SVILUPPI FUTURI 111

Per alcune di queste abbiamo proposto qualche idea nei capitoli preceden-

ti, documentando a fondo il lavoro svolto. Ma numerose problematiche ri-

guardanti il campo dell’audio binaurale restano aperte e affascinanti; sono

necessari ancora molti studi per scoprire la bellezza di essere ascoltatori.

112 CONCLUSIONI E SVILUPPI FUTURI

Appendice A

Algoritmo d’analisi, variazioni

parametriche

Approfondiamo gli effetti provocati dalle variazioni parametriche sull’algo-

ritmo d’analisi. Il punto di riferimento per queste sperimentazioni e la

configurazione usata nel capitolo 4. La richiamiamo per comodita:

Nceps = 4;

soglia = 10−3;

rid = 2;

A.1 Parametro Nceps

Le figure A.2 - A.5 riassumono i risultati per Nceps = 2, 3, 5, 7. L’effetto del-

l’incremento dei coefficienti cepstrali per la rappresentazione dell’inviluppo

della PRTFi e notevole. Questo parametro gioca il ruolo di vero e proprio

ago della bilancia nella determinazione dei rispettivi contributi (notch e ri-

sonanze), in modo da equilibrarne la stima. All’aumentare dei coefficienti

cepstrali:

il contributo dei notch diminuisce, sia in ampiezza, sia in banda pas-

sante;

113

114 APPENDICE

aumenta il dettaglio della mappa delle risonanze;

rimangono le tracce dei notch piu significativi con ampiezza ridotta ed

emergono nuove tracce di debole intensita;

Il bilanciamento tra contributo delle risonanze e contributo dei notch deve

essere guidato da considerazioni acustico-percettive. Secondo le considera-

zioni esposte nel par. 4.1.2, il valore piu indicato per il numero di coefficienti

cepstrali e Nceps = 4

A.2 Parametro soglia

Le figure A.6 - A.9 riassumono i risultati per soglia = 100, 10−1, 10−2, 10−6.

L’effetto di questo parametro e visibile solo nel passaggio da soglia = 1 a

soglia = 10−1. L’andamento delle risonanze possiede piu energia, mentre il

numero di notch individuati aumenta (compare la traccia per φ = 50 − 85

a frequenza ∼6 kHz, e la traccia per φ = 63 − 90 a frequenza ∼1.3 kHz).

Per soglia→ 0 l’effetto e sempre meno visibile, quindi possiamo scegliere

con tranquillita un valore di soglia < 10−2 che ci permetta di ottenere un

numero di iterazioni contenuto.

A.3 Parametro rid

Le figure A.10 - A.13 riassumono i risultati per rid = 1, 2, 4, 8. Ricordiamo

che questo parametro e stato introdotto per limitare l’influenza delle bande

passanti tra i notch di compensazione.

Per rid = 1 si nota in Fig. A.10 un andamento molto confuso dello spettro

riguardante i notch, con aree di giallo piuttosto estese (= bande passanti dei

notch estese). Il contributo delle risonanze presenta qualche anomalia.

Caso rid = 1, Esempio φ = 50:

sembrerebbe essere presente un’unica risonanza al centro dello spettro,

smentita dall’andamento delle risonanze per elevazioni inferiori e supe-

A.3. PARAMETRO RID 115

riori, nonche dagli altri grafici per valori di rid maggiori. Questo e un

esempio di effetto a catena provocato dalla compensazione esagerata

dei notch filter;

le ampiezze sono mediamente piu elevate (colore dei grafici piu caldo),

in quanto la sovra-compensazione dei notch provoca un innalzamento

innaturale dell’inviluppo cepstrale (v. par. 4.1.4);

Figura A.1: Subject 048 left pinna φ = −28 - rid = 8,Particolare 1 - 4 kHz.

Aumentando il valore di rid, la compensazione viene frammentata in piu

iterazioni, e il risultato finale dell’algoritmo puo produrre un andamento

segmentato nel filtro multi-notch, come mostrato nel particolare di Fig. A.1.

116 APPENDICE

−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (hz)


−40

−35

−30

−25

−20

−15

−10

−5

0

(a) contributo notch

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30

(b) contributo risonanze

−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking

Figura A.2: Subject 134 left pinna - Nceps = 2, Risultati.

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (hz)


−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


A.3. PARAMETRO RID 117

−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (hz)


−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking

Figura A.6: Subject 165 left pinna - soglia = 1, Risultati.

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking

Figura A.7: Subject 165 left pinna - soglia = 10−1, Risultati.

−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (hz)


−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (kHz)


−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


118 APPENDICE

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking

Figura A.10: Subject 048 left pinna - rid = 1, Risultati.

−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (hz)


−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5x 10

4

elev (gradi)

f (kHz)


−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (hz)

−40

−35

−30

−25

−20

−15

−10

−5

0


−50

0

50

100

00.511.522.5

x 104

elev (gradi)


f (kHz)

−60

−50

−40

−30

−20

−10

0

10

20

30


−40 −20 0 20 40 60 800

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x 104

Elevation

Fre

quen

cy (

Hz)

Notch tracking

(c) notch tracking


Appendice B

Sorgenti MATLAB

I sorgenti sviluppati per questa tesi utilizzano alcune function MatlabTM con-

tenute nel Spectral Modeling Synthesis Tools (SMS) (http://mtg.upf.edu/technologies/sms/,

The Music Technology Group (MTG), Universitat Pompeu Fabra di Barce-

lona) e negli strumenti di navigazione forniti dal CIPIC database. Per com-

pletezza, riportiamo tutto il codice necessario per poter comprendere a fondo

il lavoro svolto, lasciando i commenti originali nelle righe di terze parti.

% % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% Universita degli Studi di Padova %

% DEI − Dipartimento di Ingegneria dell'Informazione %

% % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% Author: Michele Geronazzo

% Date: 20−01−09

% Script per l'Analisi e Modellazione delle HRTF dei soggetti

% del CIPIC database

function analisi hrtf()

addpath('../');

addpath('./fz−arma');

119

http://mtg.upf.edu/technologies/sms/

120 APPENDICE

clear all;

close all;

%% −−−−−−−−−−−−−−−−−−−−−−−Inizializzazione−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

start flag = 1; % Prevents any action until data is available

log scale = 0; % Default to use a logarithmic frequency scale

do smooth = 0; % Default to use constant−Q smoothing

Q = Inf; % no smoothing Q = inf

fs = 44100; % Sampling frequency in Hz

fNy = fs/2000; % Nyquist frequency in kHz

Ia = 13; % Initial azimuth index (0 deg)

Ie = 9; % Initial elevation index (0 deg)

A = [−80 −65 −55 −45:5:45 55 65 80]; LA = length(A); % Azimuths

E = −45:(360/64):235; LE = length(E); % Elevation

LT = 200; % Number of time points

T = (0:(1/fs):((LT−1)/fs))*1000; % Times in ms

subject numbers = sub nums; % Load cell array of subject numbers from file

individuo = 1; % Count individui tracciati

NFFT = 2*LT; % More freq resolution if NFFT > LT,

% Slows computing High resolution

Fmin = 0.5; Fmax = 21.93975; % Min and max frequency in kHz

Tmin = 0.3; Tmax = 2.5; % Min and max time in ms

miT = max(find(T≤Tmin)); % Min time index

mxT = min(find(T≥Tmax)); % Max time index

dBmin = −30; % Lower limit for spectral plots in dB

dBmax = +20; % Upper limit for spectral plots in dB

% Caricamento dati delle HRIR per azimut = 0

% Find path to hrir data files

script name = 'analisi hrtf.m';

data path head = which(script name);

if isempty(data path head),

error(['Cannot find ' script name '. This should not happen!']);

APPENDICE 121

end;

num chars = length(data path head);

sep = data path head(num chars−length(script name));

sep locs = findstr(sep,data path head);

if length(sep locs) < 2,

error([script name ...

' ...cannot be at the root of the file directory tree.']);

end;

data path head = [data path head(1:sep locs(end−2)) 'subject '];

data path tail = [sep 'hrir final.mat'];

%% Caricamento e analisi di ciascuna hrir

for sub = 2 : length(subject numbers)

%subject string = subject numberssub;subject string = '134';

Fname = [data path head subject string data path tail];

load(Fname); % carica hrir l, hrir r, OnR, OnL

% Per il notch tracking in space

nNotches = 20;

ntrack n = 200;

Delta n = 1000;

previousnloc = zeros(ntrack n,1);

previousnval = zeros(ntrack n,1);

previousnband = zeros(ntrack n,1);

nloc = zeros(ntrack n,1);

nval = zeros(ntrack n,1);

nband = zeros(ntrack n,1);

emergingnloc = zeros(nNotches,1);

emergingnval = zeros(nNotches,1);

emergingnband = zeros(nNotches,1);

% Per il peack tracking in space

ncisoids = 3;

ntrack = 50;

Delta = 2000;

previousiloc = zeros(ntrack,1);

122 APPENDICE

previousival = zeros(ntrack,1);

iloc = zeros(ntrack,1);

ival = zeros(ntrack,1);

emergingPicksloc = zeros(ncisoids,1);

emergingPicksval = zeros(ncisoids,1);

minSpacePeaks = 5; % minimum space (bins) between two picked peaks

minSpacePeaks notch = 7; % minimum space (bins)

% nel filtro notch sintetico

first = 1;

fig5 = 'yes';

fig9 = 'yes';

NotchFreq = zeros(ntrack n,LE);

NotchAmp = zeros(ntrack n,LE);

NotchBand = zeros(ntrack n,LE);

PeakFreq = zeros(ntrack,LE);

PeakAmp = zeros(ntrack,LE);

AllFreq = zeros(ncisoids,LE);

% Grafici globali

El = [];

Frequenze = [];

RMagnitude = [];

NMagnitude = [];

%% Query per la ricerca del range di elevazione desiderata

query low = 0;

query up = 90;

[val,ind] = min(abs(E − query low));

j = ind;

%%

% Analisi elevazione per elevazione con azimuth = 0

while (j < LE)

Ie = j;

APPENDICE 123

Onset L = OnL(Ia,Ie)/(fs/1000);

Onset R = OnR(Ia,Ie)/(fs/1000);

onL iT = max(find(T≤Onset L)); % Index del Onset L

onR iT = max(find(T≤Onset R)); % Index del Onset R

hL = squeeze(hrir l(Ia,Ie,onL iT:mxT));

hR = squeeze(hrir r(Ia,Ie,onR iT:mxT));

% Normalizzato [−1, 1]

hL norm = 0.99*hL./max(max(abs(hL)));

hR norm = 0.99*hL./max(max(abs(hR)));

% Ricavo la HRTF

[AL,Fr] = freq resp(hL,Fmin*1000,Fmax*1000,Q,log scale, NFFT, fs);

[AR,Fr] = freq resp(hR,Fmin*1000,Fmax*1000,Q,log scale, NFFT, fs);

LF = length(Fr);

% −−−−−−−−−−−−−−−−−Pre−processing−−−−−−−−−−−−−−−−−[hL w, HL w, AL w, phiL w, Fr] = ...

window(hL norm, Fmin*1000, Fmax*1000, Q, log scale);

[hR w, HR w, AR w, phiR w, Fr] = ...

window(hR norm, Fmin*1000, Fmax*1000, Q, log scale);

% −−−−−−−−−−−−−−−−−−−−Analisi−−−−−−−−−−−−−−−−−−−

% −−−Peack detection−−−Fpeack = [];

nPeaks = 20; % number of peaks detected

nNotch = 20; % number of peaks detected

% ===Left HRTF===

%−−− peak detection

[pftloc w, pftval w] = PickPeaks(AL w,nPeaks,minSpacePeaks);

%−−− notch detection

[nftloc w, nftval w] = PickNotch(AL w,nNotch,minSpacePeaks);

ind = find(nftloc w 6= 0); nftloc w = nftloc w(ind);

124 APPENDICE

nftval w = nftval w(ind);

ind = find(pftloc w 6= 0); pftloc w = pftloc w(ind);

pftval w = pftval w(ind);

%% −−−Calcolo Notch & Contributo Risonanze−−−HRTF = HL w;

A HRTF = AL w;

H multi notch = ones(length(HRTF),1);

nftlocHRTF = nftloc w; nftvalHRTF = nftval w;

pftlocHRTF = pftloc w; pftvalHRTF = pftval w;

iterazioni = 1;

nonotch = 1;

cepstrum BEGIN = [];

cepstrum END = [];

number iterazioni = 300;

soglia = 10ˆ−3;N1 = 4; %−−cutrid = 2;

while (iterazioni ≤ number iterazioni && nonotch == 1)

% 1 − Calcolo Cepstrum

flog = A HRTF; %−− calcolo cepstrum per l'inviluppo

cep = ifft(flog);

cep cut = [cep(1);2*cep(2:N1);cep(N1+1); ...

zeros(length(cep)−N1−1,1)];flog cut = real(fft(cep cut));

if (iterazioni == 1)

cepstrum BEGIN = flog cut;

end

%−−notch[nftloc, nftval] = PickNotch(A HRTF,nNotches,minSpacePeaks);

%−−risonanze

APPENDICE 125

[rftloc, rftval] = PickPeaks(A HRTF,nNotches,minSpacePeaks);

ind = find(nftloc 6= 0); nftloc = nftloc(ind);

nftval = nftval(ind);

ind = find(rftloc 6= 0); rftloc = rftloc(ind);

rftval = rftval(ind);

nftfreq = Fr(nftloc)';

A res db = A HRTF − flog cut;

A res = 10.ˆ(A res db/20);

E res= 10*log10(A res.ˆ2); %−−in db

E res2= A res.ˆ2; %−−no db

[pftlocEres, pftvalEres] = ...

PickPeaks(E res,nPeaks,minSpacePeaks); %−−− peak detection

ind = find(pftlocEres 6= 0); pftlocEres = pftlocEres(ind);

pftvalEres = pftvalEres(ind);

%−−−−

% 4 − Costruzione del filtro multi notch

%Ricerca dei parametri corretti

% Fc − cut off frequencies

% Fb − bandwidth

ampiezza relativa notch = [];

bandwidth = [];

for no = 1 : length(nftval)

amp = A res db(nftloc(no));

% Estremi dell'intervallo di ricerca di fb

target = 0;

nEnergia = E res2(nftloc(no)); %−− energia no in db

nEnergia db = E res(nftloc(no));

if (amp ≤ −5) %−− Caso 1

target = nEnergia db + 3;

126 APPENDICE

else

target = (1 − nEnergia)/2 + nEnergia; %−− no in db,


end

[val, index] = sort(abs(pftlocEres − nftloc(no)));

ind1 = index(1);

ind2 = index(2);

if (pftlocEres(ind1) < pftlocEres(ind2))

e peak sx = E res(pftlocEres(ind1));

e peak dx = E res(pftlocEres(ind2));

else

e peak sx = E res(pftlocEres(ind2));

e peak dx = E res(pftlocEres(ind1));

end

e max peack db = max([e peak sx e peak dx]);

[e min peack db,in] = min([e peak sx e peak dx]);

no left = 1; no right = 1; %−− Caso 3−4if (e min peack db < 0) && (e max peack db > 0)

if (in == 1)

no left = 0;

elseif (in == 2)

no right = 0;

end

end

if (e max peack db < 0) %−− Caso 2

e max peack = 10ˆ(e max peack db/10);

amp max peack db = 20*log10(sqrt(e max peack));

amp rel = amp − amp max peack db;

if (abs(amp rel) ≥ soglia) %−− se false,0db rifer

amp = amp rel;

if (amp ≤ −5)target = amp + 3;

else

APPENDICE 127

target = (e max peack − nEnergia)/2+nEnergia;


end

end

end

loc1 = 0; loc2 = 0;

% Ricerca a destra

trovato = 0;

ind = 1;

while (trovato == 0 && no right == 1)

if ( (E res(nftloc(no)+ind) ≥ target) | | ...

(nftloc(no)+ind == length(E res)) )

trovato = 1;

loc2 = nftloc(no) + ind;

else

ind = ind + 1;

end

end

% Ricerca a sinistra

trovato = 0;

ind = −1;while (trovato == 0 && no left == 1)

if ( (E res(nftloc(no)+ind) ≥ target) | | ...

(nftloc(no)+ind == 1))

trovato = 1;

loc1 = nftloc(no) + ind;

else

ind = ind − 1;

end

end

fc = Fr(nftloc(no));

if (no right == 0)

fb = (fc − Fr(loc1))*2;

elseif (no left == 0)

fb = (Fr(loc2) − fc)*2;

128 APPENDICE

else

fb = Fr(loc2) − Fr(loc1);

end

bandwidth = [bandwidth; fb];

ampiezza relativa notch = [ampiezza relativa notch; amp];

end

ind = find(abs(ampiezza relativa notch) ≥ soglia);

ampiezza relativa notch = ampiezza relativa notch(ind);

nftfreq = nftfreq(ind);

if (length(ampiezza relativa notch) 6= 0)

for k = 1 : size(ampiezza relativa notch,1)

[B,A] = peaking2(ampiezza relativa notch(k), ...

nftfreq(k), bandwidth(k)/rid, fs);

if (k == 1)

b tot = B;

a tot = A;

else

a tot = conv(a tot,A);

b tot = conv(b tot,B);

end

end

[H tot,f tot] = freqz(b tot,a tot,NFFT/2,fs);

A tot = 20*log10(abs(H tot));

ind = find(f tot == Fr(1));

A tot = A tot(ind:end);

H tot = H tot(ind:end);

A riso = A HRTF(1:length(A tot)) − A tot;

H riso = HRTF./H tot;

H multi notch = H multi notch.*H tot;

%A riso pr = 20*log10(abs(H riso)); %−−Controprova% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−% ===Aggiornamento variabili

APPENDICE 129

iterazioni = iterazioni + 1;

HRTF = H riso;

A HRTF = A riso;


[nftlocHRTF, nftvalHRTF] = ...

PickNotch(A HRTF,nNotch,minSpacePeaks);

ind = find(nftlocHRTF 6= 0); nftlocHRTF = nftlocHRTF(ind);

nftvalHRTF = nftvalHRTF(ind);

cepstrum END = flog cut;

% ===

else

nonotch = 0;

end

end

%%

%A Notch = AL w − A HRTF;

A Notch = 20*log10(abs(H multi notch));


[nftlocNotch, nftvalNotch] = ...

PickNotch(A Notch,nNotch,minSpacePeaks notch);

ind = find(nftlocNotch 6= 0); nftlocNotch = nftlocNotch(ind);

nftvalNotch = nftvalNotch(ind);

ind = find(nftlocNotch ≤ 160); nftlocNotch = nftlocNotch(ind);

nftvalNotch = nftvalNotch(ind);

%% −−− Ricostruzione − HRTF sintetica−−−

% −−−Filtro Risonanze−−−plot6 = figure(6);

set(plot6,'Position',[100,100,800,600]);

HRTF2 = 10.ˆ(A HRTF/20); %−−no in db, A HRTF contributo risonanze

HRTF simmetrica = HRTF2;

for sim = 1:length(HRTF2)

HRTF simmetrica = [HRTF simmetrica; ...

HRTF2(length(HRTF) + 1 − sim)];

end

130 APPENDICE

HRIR = real(ifft(HRTF simmetrica));

subplot(2,2,1);

plot(Fr, A HRTF);

grid on; zoom on;

axis([Fr(1) Fr(end) dBmin dBmax]);

title('Left HRTF − Contributo risonanza');

xlabel('f (Hz)');

ylabel('dB');

%−−−Steiglitz−McBride iteration method−−−NB = 0;

NA = 6;

[b,a] = stmcb(HRIR,NB,NA); % b zeri,a poli dell' AR model

ra = flipud(sort(roots(a))); % sorts roots per radius

[H pol,f pol] = freqz(b,a,NFFT/2,fs);

ind = find(f pol == Fr(1));

H pol = H pol(ind:end);

H pol simmetrica = H pol;

for sim = 1:length(H pol)

H pol simmetrica = [H pol simmetrica; ...

H pol(length(H pol) + 1 − sim)];

end

H ris res = HRTF simmetrica./H pol simmetrica;

HRIR res = real(ifft(H ris res));

alp = lpc(HRIR res,6);

rb = flipud(sort(roots(alp)));

[H lp,f lp] = freqz([1],alp,NFFT/2,fs);

ind = find(f lp == Fr(1));

H lp = H lp(ind:end);

H Rison sint = H lp.* H pol;

subplot(2,2,3);

%[H Rison sint,f Rison] = freqz(b,a,NFFT/2,fs);

%ind = find(f Rison == Fr(1));

%H Rison sint = H Rison sint(ind:end);

APPENDICE 131

A Rison = 20*log10(abs(H Rison sint));

zplane(rb,ra);

%Zeri(j) = rb; % Zeri

p = conj(ra(1:NA)); % Poli

%z = conj(rb(1:NB));

%−−−

% −−−Filtro Notch−−−a notch = []; b notch = [];

ampiezza = [];

bandwidth notch = [];

A Notch2 = 10.ˆ(A Notch/20);

E Notch2 = A Notch2.ˆ2; %−−no db

E Notch = 10*log10(E Notch2); %−−in db

for no = 1 : length(nftvalNotch)

% Estremi dell'intervallo di ricerca di fb

target = 0;

n energia = E Notch2(nftlocNotch(no)); %−−no in db

if (nftvalNotch(no) ≤ −5)target = E Notch(nftlocNotch(no)) + 3;

else

%−− energia raddoppiata, no in db

target = (1 − n energia)/2 + n energia;


end

loc1 = 0; loc2 = 0;

% Ricerca a destra

trovato = 0;

ind = 1;

tmp dx = E Notch(nftlocNotch(no));

dx = 0;

132 APPENDICE

while trovato == 0

if ( (E Notch(nftlocNotch(no)+ind) ≥ target) | | ...

(nftlocNotch(no)+ind == length(E Notch)))

trovato = 1;

loc2 = nftlocNotch(no) + ind;

elseif (tmp dx > E Notch(nftlocNotch(no)+ind))

trovato = 1;

loc2 = nftlocNotch(no) + ind − 1;

dx = 1;

else

tmp dx = E Notch(nftlocNotch(no)+ind);

ind = ind + 1;

end

end

% Ricerca a sinistra

trovato = 0;

ind = −1;tmp sx = E Notch(nftlocNotch(no));

sx = 0;

while trovato == 0

if ( (E Notch(nftlocNotch(no)+ind) ≥ target) | | ...

(nftlocNotch(no)+ind == 1) )

trovato = 1;

loc1 = nftlocNotch(no) + ind;

elseif (tmp sx > E Notch(nftlocNotch(no)+ind))

trovato = 1;

loc1 = nftlocNotch(no) + ind + 1;

sx = 1;

else

tmp sx = E Notch(nftlocNotch(no)+ind);

ind = ind − 1;

end

end

f2 = Fr(loc2);

f1 = Fr(loc1);

APPENDICE 133

fc = Fr(nftlocNotch(no));

amp = nftvalNotch(no);

if ((dx == 1) && (sx == 0))

f1 = Fr(loc1);

f2 = fc + (fc − f1);

elseif ((sx == 1) && (dx == 0))

f2 = Fr(loc2);

f1 = fc − (f2 − fc);

elseif ((sx == 1) && (dx == 1)) %−− notch poco pronunciato

min estremo = min([A Notch(loc1) A Notch(loc2)]);

amp = nftvalNotch(no) − min estremo;

end

banda dx = f2 − fc;

banda sx = fc − f1;

fb = f2 − f1;

banda min = min([banda dx banda sx]);

if (2*banda min < fb)

fb = 2*banda min;

end

bandwidth notch = [bandwidth notch; fb];

ampiezza = [ampiezza; amp];

end

nftfreq = Fr(nftlocNotch)';

for k = 1 : length(nftvalNotch)

[B,A] = peaking2(ampiezza(k), nftfreq(k), ...

bandwidth notch(k), fs);

if (k == 1)

b notch = B;

a notch = A;

else

a notch = conv(a notch,A);

134 APPENDICE

b notch = conv(b notch,B);

end

end

[H Notch sint,f Notch sint] = freqz(b notch,a notch,NFFT/2,fs);

A Notch sint = 20*log10(abs(H Notch sint));

ind = find(f Notch sint == Fr(1));

A Notch sint = A Notch sint(ind:end);

H Notch sint = H Notch sint(ind:end);

rb = flipud(sort(roots(b notch)));

ra = flipud(sort(roots(a notch))); % sorts roots per radius

% −−−

% −−−Filtro Sintetico−−−H sint = H Notch sint.*H Rison sint;

A sint = 20*log10(abs(H sint));

%−−−

%% % −−−Notch Tracking nello spazio−−−

nftbandNotch = bandwidth notch/2;

if (first==1) %for the first spatial frame

nNewNotches = 1;

padding = zeros(1,nNotches−length(nftlocNotch));emergingnloc = [nftlocNotch; padding'];

emergingnval = [nftvalNotch; padding'];

emergingnband = [nftbandNotch; padding'];

else

%creating new born tracks

for id=1:nNotches

if (emergingnloc(id) 6= 0)

previousnloc(nNewNotches) = emergingnloc(id);

previousnval(nNewNotches) = emergingnval(id);

previousnband(nNewNotches) = emergingnband(id);

nNewNotches = nNewNotches + 1;

APPENDICE 135

end

end

%Peak tracker

[nloc,nval,nband,emergingnloc,emergingnval,emergingnband, ...

relation]=notchTrack(Delta n,ntrack n,nNotches,nftlocNotch, ...

nftvalNotch,nftbandNotch,previousnloc,previousnval, ...

previousnband);

end

for par = 1 : length(nloc)

if (nloc(par) == 0 | | nloc(par) == −1)NotchFreq(par,j) = 0.;

NotchAmp(par,j) = 0.;

NotchBand(par,j) = 0.;

else

if (j > 1) && (NotchFreq(par,j−1) 6= 0)

NotchFreq(par,j−1) = Fr(previousnloc(relation(par)));

NotchAmp(par,j−1) = previousnval(relation(par));

NotchBand(par,j−1) = previousnband(relation(par));

end

NotchFreq(par,j) = Fr(nloc(par));

NotchAmp(par,j) = nval(par);

NotchBand(par,j) = nband(par);

end

end

if(exist('fig5')) % plot: the trackings of partials

figure(5); clf; hold on

PlotTracking(NotchFreq(:,1:j),NotchAmp(:,1:j),[], ...

NotchBand(:,1:j));

%PlotTracking(NotchFreq(:,1:j),NotchAmp(:,1:j));

xlabel('Elevation');ylabel('Frequency (Hz)');

title('Notch tracking');

grid on;

axis([E(1) 90 0 Fr(end)]); drawnow

end

136 APPENDICE

previousnval = nval;

previousnloc = nloc;

%% %−−−Peack Tracking nello spazio−−−f = −angle to freq(angle(p),fs);

t = −1./(log(abs(p)).*fs);

[f,inds]=sort(f);

t=t(inds);

p=p(inds);

f = f(find(f>0));

iftloc = ones(length(f),1);

iftval = ones(length(f),1);

for r = 1:length(f)

for l = 1 : length(Fr)

if Fr(l) ≥ f(r)

bin = l;

break;

end

end

iftloc(r,1) = bin;

iftval(r,1) = A Rison(bin);

end

if (first==1) %for the first spatial frame

nNewPeaks = 1;

padding = zeros(1,ncisoids−length(f));emergingPicksloc = [iftloc; padding'];

emergingPicksval = [iftval; padding'];

else

%creating new born tracks

for id=1:ncisoids

if (emergingPicksloc(id) 6= 0)

previousiloc(nNewPeaks) = emergingPicksloc(id);

previousival(nNewPeaks) = emergingPicksval(id);

nNewPeaks = nNewPeaks + 1;

end

APPENDICE 137

end

%Peak tracker

[iloc,ival,iphase,emergingPicksloc,emergingPicksval, ...

relation] = peakTrackSimple(Delta,ntrack,ncisoids, ...

iftloc,iftval,[],previousiloc,previousival);

end

for par = 1 : length(iloc)

if (iloc(par) == 0 | | iloc(par) == −1)PeakFreq(par,j) = 0.; % frequency of the partials,

PeakAmp(par,j) = 0.;

else

if (j > 1) && (PeakFreq(par,j−1) 6= 0)

PeakFreq(par,j−1) = Fr(previousiloc(relation(par)));

PeakAmp(par,j−1) = previousival(relation(par));

end

PeakFreq(par,j) = Fr(iloc(par));

PeakAmp(par,j) = ival(par);

end

end

AllFreq(:,j)=Fr(iftloc);

if(exist('fig9')) % plot: the trackings of partials

figure(9); clf; hold on

PlotTracking(PeakFreq(:,1:j),PeakAmp(:,1:j),AllFreq(:,1:j));

%PlotTracking(PeakFreq(:,1:j),PeakAmp(:,1:j));

xlabel('Elevation');ylabel('Frequency (Hz)');

%axis([E(1) 90 0 Fmax]);

grid on;

title('Peak tracking'); drawnow

end

previousival = ival;

previousiloc = iloc;

%−−−

%% −−−−−−−−−−−−−−−−−−−−−

138 APPENDICE

% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−% Aggiornamento per grafici globali

El = [El E(j)];

Frequenze = [Frequenze Fr'];

RMagnitude = [RMagnitude A HRTF];

NMagnitude = [NMagnitude A Notch];

first = 0; %−−eseguite le operazioni di inizializzazione

% =======================================

% Ricerca Upper bound della query

if ( (j 6= LE) && (E(j+1) > query up))

j = LE;

else

j=j+1;

end

% =======================================

end

Elevazione = [];

for row = 1 : size(Frequenze,1)

Elevazione = [Elevazione; El];

end

figure(10+individuo);

view([145 30]);

surfc(Elevazione,Frequenze,RMagnitude);

caxis([−60 30]);

shading flat;

title(['Risonanze LEFT−hrir − ' 'Subject: ' subject string]);

xlabel('elev (gradi)');

ylabel('f (kHz)');

zlabel('magnitude');

colormap('jet');

colorbar;

figure(11+individuo);

APPENDICE 139

view([145 30]);

surfc(Elevazione,Frequenze,NMagnitude);

caxis([−40 0]);

shading flat;

title(['Notches LEFT−hrir − ' 'Subject: ' subject string]);

xlabel('elev (gradi)');

ylabel('f (hz)');

zlabel('magnitude');

colormap('jet');

colorbar;

individuo = individuo + 2;

end

end

%% −−−−−Analysis Functions−−−−−

function [loc, val] = PickPeaks(spectrum, nPeaks, minspace)

%

%==> peaking the nPicks highest peaks in the given spectrum

% from the greater to the lowest

% OUT

% loc: bin number of peaks (if loc(i)==0, no peak detected)

% val: amplitude of the given spectrum

% IN

% spectrum: spectrum (abs(fft(signal))

% nPicks: number of peaks to pick

% minspace: minimum of space between two peaks

[r, c] = size(spectrum);

rmin = min(spectrum) − 1;

% −−−find a peak, zero out the data around the peak, and repeat

val = ones(nPeaks,c)*−100;loc = zeros(nPeaks,c);

for k=1:c %−−− find all local peaks

difference = diff([rmin; spectrum(:,k); rmin]); % derivate

iloc = find(difference(1:r)≥ 0 & difference(2:r+1) ≤ 0);

% peak locations

140 APPENDICE

ival = spectrum(iloc,k); % peak values

for p=1:nPeaks

[val(p,k),l] = max(ival); % find current maximum

loc(p,k) = iloc(l); % save value and location

ind = find(abs(iloc(l)−iloc) > minspace);

% find peaks which are far away

if (isempty(ind))

break % no more local peaks to pick

end

ival = ival(ind); % shrink peak value and location array

iloc = iloc(ind);

end

end

end %end function

function [loc, val] = PickNotch(spectrum, nNotch, minspace)

%

%==> peaking the nNotch deepest notch in the given spectrum

% from the greater to the lowest

% OUT

% loc: bin number of notch (if loc(i)==0, no notch detected)

% val: amplitude of the given spectrum

% IN

% spectrum: spectrum (abs(fft(signal))

% nNotch: number of notch to pick

% minspace: minimum of space between two notch

[r, c] = size(spectrum);

rmin = min(spectrum) − 1;

% −−−find a notch, zero out the data around the peak, and repeat

val = ones(nNotch,c)*−100;loc = zeros(nNotch,c);

for k=1:c %−−− find all local peaks

difference = diff([rmin; spectrum(:,k); rmin]); % derivate

iloc = find(difference(1:r)≤ 0 & difference(2:r+1) ≥ 0);

% notch locations

ival = spectrum(iloc,k); % notch values

for p=1:nNotch

APPENDICE 141

[val(p,k),l] = min(ival); % find current minimum

loc(p,k) = iloc(l); % save value and location

ind = find(abs(iloc(l)−iloc) > minspace);

% find notch which are far away

if (isempty(ind))

break % no more local notch to pick

end

ival = ival(ind); % shrink notch value and location array

iloc = iloc(ind);

end

end

end %end function

%% −−−−− Partial Tracking Function −−−−−

function[iloc,ival,iphase,emergingPicksloc,emergingPicksval,relation]= ...

peakTrackSimple(Delta,nSines,nPeaks,iftloc,iftval,iftphase, ...

previousiloc,previousival)

%

%==> simplest partial tracking

% OUT:

% iloc,ival,iphase: location (bin), magnitude and phase of peaks

% emergingPicksloc,emergingPicksval: loc e magnitude freq non tracciate

% relation: relazioni con il frame erecedente

% IN:

% previousiloc,previousival: info previous frame

% iftloc, iftval, iftphase: info picchi correnti

% nPeaks: number of peaks detected

% nSines: number of peaks tracked

% Delta: ∆ per il tracking delle frequenze

MindB = −100;NFFT = 200;

fs = 44100;

iloc = zeros(nSines,1);

ival = zeros(nSines,1) + MindB;

iphase = zeros(nSines,1);

emergingPicksloc = zeros(nPeaks,1);

142 APPENDICE

emergingPicksval = zeros(nPeaks,1);

relation = zeros(nSines,1);

distminindex = zeros(nSines,1);

for i=1:nSines %−−− for each sinus detected

if (previousiloc(i) 6= 0)

[closestpeakmag,closestpeakindex]=min(abs(iftloc−previousiloc(i)));if (closestpeakmag/NFFT*fs ≤ Delta) %−−− soglia

iloc(i) = iftloc(closestpeakindex); %−−− bin of the closest

ival(i) = iftval(closestpeakindex); %−−− value of the closest

%iphase(i) = iftphase(closestpeakindex);

prev=closestpeakindex − 1;

next=closestpeakindex + 1;

iftloc = [iftloc(1:prev); iftloc(next:end)];

iftval = [iftval(1:prev); iftval(next:end)];

dist = abs(previousiloc−iloc(i));[distminval, distminindex(i)] = min(dist);

relation(i) = i; %−−mantiene la relazione con dati precedenti

else

iloc(i) = −1; %−−escludo la traccia dal venir riutilizzata

ival(i) = MindB;

relation(i) = −1;end

end

end

for j=1:length(iftloc)

emergingPicksloc(j) = iftloc(j);

emergingPicksval(j) = iftval(j);

end

end %end function

function[iloc,ival,iband,emergingPicksloc,emergingPicksval, ...

emergingPicksband,relation]=notchTrack(Delta,nSines,nPeaks, ...

iftloc,iftval,iftband,previousiloc,previousival,previousiband)

APPENDICE 143

MindB = −100;NFFT = 200;

fs = 44100;

iloc = zeros(nSines,1);

ival = zeros(nSines,1) + MindB;

iband = zeros(nSines,1);

emergingPicksloc = zeros(nPeaks,1);

emergingPicksval = zeros(nPeaks,1);

emergingPicksband = zeros(nPeaks,1);

relation = zeros(nSines,1);

distminindex = zeros(nSines,1);

for i=1:nSines %−−− for each sinus detected

if (previousiloc(i) 6= 0)

[closestpeakmag,closestpeakindex]=min(abs(iftloc−previousiloc(i)));if (closestpeakmag/NFFT*fs ≤ Delta) %−−− soglia

iloc(i) = iftloc(closestpeakindex); %−−− bin of the closest

ival(i) = iftval(closestpeakindex); %−−− value of the closest

iband(i) = iftband(closestpeakindex); %−−− value of the closest

%iphase(i) = iftphase(closestpeakindex);

prev=closestpeakindex − 1;

next=closestpeakindex + 1;

iftloc = [iftloc(1:prev); iftloc(next:end)];

iftval = [iftval(1:prev); iftval(next:end)];

iftband = [iftband(1:prev); iftband(next:end)];

dist = abs(previousiloc−iloc(i));[distminval, distminindex(i)] = min(dist);

relation(i) = i; %−−mantiene la relazione con i dati precedenti

else

iloc(i) = −1; %−−escludo la traccia dal venir riutilizzata

ival(i) = MindB;

iband(i) = 0;

relation(i) = −1;end

144 APPENDICE

end

end

for j=1:length(iftloc)

emergingPicksloc(j) = iftloc(j);

emergingPicksval(j) = iftval(j);

emergingPicksband(j) = iftband(j);

end

end %end function

function PlotTracking(SineFreq,AllAmp,AllFreq,AllBand)

%

%==> plot the partial tracking

% IN

% SineFreq: Frequenze delle tracce

% AllAmp: Ampiezza delle tracce

% AllBand: Bandapassante delle tracce

hold on;

E = −45:(360/64):235;Fmax = 21.93975;

[nSines, nFrames] = size(SineFreq);

for n=1:nSines

f=1;

while (f≤nFrames)

while (f≤nFrames & SineFreq(n,f)==0)

f = f+1;

end

iStart = min(f,nFrames);

while (f≤nFrames & SineFreq(n,f)>0)

f = f+1;

end

iEnd = min(max(1,f−1),nFrames);if (iEnd > iStart)

if (nargin == 4)

errorbar(E(iStart:iEnd), SineFreq(n,iStart:iEnd), ...

AllBand(n,iStart:iEnd));

APPENDICE 145

elseif (nargin == 3)

plot(E(iStart:iEnd), SineFreq(n,iStart:iEnd));

z = AllAmp(n,iStart:iEnd);

scatter(E(iStart:iEnd),SineFreq(n,iStart:iEnd),50,z,...

'filled');

caxis([0 20]);

%axis([−50 95 0 Fmax]);

colorbar;

elseif (nargin == 2)

plot(E(iStart:iEnd), SineFreq(n,iStart:iEnd));

z = AllAmp(n,iStart:iEnd);

scatter(E(iStart:iEnd),SineFreq(n,iStart:iEnd), ...

50,z,'filled');

caxis([−40 0]);

colorbar;

end

end

end

end

if (nargin == 3)

[nSines, nFrames] = size(AllFreq);

frame = [1:1:nFrames];

frame tot = [];

for n=1:nSines

frame tot = [frame tot; frame];

end

%plot(frame tot,AllFreq,'rs');

end;

hold off;

end %end function

function [b, a] = peaking2(G, fc, fb, fs)

d = −cos(2*pi*(fc/fs));V0 = 10ˆ(G/20);

146 APPENDICE

H0 = V0 − 1;

if(G ≤ 0)

a = (tan(pi*fb/fs) − V0)/(tan(pi*fb/fs) + V0);

else

a = (tan(pi*fb/fs) − 1)/(tan(pi*fb/fs) + 1);

end

b0 = 1 + (1 + a)*(H0/2);

b1 = d*(1 − a);

b2 = −a − (1 + a)*(H0/2);

a1 = d*(1 − a);

a2 = −a;

%return values

a = [ 1, a1, a2];

b = [ b0, b1, b2];

function [x w, H w, A w, phi w, Fr] = window(x, Fmin, Fmax, Q, log scale)

% Finestratura

% Filtraggio con una half−Hann window

NFFT = 400; %2*LT

fs = 44100;

L = fs/1000; % sample in 1 ms

w = hann(2*L);

w half = w((L+1):(length(w)+1));

w half l = [w half; zeros( (length(x)−length(w half)), 1 )];

x w = x.*w half l;

% Calcolo delle HRTF dei segnali finestrati

if (Fmin == 0 && Fmax == Inf)

[H w,A w,phi w,Fr] = freq resp tot(x w,0,Inf,Q,log scale, NFFT, fs);

else

[H w,A w,phi w,Fr] = freq resp tot(x w,Fmin,Fmax,Q,log scale, NFFT, fs);

APPENDICE 147

end

148 APPENDICE

Bibliografia

[1] B. Kapralos, M. Jenkin and E. Milios, 2008 Virtual audio systems,

Presence: Teleoperators and Virtual Environments, 17(6):527-549.

[2] Eric J. Angel, V. Ralph Algazi, and Richard O. Duda, 2002 On the desi-

gn of canonical sound localization environments, Presented at the 113th

Convention, California, USA.

[3] V. Ralph Algazi, Richard O. Duda, Dennis M. Thompson and Carlos

Avendano, 2001. The CIPIC HRTF Database, IEEE Workshop on Ap-

plications of Signal Processing to Audio and Acoustics, New Paltz, New

York, W2001 1-4.

[4] V. Ralph Algazi, Richard O. Duda, Dennis M. Thompson and Reed P.

Morrison, 2001. Structural Composition and Decomposition of HRTFS,

IEEE Workshop on Applications of Signal Processing to Audio and

Acoustics, New Paltz, New York, W2001 1-4.

[5] CIPIC Interface Laboratory, 1998. Documentation for UCD HRIR Files,

University of California at Davis.

[6] R. Roy, A. Paulraj, and T. Kailath, 1986 ESPRIT-A Subspace Rotation

Approach to Estimation of Parameters of Cisoids in Noise, IEEE Tran-

sactions on Acoustics, Speech, and Signal Processing, vol. ASSP-34, no.

5.

149

150 BIBLIOGRAFIA

[7] Jean Laroche, 1993. The use of the matrix pencil method for the spectrum

analysis of musical signals, Journal of the Acoustical Society of America.

94 (4).

[8] Roland Badeau, Gael Richard and Bertrand David, 2005. Fast Adapti-

ve ESPRIT Algorithm, 13th Workshop on Statistical Signal Processing,

IEEE/SP.

[9] Paulo A. A. Esquef, Matti Karjalainen and Vesa Valimaki, 2003.

Frequency-Zooming ARMA Modeling for Analysis of Noisy String Instru-

ment Tones, EURASIP Journal on Applied Signal Processing,953-967.

[10] Hanna Jarvelainen and Matti Karjalainen, 2005. Importance of Inhar-

monicity in the Acoustic Guitar, Proceedings Int. Computer Music Conf.,

pp.363-366, Barcelona, Spain, September 5-9.

[11] Matti Karjalainen, Vesa Valimaki and Paulo A. A. Esquef, , 2002. Ef-

ficient Modeling and Synthesis of Bell-like sounds, Proc. of the 5th Int.

Conferenze on Digital Audio Effects (DAFx-02), Hamburg, Germany.

[12] Heidi-Maria Lehtonen Analysis of Piano Tones Using an Inharmonic

Inverse Comb Filter, Proc. of the 11th Int. Conferenze on Digital Audio

Effects (DAFx-08), Espoo, Finland.

[13] Carlos Avendano, Richard O. Duda and V. Ralph Algazi, 1999. Mode-

ling the Contralateral HRTF, AES 16th International Conference Spatial

Sound Reproduction, Helsinki, Finland.

[14] Vikas C. Raykar, Ramani Duraiswami, Larry Davis and B. Yegnanara-

yana, 2003. Extracting Significant Features from the HRTF, Proceedings

of the 2003 International Conference on Auditory Display, Boston, MA,

USA.

[15] Robert J. McAulay and Thomas F. Quartieri, 1986. Speech Analysis/-

Synthesis Based on a Sinusoidal Representation, IEEE Transactions on

ACoustics, Speech, and Signal Processing, vol. ASSP-34, no. 4.

BIBLIOGRAFIA 151

[16] Mathieu Lagrange, Sylvain Marchand, Martin Raspaud, and Jean-

Bernard Rault, 2003. Enhanced Partial Tracking Using Linear Predic-

tion, Proc. of the 6th Int. Conference on Digital Audio Effects (DAFx-03),

London, UK.

[17] Corey Kereliuk, and Philippe Depalle, 2008 Improved Hidden Markov

Model Partial Tracking Through Time-frequency Analysis, Proc. of the

11th Int. Conference on Digital Audio Effects (DAFx-08), Espoo, Finland.

[18] Martin Raspaud and Gianpaolo Evangelista, 2008 Binaural Partial

Tracking, Proc. of the 11th Int. Conference on Digital Audio Effects

(DAFx-08), Espoo, Finland.

[19] Mathieu Lagrange, Bertrand Scherrer, 2008. TWO-STEP Modal Iden-

tification for Increased Resolution Analysis of Percussive Sounds, Proc.

of the 11th Int. Conference on Digital Audio Effects (DAFx-08), Espoo,

Finland.

[20] Julius O. Smith III, and Xavier Serra, 1987. PARSHL: An Analysis/Syn-

thesis Program for Non-Harmonic Sounds Based on a Sinusoidal Repre-

sentation, Proceedings of the International Computer Music Conference

(ICMC-87).

[21] Patrick Satarzadeh, V. Ralph Algazi, and Richard O. Duda, 2007. Phy-

sical and Filter Pinna Models Based on Anthropometry, Presented at the

122nd AES Convention, Vienna, Austria.

[22] Patrick Satarzadeh, V. Ralph Algazi, and Richard O. Duda, 2006. A

Study of Physical and Circuit Models of Human Pinnae, Master of Science

Thesis, University of California Davis.

[23] E.A.G Shaw, 1997. Acoustical features of human ear, In Robert H.

Gilkey and Timothy R. Anderson, editors, Binaural and Spatial Hea-

ring in Real and Virtual Environments pag 25-47, Lawrence Erlbaum

Associates,Mahwah,NJ.

152 BIBLIOGRAFIA

[24] C. Phillip Brown and Richard O. Duda, 1998. A Structural Model

for Binaural Sound Synthesis, IEEE Transactions on Speech and Audio

Processing, Vol 6, No. 5, 476-488.

[25] Vikas C. Raykar, and Ramani Duraiswami, 2005. Extracting the fre-

quencies of the pinna spectral notches in measured head related impulse

responses, J. Acoust. Soc. Am. 118 (1), July 2005

[26] Enrique A. Lopez-Poveda and Ray Meddis, 1996. A physical model of

sound diffraction and reflections in the human concha, J. Acoust. Soc.

Am. 100 (5), November 1996

[27] Sergio G. Rodrıguez, and Miguel A.Ramırez, 2005. Extracting and Mo-

deling Approximated Pinna-related Transfer Functions from HRTF Data,

Proceedings of ICAD 05th Meeting of the International Conference on

Auditory Display, Limerick, Ireland

[28] V. Ralph Algazi, Richard O. Duda, Reed P. Morrison and Dennis M.

Thompson, 2001. Structural Composition and Decomposition of HRTFs,

IEEE Workshop on Applications of Signal Processing to Audio and

Acoustics, New Paltz, New York, 103-106.

[29] V. Ralph Algazi, Carlos Avendano and Richard O. Duda, 2001. Estima-

tion of a Spherical-Head Model from Anthropometry, JAES Volume 49,

Issue 6, 472-479.

[30] Richard O. Duda, Carlos Avendano, and V. Ralph Algazi, 1999. An

Adaptable Ellipsoidal Head Model for the Interaural Time Difference,

Proc. IEEE Int. Conf. Acoustics Speech and Signal Processing, Phoenix,

AZ, 965-968.

[31] V. Ralph Algazi, Carlos Avendano and Richard O. Duda, 2000. Ele-

vation localization and head-related transfer function analysis at low

frequencies, J. Acoust. Soc. Am. 109 (3), March 2001.

BIBLIOGRAFIA 153

[32] V. Ralph Algazi, Richard O. Duda, and Dennis M. Thompson 2002.

The Use of Head-and-Torso Models for Improved Spatial Sound Synthesis,

AES 113th Convention, Los Angeles, CA, USA, 2002.

[33] V. Ralph Algazi, Richard O. Duda, Ramani Duraiswami, Nail A. Gu-

merov, and Zhihui Tang, 2002. Approximating the head-related transfer

function using simple geometric models of the head and torso, J. Acoust.

Soc. Am. 112 (5), Pt. 1, Nov. 2002.

[34] Kenneth John Faller II, Armando Barreto, Navarun Gupta and Naph-

tali Rishe, 2005. Enhanced Modeling of Head-Related Impulse Responses

Towards the Development of Customizable Sound Spatialization, Procee-

dings of the WSEAS / IASME International Conference on Computa-

tion Intelligence, Man-Machine System and Cybernetics (CIMMACS’05),

Miami, Florida, USA.

[35] Kenneth John Faller II, Armando Barreto, Navarun Gupta and Naphtali

Rishe, 2006. Time and Frequency Decomposition of Head-Related Impulse

Responses for the Development of Customizable Spatial Audio Models,

WSEAS Transactions on Signal Processing, Vol. 2 (11), pp. 1465-1472.

[36] Kenneth John Faller II, Armando Barreto, Navarun Gupta and Naphtali

Rishe, 2006. Accelerated Method for the Reduced-Parameter Modeling of

Head-Related Transfer Function for Customizable Spatial Audio, Procee-

dings of the 5th WSEAS International COnference on Circuits, System,

Electronics, Control and Signal Processing (CSECS’06), Dallas, Texas,

USA.

[37] J. Huopaniemi, 1999. Virtual Acoustics and 3-D Sound in Multime-

dia Signal Processing, Thesis for the degree of Doctor of Science in

Technology, Helsinki University of Technology, Faculty of Electrical and

Communications Engineering, Laboratory of Acoustics and Audio Signal

Processing.

154 BIBLIOGRAFIA

[38] F. Avanzini, 2008. Sound in Space, in Algorithms for Sound and Music

Computing, for Computer Science class in Sound and Music Computing.

[39] Durand R. Begault, 2000. 3-D Sound for Virtual Reality and Multimedia,

Ames Research Center, Moffett Field, California.

[40] Edited by Udo Zolzer, 2002. Digital Audio Effects, Copyright ©2002 by

John Wiley and Sons, Ltd.

[41] Udo Zolzer, 2008. Digital Audio Signal Processing, Second Edition, John

Wiley and Sons, Ltd.

MODELLO STRUTTURALE DI ORECCHIO ESTERNO PER IL...

Documents

Transcript of MODELLO STRUTTURALE DI ORECCHIO ESTERNO PER IL...