add monochromatic version of LWFS-SBL

[fietew]

• circular basis expansion of monochromatic plane wave and point source
• extension of driving_function_mono_wfs_vss to handle plane wave and focused source as virtual secondary sources (similar to time domain implementation)
• add driving_function_mono_wfs_pwd to handle plane wave decompositions
• add pwd_mono_circexp to convert circular basis expansion into plae wave decomposition

[hagenw]

So far we have always used 1i as the imaginary number in the toolbox. In your functions you have used 1j. If you don't mind, you could switch to 1i as well.

[hagenw]

If have tested it so far by running:

conf = SFS_config;
conf.plot.normalisation = 'center';
sound_field_mono_localwfs_sbl([-2 2],[-2 2],0,[0 2.5 0],'ps',1000,conf)
sound_field_mono_localwfs_sbl([-2 2],[-2 2],0,[0 1 0],'pw',1000,conf)

which resulted in those figures, that look fine to me:

Would it be worth to add a test_localwfs_sbl() function parallel to test_localwfs_vss()? In addition, we have test_exp_imp(), but nothing for the monaural case.

[fietew]

I fixed the point source. Could you please test again?

[hagenw]

The differences are very subtle, but now the figure is exactly symmetric, which was not the case before.

[fietew]

So far we have always used 1i as the imaginary number in the toolbox. In your functions you have used 1j. If you don't mind, you could switch to 1i as well.

fixed in 7dc3abf

[hagenw]

Both testing functions are running, but of course I cannot guarantee that all of the resulting figures look like they should look ;) For test_localwfs_sbl(0) I get a lot of the following warnings, which doesn't happen for test_localwfs_vss(0):

Warning: SECONDARY_SOURCE_SELECTION: 0 secondary sources were selected. 
> In secondary_source_selection (line 176)
  In driving_function_mono_wfs_vss>@(X0,XS)secondary_source_selection(X0,XS(1:3),'pw') (line 98)
  In driving_function_mono_wfs_vss (line 108)
  In driving_function_mono_wfs_pwd (line 78)
  In driving_function_mono_localwfs_sbl (line 96)
  In sound_field_mono_localwfs_sbl (line 92)
  In test_localwfs_sbl (line 159) 

[fietew]

Both testing functions are running, but of course I cannot guarantee that all of the resulting figures look like they should look ;) For test_localwfs_sbl(0) I get a lot of the following warnings, which doesn't happen for test_localwfs_vss(0):

The plane wave decomposition used for LWFS-SBL contains plane waves with propagation angles uniformly distributed on the unit circle. For some (non-enclosing) secondary source distributions, e.g. linear, the propagation directions of some plane waves lead to an empty x0 after the secondary source selection.

[hagenw]

Should we then disable that warning in the test function or maybe directly in driving_function_mono_wfs_pwd()?

[hagenw]

Under SFS_analysis we had some outdated local WFS functions as well. I updated them to

freq_response_localwfs_sbl()
freq_response_localwfs_vss()
time_response_localwfs_sbl()
time_response_localwfs_vss()

If I run the following code, it results in the attached figures, presented in the same order as the execution of the code. Could you have a look at them if they look all right to you.

conf = SFS_config;
freq_response_localwfs_sbl([0 0 0],[0 2.5 0],'ps',conf)
freq_response_localwfs_vss([0 0 0],[0 2.5 0],'ps',conf)
time_response_localwfs_sbl([0 0 0],[0 2.5 0],'ps',conf)
time_response_localwfs_vss([0 0 0],[0 2.5 0],'ps',conf)

[fietew]

The results for freq_response_localwfs_sbl([0 0 0],[0 2.5 0],'ps',conf) do look as expected, but the large magnitude at low frequency is caused by numerical instabilities. For the time-domain implementation of the point source driving function, we used crossover between conventional WFS (low frequencies) and LWFS-SBL (high frequencies) to circumvent this instabilities. Shall this also be implemented for the mono version?

[hagenw]

I would say yes. It's always a good ideas if both domains are implemented as similar as possible.

[fietew]

What is the purpose of the functions in SFS_analysis? Are they meant to just show a good parametrised example? At the moment, the number of plane waves conf.localwfs_sbl.Npw is estimated based on the given f, if not given explicitly. It grows linearly with frequency, so it might very high for some iterations in freq_response_localwfs_sbl. Should we fix the parameters inside freq_response_localwfs_sbl?

[hagenw]

If you think there could be better parameter it would be a good idea to fix this in directly freq_response_localwfs_sbl().

For me the functions in SFS_analysis have two purposes: 1) They show how you could also calculate impulse and frequency response of your system, without using dummy_irs(). 2) The freq_response_* functions are a good way to compare the imp and mono implementations of the toolbox as they calculate the frequency response using the mono functions and you can do the same with dummy_irs() using the imp functions.

