Keep wfs delay

[trettberg]

The current implementation of time-domain WFS minimizes the delay: At least one loudspeaker starts playing at t = 0. This is very convenient behaviour for a single virtual source.

For multiple virtual sources this seems very unconvenient: You have to figure out how much delay was removed or added. So you basically calculate all delays again.

Proposed Optional Behaviour

I'd like to propose a flag conf.wsf.removedelay:

• true (default) results in the previous behaviour, suitable for a single virtual source
• false does the following:
• A common delay offset t0 is introduced for all virtual sources.
• t0 ensures causality for all plane waves and focused sources. Restrictions discussed below!
• t0 := diam(x0) / c corresponds to the maximum travel time through the array, where
  • diam(x0) denotes the diameter of the SSD, i.e. the maximum distance between two secondary sources. (agrees with conf.secondary_sources.size for 'line', 'circle' and 'sphere').

    Example

Three virtual sources (ps,pw and fs), shortly after t0:

conf = SFS_config();
conf.secondary_sources.geometry = 'circle';
x0 = secondary_source_positions(conf);
conf.wfs.usehpre = 0;
conf.wfs.removedelay = 0;

xs_ps = [1.6,0,0];
xs_pw = [1,-2,0];
xs_fs = [0,-1,0,0,1,0];

[d_ps,delay_ps] = driving_function_imp_wfs(x0,xs_ps,'ps',conf);
[d_pw,delay_pw] = driving_function_imp_wfs(x0,xs_pw,'pw',conf);
[d_fs,delay_fs] = driving_function_imp_wfs(x0,xs_fs,'fs',conf);
[~,sss_ps] = secondary_source_selection(x0,xs_ps,'ps');
[~,sss_pw] = secondary_source_selection(x0,xs_pw,'pw');
[~,sss_fs] = secondary_source_selection(x0,xs_fs,'fs');

d_sum = zeros(conf.N,conf.secondary_sources.number);
d_sum(:,sss_ps) = d_sum(:,sss_ps) + d_ps(:,sss_ps);
d_sum(:,sss_pw) = d_sum(:,sss_pw) + d_pw(:,sss_pw);
d_sum(:,sss_fs) = d_sum(:,sss_fs) + d_fs(:,sss_fs);
sss_sum = sss_ps | sss_pw | sss_fs;

conf.wfs.usehpre = 1;
d_sum = wfs_preequalization(d_sum,conf);

%% Plot
t0 = secondary_source_maximum_distance(x0)/conf.c; % reference time
t0 = round(t0*conf.fs) + 64 + 1; % in samples + prefilter delay 
t = t0 + 20;            % some samples more to see the point source

conf.plot.usedb = 1;
conf.plot.usenormalisation = 0;
sound_field_imp([-2,2],[-2,2],0,x0(sss_sum,:),'ps',d_sum(:,sss_sum),t,conf);
hold;
scatter(0,0,'kx');               % origin 
scatter(xs_ps(1),xs_ps(2),'ko'); % point source   
scatter(xs_fs(1),xs_fs(2),'ko'); % focused source
hold off;

Restrictions

Focused sources may only be placed within the ball B with radius diam(x0)/2, that encloses the SSD.

