Arsén Ivanóv: EEN 502 Project II (Fall 06)

Abstract:

The purpose of the project is to demonstrate understanding and ability to use Head-Related Transfer Function (HRTF) data. To achieve that goal, the project is assembled from three distinct tasks:
  1. Simulate a sound source rotating clockwise around the observer in the azimuth plane (zero degrees elevation relative to the eye level), starting and ending in front of the observer.
  2. Simulate a sound source rotating upward from the point 40 degrees below azimuth plane in fron of the observer to the point 90 degrees above the azimuth plane in fron of the observer.
  3. Simulate the source tyraveling in the strate line at a certain distance from the observer using:
    1. Interaural time and distance-caused intensity attenuation with Doppler effect,
    2. Interaural HRTFs without Doppler effect,
    3. Both interaural HRTF and interaural Doppler effect

Table of Contents:

  1. Simulating a roating sound source in the azimuth plane using HRTFs.
  2. Simulating a roating sound source in the elevation plane using HRTFs.
  3. Combining Doppler effect and HRTFs for interaural simulation of a sound source traveling in a straight line at a distance from the observer:
    1. Interaural simulation using Doppler effect.
    2. Interaural simulation using HRTFs.
    3. Interaural simulation using both Doppler effect and HRTFs.
  4. Conclusion

Simulating a sound source rotating in the azimuth plane using HRTFs.

Objective:

Simulate using HRTFs a source in a circular orbit with the head in the center of the circle. The source moves on the azimuth plane (the elevation angle is zero), with constant speed. Assume that Doppler effects are negligible.

Discussion:

The task requires several steps.
  1. the 360 degrees around the head must be partitioned into chunks of 5 degrees, since the set of HRTFs available are defined only for 5 degrees each.
  2. during the time that source is within a certain 5 degree slice, its sound output must be convolved with the particular HRTF for that spot, separately for each ear.
  3. the resulting vectors must be combined appropriately, compensating for the fact that convolution operation produces more samples than there are in the input sound.
The most difficult part of the design process was actually understanding the format of the HRTF recordings. It became apparent that the HRTFs in the set were actually recorded using two different pinnae - one normal and one enlarged. The enlarged pinna was attached to the left of the KEMAR manequin, while enlarged pinna was attached to the right. The confusion stemmed from the fact that any given HRTF was labeled either 'R' or 'L' for left or right pinnae, referring to the type of the pinna, but themselves the HRTFs were symmetrical, i.e. each HRTF held data simulating the same pinna hearing the same sound on both sides of the manequin's head. Once that confusion was resolved, the process was straightforward.
The code included below was used to generate the simulation, that can be heard here. The code iterates through 36 slices of 5 degrees each on the right hemisphere, and then repeats iteration again, flipping certain parts to account for the continuous motion of the sound source on the other side of the head. 
fs=44100;
elevation=0;
T=10;
num_samples=T*fs;
num_samples_5deg=floor(num_samples/72);
snd=zeros(num_samples,2);
for i=0:36
    hrtf=readhrtf(elevation,i*5,'L');
    src=randn(num_samples_5deg,1);
    R=conv(hrtf(1,:),src);
    L=conv(hrtf(2,:),src);
    [r c]=size(R);
    snd((i*num_samples_5deg+1):(i*num_samples_5deg+r),1)=R;
    snd((i*num_samples_5deg+1):(i*num_samples_5deg+r),2)=L;
end
for i=0:36
    hrtf=readhrtf(elevation,(36-i)*5,'L');
    src=randn(num_samples_5deg,1);
    R=conv(hrtf(1,:),src);
    L=conv(hrtf(2,:),src);
    [r c]=size(R);
    snd(((36+i)*num_samples_5deg+1):((36+i)*num_samples_5deg+r),1)=L;
    snd(((36+i)*num_samples_5deg+1):((36+i)*num_samples_5deg+r),2)=R;
end
soundsc(snd,fs);

Simulating a sound source rotating in the elevation plane from -40 to 90 degrees using HRTFs.

Objective:

Simulate using HRTFs a source in a circular orbit with the head in the center of the circle. The source moves on the elevation plane (the azimuth angle is zero), with constant speed. Assume that Doppler effects are negligible.

Discussion:

