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.
8 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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
10 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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.
12 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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:
1.2. LOCALIZZAZIONE DEL SUONO: L’AUDIO BINAURALE 13
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.
14 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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:
1.2. LOCALIZZAZIONE DEL SUONO: L’AUDIO BINAURALE 15
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.
16 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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 + θ
1.2. LOCALIZZAZIONE DEL SUONO: L’AUDIO BINAURALE 17
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.
18 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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.
1.2. LOCALIZZAZIONE DEL SUONO: L’AUDIO BINAURALE 19
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
20 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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
1.2. LOCALIZZAZIONE DEL SUONO: L’AUDIO BINAURALE 21
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
22 CAPITOLO 1. BACKGROUND E STATO DELL’ARTE
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.
26 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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)
2.1. HRIR E HRTF PERSONALI 27
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.
28 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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)
2.1. HRIR E HRTF PERSONALI 29
∇ ·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.
30 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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;
2.1. HRIR E HRTF PERSONALI 31
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.
32 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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-
2.1. HRIR E HRTF PERSONALI 33
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.
34 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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)
2.1. HRIR E HRTF PERSONALI 35
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, λ)
36 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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;
38 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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-
2.2. MODELLI STRUTTURALI 39
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)
40 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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.
2.2. MODELLI STRUTTURALI 41
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-
42 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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.
44 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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
2.3. PRTF E ANTROPOMETRIA 45
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.
46 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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].
2.3. PRTF E ANTROPOMETRIA 47
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)
48 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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;
2.3. PRTF E ANTROPOMETRIA 49
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;
50 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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.
2.3. PRTF E ANTROPOMETRIA 51
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
52 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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)
2.3. PRTF E ANTROPOMETRIA 53
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;
54 CAPITOLO 2. MODELLAZIONE HRIR E HRTF
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.
58 CAPITOLO 3. METODOLOGIE E STRUMENTI
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
60 CAPITOLO 3. METODOLOGIE E STRUMENTI
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
62 CAPITOLO 3. METODOLOGIE E STRUMENTI
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.
3.3. MODELLI SINUSOIDALI 63
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.
64 CAPITOLO 3. METODOLOGIE E STRUMENTI
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.
3.3. MODELLI SINUSOIDALI 65
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
66 CAPITOLO 3. METODOLOGIE E STRUMENTI
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.
68 CAPITOLO 3. METODOLOGIE E STRUMENTI
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.
3.4. CIPIC DATABASE 69
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.
70 CAPITOLO 3. METODOLOGIE E STRUMENTI
Figura 3.7: Parametri antropometrici.
3.4. CIPIC DATABASE 71
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.
72 CAPITOLO 3. METODOLOGIE E STRUMENTI
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;
76 CAPITOLO 4. ANALISI DELLA PRTF
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
4.1. ALGORITMO DI ANALISI 77
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
78 CAPITOLO 4. ANALISI DELLA PRTF
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)
4.1. ALGORITMO DI ANALISI 79
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)
80 CAPITOLO 4. ANALISI DELLA PRTF
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.
4.1. ALGORITMO DI ANALISI 81
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
82 CAPITOLO 4. ANALISI DELLA PRTF
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;
target = 10*log10(target);
end
end
end
4.1. ALGORITMO DI ANALISI 83
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.
84 CAPITOLO 4. ANALISI DELLA PRTF
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.
4.1. ALGORITMO DI ANALISI 85
(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.
86 CAPITOLO 4. ANALISI DELLA PRTF
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):
4.1. ALGORITMO DI ANALISI 87
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.
88 CAPITOLO 4. ANALISI DELLA PRTF
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)
Notches LEFT−hrir − Subject: 010
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
90 CAPITOLO 4. ANALISI DELLA PRTF
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;
4.2. ANALISI DEI DATI 91
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (kHz)
Risonanze LEFT−hrir − Subject: 010
−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.
92 CAPITOLO 4. ANALISI DELLA PRTF
gliamo far notare come per il soggetto 165 (testa KEMAR) siano accentuati
i tratti che riguardano la sovrapposizione degli effetti di risonanze vicine.
4.2. ANALISI DEI DATI 93
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Risonanze LEFT−hrir − Subject: 048
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.8x 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.9: Subject 048 left pinna - Contributo delle risonanze in relazione
all’elevazione.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Risonanze LEFT−hrir − Subject: 134
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.8x 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.10: Subject 134 left pinna - Contributo delle risonanze in relazione
all’elevazione.
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (kHz)
Risonanze LEFT−hrir − Subject: 165
−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.8x 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.11: Subject 165 left pinna - Contributo delle risonanze in relazione
all’elevazione.
94 CAPITOLO 4. ANALISI DELLA PRTF
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 048
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.12: Subject 048 left pinna - Contributo dei notch in relazione
all’elevazione.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 134
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.13: Subject 134 left pinna - Contributo dei notch in relazione
all’elevazione.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 165
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.14: Subject 165 left pinna - Contributo dei notch in relazione
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.
98 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA PINNA
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-
5.2. FILTRO SINTETICO 99
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.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;
Subject 048 tra φ = 22 e φ = 62;
Subject 134 tra φ = 40 e φ = 84;
Subject 165 tra φ = 12 e φ = 28;
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.
100 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA 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
(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;
5.2. FILTRO SINTETICO 101
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
102 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA PINNA
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
Left HRTF − Risonanze − Z−plane
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:
5.2. FILTRO SINTETICO 103
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
Left HRTF − Risonanze − Z−plane
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
104 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA PINNA
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
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
(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
Left HRTF − Risonanze − Z−plane
(b) Contributo Notch.
Figura 5.9: Subject 134 left pinna φ = 0, filtri sintetici
5.2. FILTRO SINTETICO 105
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.
106 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA PINNA
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.
5.2. FILTRO SINTETICO 107
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: 165 − Elevation: 79
PRTF originalePRTF sintetica
(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
PRTF originalePRTF sintetica
(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
20Left HRTF − Filtro Notch sintetico Subject: 134 − Elevation: 6
PRTF originalePRTF sintetica
(c)
Figura 5.11: Alcuni esempi di filtri sintetici ottenuti nella fase di
modellazione, messi a confronto con le corrispettive PRTF reali.
108 CAPITOLO 5. UN MODELLO STRUTTURALE PER LA PINNA
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)
Notches LEFT−hrir − Subject: 134
−40
−35
−30
−25
−20
−15
−10
−5
0
(a) contributo notch
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Risonanze LEFT−hrir − Subject: 134
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)
Notches LEFT−hrir − Subject: 134
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)
Risonanze LEFT−hrir − Subject: 134
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.3: Subject 134 left pinna - Nceps = 3, Risultati.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 134
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)
Risonanze LEFT−hrir − Subject: 134
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.4: Subject 134 left pinna - Nceps = 5, Risultati.
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (hz)
Notches LEFT−hrir − Subject: 134
−40
−35
−30
−25
−20
−15
−10
−5
0
(a) contributo notch
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Risonanze LEFT−hrir − Subject: 134
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.5: Subject 134 left pinna - Nceps = 7, Risultati.
A.3. PARAMETRO RID 117
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (hz)
Notches LEFT−hrir − Subject: 165
−40
−35
−30
−25
−20
−15
−10
−5
0
(a) contributo notch
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Risonanze LEFT−hrir − Subject: 165
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.6: Subject 165 left pinna - soglia = 1, Risultati.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 165
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)
Risonanze LEFT−hrir − Subject: 165
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.7: Subject 165 left pinna - soglia = 10−1, Risultati.
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (hz)
Notches LEFT−hrir − Subject: 165
−40
−35
−30
−25
−20
−15
−10
−5
0
(a) contributo notch
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (kHz)
Risonanze LEFT−hrir − Subject: 165
−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.8: Subject 165 left pinna - soglia = 10−2, Risultati.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 165
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)
Risonanze LEFT−hrir − Subject: 165
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.9: Subject 165 left pinna - soglia = 10−6, Risultati.
118 APPENDICE
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 048
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)
Risonanze LEFT−hrir − Subject: 048
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.10: Subject 048 left pinna - rid = 1, Risultati.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 048
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)
Risonanze LEFT−hrir − Subject: 048
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.11: Subject 048 left pinna - rid = 3, Risultati.
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (hz)
Notches LEFT−hrir − Subject: 048
−40
−35
−30
−25
−20
−15
−10
−5
0
(a) contributo notch
−50
0
50
100
00.511.522.5x 10
4
elev (gradi)
f (kHz)
Risonanze LEFT−hrir − Subject: 048
−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.12: Subject 048 left pinna - rid = 4, Risultati.
−50
0
50
100
00.511.522.5
x 104
elev (gradi)
Notches LEFT−hrir − Subject: 048
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)
Risonanze LEFT−hrir − Subject: 048
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.13: Subject 048 left pinna - rid = 8, Risultati.
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
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,
target = 10*log10(target);
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;
target = 10*log10(target);
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;
%−−− notch detection
[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));
%−−− notch detection
[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;
target = 10*log10(target);
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.