• (Note that for any closed SSD, the entire listening area lies within B, so that's fine.)
• Rationale:
  • In an unbounded listening area, a focused source could be arbitrarily far away from the SSD, resulting in arbitrarily large negative delay. The admissable region for focussed sources has to be restricted in some way or another.
  • Personally, I like this simple rule. (It also suggests a work-around for itself.)
  • The alternative would be a somewhat obscure conf.wfs.maximum_focused_source_distance that is only relevant for this particular case.

    Changes

• additional flag conf.wsf.removedelay
• internal changes in driving_function_imp_wfs()
• additional function secondary_sources_maximum_distance which computes the diameter and the center of B.

Thoughts?

[fietew]

I discussed a similar issue with Hagen related to the additional delay added by the fractional delay filters. Imho, it would be easier to return the pre-delay which has been added to the driving function by a particular processing step, e.g. fractional delay, wfs prefilter, focused source, plane wave, etc.. Furthermore, I would prefer that t=0 corresponds to the t=0 of the virtual source, i.e. where the sound "starts" to radiate from a point/focused source or when the plane wave passed its point of zero phase.

[trettberg]

Returning the introduced pre-delay would work. Maybe it's a better strategy in the bigger picture, I don't know. But I would have to correct the delay after each call to driving_function_imp_wfs. I would prefer thatdriving_function_imp_wfs simply "does it right in the first place." (Context: The use case I have in mind is simulation of databased SFS, i.e. many virtual sources.)

~~I'm not sure if understand your second point correctly. t0 is the shifted t = 0 for all virtual sources. Or put differently, each virtual source has the time signal delta(t - t0). The plot merely depicts a time instant short after t0, because I wanted the point source to be visible.~~ If no prefilter is used then t0 happens to be the overall delay. In this case, the virtual sources have the time signal delta(t - t0).

[trettberg]

Upon re-reading: My motivation was not very clear and parts of the code example were misleading. Also my response to Fiete was partly wrong.

An attempt to clarify:

What it's about: a bit more superposition

If I have a sum of simple virtual sources (e.g. an image source model), I just want to sum up the driving functions. (For a sum of plane waves, I'd still have to take care for the arrival times. But I think it's step in the right direction.)

What it's not about: reference time

I make no assumptions about the total delay at the end of processing. The user is still responsible to keep track of that.

The corresponding part in the code example should have read like this:

%% Plot
t0 = secondary_source_maximum_distance(x0)/conf.c; % max. travel time
total_delay = round(t0*conf.fs) + conf.wfs.hpreFIRorder/2  % total delay in samples 
t = total_delay + 20;            % some samples more to see the point source

[fietew]

I am totally aware, that its about the correct superposition of driving functions, but let me phrase the problem(s) the other way around:

1. "there are some issues with the time-correct superposition of driving functions"
2. "there are some issues with the reference time"

Solving 1. does not necessarily solve 2. and might complicate solving 2. afterwards. I think that, if we solve 2., 1. should be solved automatically or 1. should be solvable with moderate effort. I would therefore prefer to find a solution for the bigger picture .

[trettberg]

I agree that if driving_function_imp_wfs returns the pre-delay as you proposed, the superposition would be simple enough: just use delayline again. Would definitely work for me. Still I do not consider this as very elegant for the particular case: If (for whatever reason) fractional delays are used with a suboptimal method, you'd have twice the artifacts. (Though I think this probably doesn't matter much in practice?)

I also agree that it would be cool to have the reference time calculated automatically.

[hagenw]

I started testing the current code:

1.) The secondary_source_maximum_distance() function is a nice idea and ready for integration from my perspective. The only think I would like to mention here, is that you should have in mind that for a circle or sphere it will not return the diameter of the circle on which the secondary sources are located, but the maximum distance between the secondary sources, which will be a little bit smaller in most cases. Also the center will then be a little of the actual center of the circle. But in my opinion this is ok, as we are interested in the actual sampled secondary sources.

2.) driving_function_imp_wfs(): I get the point of not applying the automatic delay adjustment and think the implementation via conf.wfs.removedelay is a good idea. For the current implementation of the else case I'm not 100% convinced. The calculation of the diameter is probably not doing what you intended to achieve. As the x0 you feed into the function are the already selected secondary sources, they are no longer a full circle, sphere, or box, but smaller parts of those. This also means that the checking for focused sources fails in most cases. For example

>> conf = SFS_config;
>> conf.wfs.removedelay = false;
>> sound_field_imp_wfs([-2 2],[-2 2],0,[0.5 0.5 0 -1 -1 0],'fs',200,conf)
error: DRIVING_FUNCTION_IMP_WFS: Using 'config.wfs.removedelay == 0', focused source positions are restricted to the ball spanned by the array diameter.
error: called from
    driving_function_imp_wfs at line 132 column 13
    sound_field_imp_wfs at line 98 column 3

As the actual calculated diameter is not 3m, but only 1.54m for the selected secondary sources.

Before fixing these issue I would like to ask, if this correction is needed at all? I think if the users set conf.wfs.removedelay = false; they have to play around with t anyhow to find the correct time step.

[hagenw]