This task is very similar to previous one, with one important exception: the HRTF set only defines elevation angles from -40 to 90 degrees. Therefore, it is impossible to simulate the source rotating in a complete circle with the head in the center. It is only possible to simulate a source rising (or falling) slowly from -40 to 90 degrees in fron of the observer. The code below is similar to the previous code. Here, the goal is to segment the samples into 14 slices 10 degrees each. Then the sound is convolved with HRTF of the appropriate elevation and azimuth zero. The results can be heard here. 
fs=44100;
T=10;
num_samples=T*fs;
num_samples_10deg=floor(num_samples/14);
snd=zeros(num_samples,2);
for i=0:13
    hrtf=readhrtf((i*10-40),0,'L');
    src=randn(num_samples_10deg,1);
    R=conv(hrtf(1,:),src);
    L=conv(hrtf(2,:),src);
    [r c]=size(R);
    snd((i*num_samples_10deg+1):(i*num_samples_10deg+r),1)=R;
    snd((i*num_samples_10deg+1):(i*num_samples_10deg+r),2)=L;
end
soundsc(snd,fs);

Simulating interaural time-difference using Doppler effect and HRTFs:

Simulation using Doppler Effect, intensity attenuation due to distance, and interaural cues.

Basic Setup


The setup is as follows: a source is moving at a constant speed along a trajectory, which is a traioght line. The stationary observer is at a certain constant distance away from the trajectory. Here, different techniques can be used separately or in combination with each other to produce the sound simulating the setup scenario.

Simulation using Doppler Effect, Interaural Time-Difference and Instensity Cues.

Objective:

Simulate the sound perceived by an observer that watches a source traveling at constant speed along a straight line trajectory, while observer is standing a constant distance away from the trajectory. Use Doppler Effect, Interaural-Time Difference and Intensity Cues.

Discussion:

The task required some optimization and reworking of the code from the previous project. First, it became necessary to re-design the function that produces the Doppler sound in such a way so that it receives start distance, end distance of the source relative to the normal from the observer (shown in the diagram above), in addition to the distance of the observer from the trajectory and the true speed of the source. This allows us to create a third function, called doppler_binaural, which accepts start distance, end distance, distance from trajectory, and the speed of the source, but computes two different Doppler sounds, combined in a stereo sound vector. Since the ears of the human being are usualy 17cm apart, doppler_binaural creates two Doppler sounds for each ear separately.

Function that computes the speed perceived by the observer:

This function accepts distance along the trajectory between the normal to the observer and the source position(D), the distance between the trajectory and the observer along the normal to the observer(dist) and the true speed of the source (true_speed). It outputs the speed of the source relative to (or perceived by) observer. 
function perceived_speed = observed_speed(D,dist,true_speed)
if D~=0
    theta = atan(dist/D);
else
    theta = pi/2;
end
perceived_speed = true_speed*cos(theta);

Function that computes the Doppler shift on the frequency:

This function calculates the shift in original frequency (freq) of the source traveling at a perceived speed (v_src). 
function perceived_freq = doppler(freq,v_src)
c = 341;
perceived_freq = freq.*(c./(c-v_src));

Function that generates the sound simulation of the Doppler Effect with intensity attentuation:

This function outputs a mono sound representing the tone received by the observer from the source traveling at a speed (dist) on the straight line trajectory that is a distance Dstart before and Dend after the normal to the observer, while observer is located at a distance (dist) away from trajectory. The function also attenuates the amplitude of the perceived tone by 1/r factor, r being the direct distance between the source and the observer. 
function snd = doppler_sound(Dstart,Dend,dist,speed)
f=440;
fs=8000;
Tstart=Dstart/speed;
Tend=Dend/speed;
t=-Tstart:(1/fs):Tend;
[row_t col_t] = size(t);
vel = zeros(1,col_t);
freq = zeros(1,col_t);
snd = zeros(1,col_t);
%observed velocity
for i=1:col_t
    if t(i)<=0
        vel(i)=observed_speed(abs(t(i)*speed),dist,speed);
    else
        vel(i)=-observed_speed(t(i)*speed,dist,speed);
    end
end
freq=doppler(f,vel);
for i=1:col_t
    if dist~=0
        snd(i)=sin(2*pi*freq(i)*t(i))/sqrt((dist)^2+(t(i)*speed)^2);
    else
        if t(i)==0
            snd(i)=sin(2*pi*freq(i)*t(i));
        else
            snd(i)=sin(2*pi*freq(i)*t(i))/abs(t(i)*speed);
        end
    end
end

Function that creates binaural simulation of the Doppler Effect (intensity attenuation due to distance included):

This function accepts the same set of inputs as doppler_sound function, and actually calls doppler_sound twice, both times for each ear. It alters the received parameters, however, so that the output of the each call to doppler_sound represent the tone received by one of the ears of the observer. The output of the function is a stereo sound, left and right channels of which contain the tone that would be received by the left and right ears of the observer, respectively. 
function snd=doppler_binaural(Dstart,Dend,dist,speed)
fs=8000;
half_ear=0.17/2;
snd_left=doppler_sound(Dstart-half_ear,Dend+half_ear,dist,speed);
snd_right=doppler_sound(Dstart+half_ear,Dend-half_ear,dist,speed);
snd_left=snd_left';
snd_right=snd_right';
snd=[snd_left snd_right];

