Correct Green's function amplitude

[VeraE]

This change corrects the amplitude of an SSR-BRS auralisation for the factor 1/(4*pi) of the Green's function. In order to work for calculation of WFS impulse responses in the time domain as well, this factor had to be removed from dummy_irs() (where it had been added in #112).

[hagenw]

Hi @VeraE,

Thanks for trying to get the amplitude correct, I always ignored that ;) I have three questions: 1) In get_ir(), ir = ir/(4*pi); is written at two positions at the moment. I guess we can reduce it to one time, if we place it after line 210? 2) I use the get_ir() function not only for sound field simulations, but also for just getting a HRTF out of a SOFA file. Before, if I specified a measured direction and distance, it returned indeed the measured one without modifications. Now, there is an mismatch in amplitude. But I also understand your point of getting the amplitude correct. Do you have further thoughts on this? 3) If I see it correctly, the amplitude should be correct inside greens_function_imp(). Does that mean, if I calculate the sound field amplitude at one position, one time with ir_wfs() and one time with sound_field_imp_wfs(), that the resulting amplitude is the same?

[hagenw]

For point 3) I mean to use ir_wfs() together with dummy_irs() and compare this to sound_field_imp_wfs().

Regarding point 2), it might be a solution to move ir = ir/(4*pi); to ir_generic() directly after line 87.

[VeraE]

Regarding point 1) and 2): I moved the factor 1/(4*pi) to ir_generic(), that's fine with me. So get_ir() can still be used for loading from SOFA files without changing the amplitude.

Regarding point 3): I tried this for one example and compared the amplitude of the first pulse of the WFS impulse response at one listener position for calculation with ir_wfs() and with sound_field_imp(), which are the same. It is very tedious though to test this further as one has to step through sample-by-sample iterations of sound_field_imp(). Are you in favour of testing this more thoroughly?

[hagenw]

I would say that it should be fine to integrate your pull request as it is at the moment. If other users are really, really interested in getting the amplitude correct, they will test and complain ;) But could you shortly post the code you used to test it?

[VeraE]

This is my code for testing:

SFS_start

%% conf struct for SFS Toolbox
conf = struct;

%basics
conf.fs = 44100;
conf.c = 343;
conf.N = 8000;
conf.showprogress = false;
conf.debug = 0;
conf.tmpdir = '';
conf.resolution = 300;

%secondary sources
conf.secondary_sources.geometry = 'line';
conf.secondary_sources.size = 4;
conf.secondary_sources.number = 41;
conf.secondary_sources.center = [0 0 0];
conf.secondary_sources.x0 = secondary_source_positions(conf);

%wfs
conf.dimension = '2.5D';
conf.driving_functions = 'default';
conf.wfs.usehpre = true;
conf.wfs.hpretype = 'FIR';
conf.wfs.hpreFIRorder = 512;
conf.wfs.hpreflow = 50;
conf.wfs.hprefhigh = conf.fs/2;
conf.xref = [0 -1 0];
conf.wfs.t0 = 'source';
conf.usetapwin = true;
conf.tapwinlen = .3;

%interpolation, delaying, weighting, interpolation
conf.ir.useinterpolation = true;
conf.ir.interpolationmethod = 'freqdomain';
conf.ir.hrirpredelay = round(1/conf.c*conf.fs);
conf.delayline.resampling = 'matlab';
conf.delayline.resamplingfactor = 100;
conf.delayline.filter = 'integer';

%bandpass for impulse
conf.usebandpass = true;
conf.bandpassflow = 10;
conf.bandpassfhigh = 20000;

%use dirac as hrir for ir_wfs() and no headphone compensation
hrtf = dummy_irs(conf);
conf.ir.usehcomp = false;

%plotting
conf.plot.useplot = true;
conf.plot.usenormalisation = false;
conf.plot.normalisation = 'max'; 
conf.plot.cmd = '';
conf.plot.usedb = false;
conf.plot.mode = 'monitor';
conf.plot.size_unit = 'px';
conf.plot.size = [540 404];
conf.plot.caxis = [-.1 .1];
conf.plot.colormap = 'default';
conf.plot.loudspeakers = true;
conf.plot.realloudspeakers = false;
conf.plot.lssize = 0.16;
conf.plot.usefile = false;
conf.plot.file = '';

%% Calculate WFS impulse response with ir_wfs()
X = [0 -1 0]; %listener position
phi = 0; %listener direction in rad
xs = [0 2 0]; %source position

ir = ir_wfs(X,phi,xs,'ps',hrtf,conf);

%Plot
figure
    plot(0:conf.N-1,ir)
    grid
    xlabel('samples')

%% Calculate sound field with sound_field_imp()
%listener positions
X = [-2 2];
Y = [0 -1.5];
Z = 0;
%source position
xs = [0 2 0];

%time instant in samples
t = 385;

sound_field_imp_wfs(X,Y,Z,xs,'ps',t,conf)

The amplitude is not exactly the same but quite close (0.032 vs. 0.027 for the first pulse at listener position (0,-1) for ir_wfs() and sound_field_imp(), respectively).

[hagenw]

There is the function time_response_wfs() that should be suitable for the comparison. I also did a small test script:

conf = SFS_config;
X = [0 0 0];
phi = 0;
xs = [2.5 0 0];
src = 'ps';
hrtf = dummy_irs(conf);
% Using an impulse response
ir = ir_wfs(X,phi,xs,src,hrtf,conf);
% Using sound_field_imp() for every time step
[s,t] = time_response_wfs(X,xs,src,conf);
figure;
plot(0:700,ir(1:701,:),'-g');
hold on
plot(t,s,'-b');
grid;
xlabel('samples');
hold off;