I tested a little bit further and I have to admit, that the code under the else statement in driving_function_imp_wfs() helps if you want to create a plot as you have shown above.

[trettberg]

The calculation of the diameter is probably not doing what you intended to achieve. As the x0 you feed into the function are the already selected secondary sources, they are no longer a full circle, sphere, or box, but smaller parts of those.

Sad but true: I completely missed that. Somehow I got into the habit of getting the driving function first and doing the selection afterwards.

A solution would be to pass conf to secondary_sources_maximum_distance instead of x0.

• For the preset geometries, the diameter and center can be directly taken from the size and center entries. That would also solve the problem for spheres and circles that you pointed out.
• For custom geometries, if x0 is present in conf, use that. Otherwise revert to size and calculate some upper bound.

Before fixing these issue I would like to ask, if this correction is needed at all? I think if the users set conf.wfs.removedelay = false; they have to play around with t anyhow to find the correct time step.

Do you mean the restriction of possible focused source locations? I think it's nice to ensure that no driving function is cut off in the beginning, which would otherwise happen with open arrays.

[hagenw]

I changed secondary_source_maximum_distance() to use now conf, which is an elegant solution as this avoids the problems we have with the secondary source selection in both cases, predefined geometries and custom geometries.

In your example plot, you just have to change one line:

t0 = secondary_source_maximum_distance(conf)/conf.c; % max. travel time

With "if this correction is needed at all" I meant the whole statement after else. But in the meantime I figured out, that it seems to be needed as otherwise it is not so easy to get the result of your figure from above.

[hagenw]

Maybe it would be better to rename secondary_source_maximum_distance() to secondary_source_diameter()?

[trettberg]

Agreed, that's a better name.

Also maybe minimizedelay would be more precise than removedelay?

One question about array geometries: Is there a consistent definition of size? (If not, that's fine and probably more convenient for adding other predefined setups.)

[hagenw]

The definition of size is a bit tricky. At the moment we have linear, circular, box, rounded-box, and spherical as predefined secondary source distributions, where size is used.

1.) Simple cases: a) linear: size == length of array b) circular: size == diameter of array c) spherical: size == diameter of array 2.) More complicated cases: a) box: size == length of one linear array element (which is non-identical to the distance between two parallel linear array elements, this one is size + 2*dx0 if I remember correctly) b) rounded-box: size == remove the round edges and use the definition from box

[hagenw]

I prefer conf.wfs.removedelay as the explanation also starts with "Remove the leading delay in WFS-time domain ...". Something like adjustt0 should also be fine, but I'm not sure if everyone knows what t0 should be.

Regarding the renaming of secondary_source_maximum_distance(), what do you prefer: 1) secondary_source_diameter() 2) secondary_source_size()

I think, I would go for 2) as the value is more or less identical to conf.secondary_sources.size.

[fietew]

rounded-box: size == remove the round edges and use the definition from box

I am not sure, what you mean by "remove" in that context, but size in this case means the size of a square bounding box around the rounded-box.

[hagenw]

@trettberg you are right, there is a small difference between box and rounded-box regarding its diameter, which I was not aware of before. That is the reason why I renamed secondary_source_maximum_distance() to secondary_source_diamter() and not to *_size() as it returns different values in the case of box or rounded-box.

In my opinion the code is now ready for integration into master. I also see no problem to change it afterwards to be more in line with the feature @fietew requested. For this I created a follow up at #61.

If you have no objections I would merge this pull request and close it.

[trettberg]

Should be ready to merge.

[hagenw]

One last question. I added a test function test_driving_functions_imp_with_delay in order to check the new delay behavior for conf.wfs.removedelay = false; and used your proposed calculation of t with the only difference that I added 30 extra samples instead of 20. If you run the test function with test_driving_functions_imp_with_delay(1) a lot of figures will pop up and you will realize that there is a larger difference between focused sources and point source. Is this ok? I guess yes, as you would argue that the impulse has traveled roughly the same time from the position of the virtual source? Here, an example:

[trettberg]

I think this is okay. To observe the timing, I think it helps to place the point source unreasonably close to the array. Otherwise, the focussed source is already "all over the place", before the point source even starts. (For the circular array, the point source distance in the test scenario is 1m, while my initial code example had 0.1m).