[hagenw]

For the driving functions it would be a good idea to add a reference.

For driving_function_imp_localwfs_sbl_ps() we could add your EUSIPCO paper and for driving_function_imp_localwfs_sbl_pw() the one from Nara, correct? What about the mono-frequent driving functions?

[fietew]

For driving_function_imp_localwfs_sbl_ps() we could add your EUSIPCO paper and for driving_function_imp_localwfs_sbl_pw() the one from Nara, correct?

No, Nara's paper takes another approach. I will add a reference.

What about the mono-frequent driving functions?

The monochromatic implementation covers mostly all aspects of the time-domain implementation, so it can be regarded as a simple inverse fourier transform of the latter.

[fietew]

Could we do a minor release after merging the PR?

[hagenw]

Yes, this shouldn't be a problem. Do you have a few days left as there are a few problems under Octave currently with this PR and I would like to include the Voronoi PR as well.

[fietew]

Do you have a few days left as there are a few problems under Octave currently with this PR and I would like to include the Voronoi PR as well.

I will run the examples in Octave, too. Is there a specific issue with some of the functions? The Voronoi PR will still take a while. We have to cover all the (pathological) cases as done for findconvexcone.

[hagenw]

The problems under Octave are related to this Ocatve bug. The output of linkwitz_riley() has not completely symmetric values, which will raise an could not pair all complex numbers error in cplxpair() which is called by zp2sos(). The asymmetry can look like this:

-4.64393188087332e+03 - 5.65864736943455e+03i
-4.64393188087333e+03 + 5.65864736943455e+03i

where the very last value of the real part of the complex number is not identical.

One workaround would be to add the following at the end of linkwitz_riley()

if isoctave
    % Avoid cplxpair error due to numerical problems in Octave,
    % see http://savannah.gnu.org/bugs/?49996
    z = round(1e9*z)*1e-9;
    p = round(1e9*p)*1e-9;
    k = round(1e9*k)*1e-9;
end

but I was not able to test if this would still give the same result at the end as under Matlab (as I have no Matlab today).

[fietew]

There is an additional octave bug. Try:

zp2sos([1],[-1; 1], 1)
zp2sos([1; -1],[-1; 1], 1)
zp2sos([1; 0],[-1; 1], 1)

The first means H(z) = (z-1)/((z-1)*(z+1)) = (z^-1 - z^-2)/(1-z^-2) The second means H(z) = ((z-1)*(z+1))/((z-1)*(z+1)) = (1-z^-2)/(1-z^-2) The third means H(z) = ((z-1)*z)/((z-1)*(z+1)) = (1 - z^-1)/(1-z^-2)

Octave:

ans =

   1  -1   0   1   0  -1
ans =

   1   0  -1   1   0  -1
ans =

   1  -1   0   1   0  -1

MATLAB:

ans =

     0     1    -1     1     0    -1
ans =

     1     0    -1     1     0    -1
ans =

     1    -1     0     1     0    -1

It is not the case, that Octave just uses a different convention for the output, because the second output should then also follow this convention. The third output is equal to the first one, which should not be the case.

[hagenw]

Your changes seem to have fixed the problems.

conf = SFS_config;
freq_response_localwfs_sbl([0 0 0],[0 2.5 0],'ps',conf)

Matlab: Octave:

[hagenw]

Regarding the zpsos() problem. I haven't looked into the details how you have fixed it, but do you think this zpsos() call in `bilinear_transform() is affected as well?

[hagenw]

I ran the test_localwfs_sbl() under Octave and Matlab and it worked well in both cases. I added another commit to disable the zero secondary source selection warning in the pwd driving functions. If you think that is not a good idea we can revert this. Otherwise this PR is ready in my opinion and I would clean up the commit history and merge if this is fine with you.

[fietew]

FYI: https://savannah.gnu.org/bugs/?51936

Regarding the zpsos() problem. I haven't looked into the details how you have fixed it, but do you think this zpsos() call in `bilinear_transform() is affected as well?

I think any usage for zp2sos() is affected by this problem, because it seems to occur whenever you have an odd number of zeros.

[hagenw]

I have updated the commits, if you agree with them I will merge it into master. I would like to put also some entries in the NEWS file. Would the following fit?

- add monochromatic implementation of LWFS using spatial bandwidth-limitation
- add monochromatic circular expansion functions for ps and pw
- add function for conversion from circular to plane wave expansion
- add freq_response_* and time_response_* for all LWFS methods

[fietew]

I have updated the commits, if you agree with them I will merge it into master. I would like to put also some entries in the NEWS file. Would the following fit?

- add monochromatic implementation of LWFS using spatial bandwidth-limitation
- add monochromatic circular expansion functions for ps and pw
- add function for conversion from circular to plane wave expansion
- add freq_response_* and time_response_* for all LWFS methods

Yes