It looks like there is still a slight mismatch in amplitude as well as in delay. At least the amplitude is closer matched then before. For comparison I run the script for version fc70479 (before you changed the amplitude the first time):

Interestingly, the delay changed as well and I don't know why yet.

[hagenw]

I thought a little bit about the time delay. In my example the point source is 2.5m away from the listener. The first peak of the simulation done with ir_wfs() in the first figure (green line) is at 385 samples, the driving function adds an additional delay of 60 samples (It can be returned by driving_function_imp_wfs() if one wants to know it), resulting in:

(385-60)/fs*c = 2.53

which is near the desired 2.5 m.

Now, I will try to see why this is not the case for time_response_wfs() (blue line in the first figure).

[hagenw]

I moved the discussion of the delay to #115

[VeraE]

Ah, I did not know about time_response_wfs().

I tested your example a bit. I think in the time domain the peaks just have different sizes because time_response_wfs() uses fractional delay and ir_wfs() with the configuration of SFS_config() does not. If you plot the corresponding frequency responses with

conf = SFS_config;
X = [0 0 0];
phi = 0;
xs = [2.5 0 0];
src = 'ps';
hrtf = dummy_irs(conf);
% Using an impulse response
ir = ir_wfs(X,phi,xs,src,hrtf,conf);
% Using sound_field_imp() for every time step
[s,t] = time_response_wfs(X,xs,src,conf);
figure;
plot(0:700,ir(1:701,:),'-g');
hold on
plot(t,s,'-b');
grid;
xlabel('samples');
hold off;

%frequency vectors
f_ir = (0:length(ir)-1)/length(ir)*conf.fs;
f_s  = (0:length(s)-1)/length(s)*conf.fs;
%plot
figure
    semilogx(f_ir,db(abs(fft(ir(:,1))))), hold on
    semilogx(f_s,db(abs(fft(s))))
    grid
    xlabel('frequency / Hz'), ylabel('amplitude / dB')
    set(gca,'XLim',[10 conf.fs/2])
    set(gca,'YLim',[-60 -20])
    legend('ir wfs','time response wfs','location','SW')

You can see that the level of the two examples is the same. I don't know where the difference at very high frequencies comes from, though.

[hagenw]

Your suspicion of the fractional delay seems to be correct. I'm still not completely sure why the amplitude is that much reduced as using simply resampling as a fractional delay filter does not give such a deviation. But I did another way to test this. I added time_response_point_source() and did a comparison with ir_point_source():

conf = SFS_config;
X = [0 0 0];
phi = 0;
xs = [2.5 0 0];
hrtf = dummy_irs(conf);
ir = ir_point_source(X,phi,xs,hrtf,conf);
[s,t] = time_response_point_source(X,xs,conf);
figure;
plot(0:700,ir(1:701,:),'-g');
hold on
plot(t,s,'-b');
grid;
xlabel('samples');
hold off;

The point source is placed 2.5 m away, so its expected amplitude is 1/(4*pi*2.5) = 0.032, which is correctly represented by ir_point_source(). As time_response_point_source() implements exactly the same equation and I can't see any error, I would also suspect that the amplitude deviation is solely given by fractional delay. @fietew: do you agree with this?

[fietew]

@hagenw: Well yes and no. First of all you'll always lose amplitude/power, if fractional delay filters are applied since all of them have lowpass characteristic (depending on the fractional part of the delay). See:

conf = SFS_config;
conf.delayline.filter = 'lagrange';

tfrac = 0.5;

d = [zeros(1,1023), 1, zeros(1,1024)].';
dspline = interp1(1:length(d), d, (1:length(d))-(tfrac+100), 'spline');
conf.delayline.filterorder = 3;
dfrac1 = delayline(d, tfrac+200, 1, conf);
conf.delayline.filterorder = 20;
dfrac2 = delayline(d, tfrac+300, 1, conf);
conf.delayline.filterorder = 100;
dfrac3 = delayline(d, tfrac+400, 1, conf);

figure;
plot(d), hold on
plot(dspline);
plot(dfrac1);
plot(dfrac2);
plot(dfrac3);
hold off

% overall power
sum(d.^2)
sum(dspline.^2)
sum(dfrac1.^2)
sum(dfrac2.^2)
sum(dfrac3.^2)

regarding your "error", I first tried to adapt the distance such that it results in an integer delay, see:

conf = SFS_config;
X = [0 0 0];
phi = 0;
r = round(2.5/343*44100) * 343/44100;
xs = [r 0 0];
hrtf = dummy_irs(conf);
ir = ir_point_source(X,phi,xs,hrtf,conf);
[s,t] = time_response_point_source(X,xs,conf);
figure;
plot(0:700,ir(1:701,:),'-g');
hold on
plot(t,s,'-b');
grid;
xlabel('samples');
hold off;

However, time_response_point_source() does still not yield a dirac impulse. After some headaches I finally discoverd that time_response_point_source() applies a bandpass filtering (see https://github.com/sfstoolbox/sfs-matlab/blob/master/SFS_time_domain/sound_field_imp.m#L155 ). After setting conf.usebandpass = false, the above script yields two dirac impulses.

[hagenw]

Ah, cool. You are right, now it looks more like I had imagined. In the figure below you find again the WFS example from above, but setting conf.usebandpass = false. The fractional delay still has an influence, but the difference in amplitude is now much smaller.

I will first merge this into the master and after that I will create a pull request/issue to discuss where the bandpass is needed at all and if we can maybe change the conf.usebandpass to another default.

[hagenw]

The changes of this branch that we agreed to integrate are now part of the master as commit aa2dfce and d3b9b60, the other changes of this pull request will be ignored.