Resulting Simulations:

Simulating a source traveling on a straight line using HRTFs

Objective:

Simulate a source traveling in a straight line, just like previous simulations, using HRTFs only.

Discussion:

The most difficult part of this assignment is conceptual understanding of how to partition the path of the source into chunks covered by 5-degree wide swaths of angles, in which each particular HRTF is defined.
The key to the solution is first calculating the angle vector, which holds the value if the 5-degree wide angle at each respective instance in time when the source travels. Next, find the range of indices in which the angle is same. That range will be indentical to the time vector indices, representing the time period in which the source is within a certain 5-degree angle. Generate the source for that time period and then convolve it with the HRTF, combining the result. It is also important not to forget that convolution produces more samples than the input, and that must be accounted for.
A signifcant problem however exists with this simulation. It does not:
  1. realistically attenuate the signal doe to distance
  2. ingores the Doppler shift
On other hand, this is a very versatile solution. It allows the user to set the starting and ending distance, distance from the trajectory and the speed. This allows a great control over the sound produced and can realistic simulate real-world (given a "real" world without attenuation and Doppler shift). 
function snd = straight_line_hrtf(Dstart,Dend,dist,speed)
fs=44100;
f0=440;
Tstart=Dstart/speed;
Tend=Dend/speed;
t=-Tstart:(1/fs):Tend;
[row_t col_t] = size(t);
snd=zeros(2,col_t);
angle=zeros(1,col_t);
elevation=0;
%get angle slices for the path
for i=1:col_t
    if t(i)<=0
        angle(i)=-floor(src_angle(abs(t(i)*speed),dist)/5)*5;
    else
        angle(i)=floor(src_angle(abs(t(i)*speed),dist)/5)*5;
    end
end

i=1;
j=1;
[r_angle col_angle]=size(angle);
while i <= col_angle
    %find indices for one angle slice
    while angle(1,i)==angle(1,j)
        j=j+1;
        if j>=col_angle
            break;
        end
    end
    %convolve source for the time duration indicated by the indices, with HRTF
    %at the azimuth angle specified:
    if angle(1,i)<0
        azimuth=-angle(1,i);
    else
        azimuth=angle(1,i);
    end
    hrtf=readhrtf(elevation,azimuth,'L');
    %generate the src
    src=randn(1,j-i);
    %src=sin(2*pi*f0.*t(i:j-1));
    %convolve
    R=conv(hrtf(1,:),src);
    L=conv(hrtf(2,:),src);
    [r c]=size(R);
    if t(i)<=0
        snd(1,i:(i+c-1))=R;
        snd(2,i:(i+c-1))=L;
    else
        snd(1,i:(i+c-1))=L;
        snd(2,i:(i+c-1))=R;
    end
    %start with next slice
    i=j;
end
R=snd(1,:);
L=snd(2,:);
snd=[R' L'];
The examples of the output:

Simulating a source traveling on a straight line using HRTFs, Doppler Effect and intensitu attenuation

Objective:

Simulate a source traveling in a straight line, just like previous simulations, using combination of HRTFs, Doppler Effect and intensity attenuation.

Discussion:

The approach is simple: generate two Doppler/Intensity effect sounds, for both ears. Next, alter these sounds with HRTFs, according to space position of source at the time when it emits this particular sound. Thus, HRTF function not only receives the distances and speed, from which it can calculate the position of the source, it also receives the sound produced by the source at that position and speed.
Accomplishing this teask required rewriting several functions. The key was making sure that functions work for mono only, in other words, that functions work only for one, specific ear. Since there is binaural simulation happening for both Doppler/Intensity effects and HRTFs, it is necessary to differentiate the two ears in space.
Since Doppler function already returned the single-row sound vector, it was used almost unchanged: 
function snd = doppler_sound2(Dstart,Dend,dist,speed)
%This function outputs a mono sound representing the tone received 
%by the observer from the source traveling at a speed (dist)  on 
%the straight line trajectory that is a distance Dstart before and 
%Dend after the normal to the observer, while observer is located at 
%a distance (dist) away from trajectory. The function also attenuates 
%the amplitude of the perceived tone by 1/r factor, r being the direct 
%distance between the source and the observer.
f=440;
fs=44100;
Tstart=Dstart/speed;
Tend=Dend/speed;
t=-Tstart:(1/fs):Tend;
[row_t col_t] = size(t);
vel = zeros(1,col_t);
freq = zeros(1,col_t);
snd = zeros(1,col_t);

