sfstoolbox / sfs-matlab

SFS Toolbox for Matlab/Octave
https://sfs-matlab.readthedocs.io
MIT License
97 stars 39 forks source link

Correct amplitude wfs time domain #112

Closed VeraE closed 7 years ago

VeraE commented 7 years ago

I suggest the following changes to achieve a correct amplitude of a point source in WFS when in the time domain:

  1. Correction of a mistake in the default 2.5D WFS driving function driving_function_imp_wfs_ps().

  2. Adding the missing amplitude terms to the FIR prefilter in wfs_fir_prefilter(), both for 2.5D and 3D.

  3. Change amplitude and delay of dummy_irs() to 1/(4*pi) in 1 m distance.

These changes can be found in the three commits. What do you think of this? Do these changes conflict with some things I missed?

hagenw commented 7 years ago

Thanks for the corrections, here are my thoughts on them:

  1. This was indeed a mistake before and can be replaced by your code
  2. + 3. This looks nice, especially as it is now conform with the cited equations.

I have one question regarding 2.+3.: I checked it by calculating the frequency response of WFS, after http://matlab.sfstoolbox.org/en/2.2.1/binaural-simulations/#frequency-response-of-your-spatial-audio-system and run the usual test_driving_functions(1). The results is still the same. Do you have an example or an idea where your changes might be critical in the time domain?

VeraE commented 7 years ago

Both your examples contain normalisations of the sound field at some point. You can try this example which demonstrates that the level of a real point source is achieved for the virtual source as well:

clear; clc

SFS_start
SOFAstart

%% conf struct
conf = struct;

%basics
conf.fs = 44100;
conf.c = 343;
conf.showprogress = false;
conf.debug = 0;

%secondary sources
conf.secondary_sources.geometry = 'line';
conf.secondary_sources.size = 20;
conf.secondary_sources.number = 160;
conf.secondary_sources.center = [0 0 0];
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 = 128;
conf.wfs.hpreflow = 50;
conf.wfs.hprefhigh = conf.c/2/x0(1,7);
conf.xref = [0 -1 0];
conf.wfs.t0 = 'source';
conf.usetapwin = false;
conf.tapwinlen = .3;

%interpolation, delaying, weighting
conf.ir.useinterpolation = true;
conf.ir.interpolationmethod = 'freqdomain';
conf.ir.hrirpredelay = 0;
conf.delayline.resampling = 'matlab';
conf.delayline.resamplingfactor = 100;
conf.delayline.filter = 'integer';

hrtf = dummy_irs(conf);
conf.ir.usehcomp = false;
conf.N = 8000;

%% Listener and source position
X = conf.xref; %listener position
xs = [0 1 0]; %source position

%% Calculate IR
[ir,x0] = ir_wfs(X,0,xs,'ps',hrtf,conf);

%% Plot
IR = fft(ir);
N = length(ir);
f = (0:N-1)/N*conf.fs;

figure
    semilogx(f,db(abs(IR)))
    grid
    set(gca,'Xlim',[10 conf.fs/2])
    xlabel('frequency / Hz'), ylabel('dB')

%% Real point source at 1 kHz
r = vector_norm(X-xs,2);
Pps = exp(-1i*2*pi*1000/conf.c*r)/(4*pi*r);
disp(['Level real point source at listener position in dB: ' num2str(db(abs(Pps)))])
hagenw commented 7 years ago

That's true, as usual I forgot the normalization. Your example worked and I will merge your changes. Thanks again.