Menu

#14 Low frequency roll-off

2.0
closed
nobody
None
2020-05-05
2019-02-14
Isaac
No

Hello,

First of all, apologies for posting so many tickets! I really appreciate the work you put on this project.

A problem that I am finding is that the HRTFs that I obtain from the BEM calculation have a consistent roll-off towards the lower frequencies of approximately 20 dB/decade, for all directions. For instance, the sound pressure level at 100 Hz is 20 dB lower than at 1 kHz, for all azimuths and elevations. This seems to happen for any evaluation grid or 3D model.

For example, following the tutorial of this wiki, I created an Icosphere with subdivisions=6 and size=100 in Blender, at coordinates (0,0,0), and with "Left ear" and "Right ear" materials assigned to the leftmost and rightmost triangle of the mesh (along the Y axis). I exported it using the Mesh2HRTF plugin (100Hz-20kHz, ARI evaluation grid). Find attached some pictures with the result of plotting the HRTF for the left "ear". You can see how the magnitude steadily increases as you go towards the high frequencies.

Maybe it is something wrong in my setup? Is it necessary to normalize the data across frequencies?

Thank you very much!

Isaac

PS: Find below the Matlab code that I used to generate the figures. Basically, I take the real and imaginary part of the data, add them together, calculate the absolute value, and plot it in dB.

%% Extract SOFA data
sofaFile = SOFAload('EvaluationGrid_GeneralTF.sofa');
Real                = sofaFile.Data.Real;
Imag                = sofaFile.Data.Imag;
Value               = Real + Imag*1i;
Mag                 = abs(Value);
Phase               = angle(Value); 
Freqs               = sofaFile.N;
SourcePosition      = sofaFile.SourcePosition;
X                   = SourcePosition(:,1);
Y                   = SourcePosition(:,2);
Z                   = SourcePosition(:,3);
[Azim,Elev,Rad]     = cart2sph(X,Y,Z);
Azim                = -Azim; % azimuth sign is flipped

%% Plot horizontal plane HRTFs as image (left ear)
pos = find(abs(Elev)<0.01); % find HRTFs with 0 elevation
azim = Azim(pos);
[azim,ind] = sort(azim); % sort ascending
azim = azim*180/pi; % rad to deg
if max(azim) > 180
    azim = azim-180; % bring it to the -180/180 range
end
pos = pos(ind);
mag = Mag(pos,:,:);

figure
s = surf(azim,Freqs,db(squeeze(mag(:,1,:))')); % db(x) == 20*log10(x)
view(2) % vertical view
s.EdgeColor = 'none'; % remove surface lines
set(gca,'yscale','log') % log scale for frequency axis
xlim([min(azim),max(azim)]),ylim([min(Freqs),max(Freqs)])
colorbar
set(gca,'YDir','normal')
title(['HRTF Spectrum (elevation = 0deg): left'])
ylabel('Frequency [Hz]'), xlabel('Azimuth [deg]')

%% Plot frontal HRTF (left ear)
pos = find(abs(Azim)<0.01&abs(Elev)<0.01); % find HRTFs with 0 azimuth
mag = squeeze(Mag(pos,:,:));
phase = Phase(pos,:,:);
uw_phase = unwrap(phase,[],3);

figure
subplot(2,1,1)
semilogx(Freqs,db(mag))
xlim([Freqs(1),Freqs(end)])
legend('Left','Right')
xlabel('f [Hz]'), ylabel('mag [dB]')
title('Frontal left HRTF - Magnitude')
grid on
subplot(2,1,2)
semilogx(Freqs,squeeze(uw_phase))
xlim([Freqs(1),Freqs(end)])
legend('Left','Right')
xlabel('f [Hz]'), ylabel('phase [rad]')
title('Frontal left HRTF - Phase')   
grid on
3 Attachments

Discussion

  • Piotr Majdak

    Piotr Majdak - 2019-02-21

    Dear Isaac,
    thanks for report and interest in Mesh2HRTF. Unfortunately, I currently have no resource to provide extensive support.
    As for your observations about increasing pressure with frequency: you probably have to divide by the frequency. (the velocity condition seems to be mixed with the pressure condition). Hope that helps so far,
    Piotr

     
  • Isaac

    Isaac - 2019-02-26

    Hi Piotr,
    Thanks for the reply and for the honesty!
    I will try your suggestion.
    Isaac

     
  • Fabian Brinkmann

    Hi Isaac,

    if you use reciprocity, you have to set reference=true in Output2HRTF. Otherwise the frequency response of the volume velocity source is not compensated from the results.

    Best, Fabian

     
  • Fabian Brinkmann

    • status: open --> closed
     

Log in to post a comment.

MongoDB Logo MongoDB