%observed velocity
for i=1:col_t
    if t(i)<=0
        vel(i)=observed_speed(abs(t(i)*speed),dist,speed);
    else
        vel(i)=-observed_speed(t(i)*speed,dist,speed);
    end
end

freq=doppler(f,vel);

for i=1:col_t
    if dist~=0
        snd(i)=sin(2*pi*freq(i)*t(i))/sqrt((dist)^2+(t(i)*speed)^2);
    else
        if t(i)==0
            snd(i)=sin(2*pi*freq(i)*t(i));
        else
            snd(i)=sin(2*pi*freq(i)*t(i))/abs(t(i)*speed);
        end
    end
end
The HRTF function had to be reworked. Several issues had to be dealt with:
  1. The function has to receive the source sound vector that it mjust alter by HRTFs, it cannot generate sound on its own anymore. It receives the sound that the source produces at the given position and speed.
  2. The function must return only one channel of HRTF, depending on whether left ear or right ear HRTF response is needed.
function snd = straight_line_hrtf2(Dstart,Dend,dist,speed,full_src,ear)
%snd = straight_line_hrtf2(Dstart,Dend,dist,speed,full_src,ear): generates
%the HRTF for the ear specified, assuming the the source emits full_src (mono sound) while traveling in 
%straight line with constant speed, with observer (the ear) being a
%constant distance (dist) away. "ear" must be R for right and L for
%left.
fs=44100;
f0=440;
%##############
Tstart=Dstart/speed;
Tend=Dend/speed;
t=-Tstart:(1/fs):Tend;
[row_t col_t] = size(t);
snd=zeros(1,col_t);
angle=zeros(1,col_t);
elevation=0;
%get angle slices for the path
for i=1:col_t
    if t(i)<=0
        angle(i)=-floor(src_angle(abs(t(i)*speed),dist)/5)*5;
    else
        angle(i)=floor(src_angle(abs(t(i)*speed),dist)/5)*5;
    end
end

i=1;
j=1;
[r_angle col_angle]=size(angle);
while i <= col_angle
    %find indices for one angle slice
    while angle(1,i)==angle(1,j)
        j=j+1;
        if j>=col_angle
            break;
        end
    end
    %convolve source for the time duration indicated by the indices, with HRTF
    %at the azimuth angle specified:
    if angle(1,i)<0
        azimuth=-angle(1,i);
    else
        azimuth=angle(1,i);
    end
    hrtf=readhrtf(elevation,azimuth,'L');
    %generate the src
    src=full_src(1,i:j-1);
    %convolve
    R=conv(hrtf(1,:),src);
    L=conv(hrtf(2,:),src);
    [r c]=size(R);
    if ear=='R'
        if t(i)<=0
            snd(1,i:(i+c-1))=R;
        else
            snd(1,i:(i+c-1))=L;
        end
    elseif ear=='L'
        if t(i)<=0
            snd(1,i:(i+c-1))=L;
        else
            snd(1,i:(i+c-1))=R;
        end 
    else
        error('Ear must be either an R for right or L for left');
    end
    %start with next slice
    i=j;
end
The final function is the one that combines both Doppler/Intensity and HRTF functions, for left and right ears separately, and then combines these two into one two-channel audio file: 
function snd = complete_doppler3(Dstart,Dend,dist,speed)
fs=44100;
f0=440;
%##############
% SOURCE IS ASSUMED TO TRAVEL FROM TO LEFT TO RIGHT
%##############
%get left and right ear Doppler/Instensity effects
half_ear=0.17/2;
snd_left=doppler_sound2(Dstart-half_ear,Dend+half_ear,dist,speed);
snd_right=doppler_sound2(Dstart+half_ear,Dend-half_ear,dist,speed);
%snd=[snd_left' snd_right'];
L=straight_line_hrtf2(Dstart-half_ear,Dend+half_ear,dist,speed,snd_left,'L');
R=straight_line_hrtf2(Dstart-half_ear,Dend+half_ear,dist,speed,snd_left,'R');
snd=[L' R'];
soundsc(snd,fs);

Examples:

Conclusion

This project is an exciting experiment with real mathematical data used in audio-related research. The 3D sound localization and HRTFs are an extremely interesting area and experience in this is definetly something to be proud of. The project seemed to be much more daunting at the beginning, and one has to admit that it is a great feeling to know how difficult it seemed at first and experience actually completing the project successfully.