More stable algorithm to calculate zeros of spherical hankel function

[fietew]

Follow up on #71:

Based on the publication of @narahahn and @spors.

This is basically a port of a scipy routine.

@narahahn: Could you try out some parametrisation and check, if it is working as expected?

[fietew]

@hagenw: Is it okay to remove the other methods?

[hagenw]

The old version of the function returned z and p both as column vectors. Your new implementation returns z as row vector and p as column vector. Is this desired or a bug?

If the new method works better than the old ones, I see no big deal in replacing them as all our methods we implemented so far were buggy.

[hagenw]

Is the code ready for me to have a look at it?

It would be cool if you could add a test function for this, maybe doing something like the figures in #57 for example those two or Figure 3 in @narahahn publication

[fietew]

It is ready. I will add a test function.

[hagenw]

The test function seems to work ok under Octave, but I had to disable the legend as it covered the whole figure. I guess you have split it to six different subplots in order to not have to many overlapping points, or are the subplots showing different things? For me they looked very similar, here is the first one of them, I guess this is the desired result?

[hagenw]

I was looking to update the reference as the Campos and Calderón paper was published in 2012, see http://doi.org/10.1155/2012/873078 The problem is that they reformulated some of their equations and I was not able to fit them easily to the ones in https://arxiv.org/pdf/1105.0957.pdf, which should be the paper you referenced. Could you have a look at it and propose if we should stay with the arxiv version or switch to the other one? I would then do the update to the comments.

I'm also wondering if we should add a link to the scipy routine as you used it as a template.

Otherwise, I'm really happy that this seems to be finally solved as it bothered me since the first version that included NFC-HOA.

[fietew]

The test function seems to work ok under Octave, but I had to disable the legend as it covered the whole figure. I guess you have split it to six different subplots in order to not have to many overlapping points, or are the subplots showing different things?

Each subplot shows a subset of orders. Unfortunately the colors are missing in octave. The graphs shouldn't be all blue, but the automatic coloring is a MATLAB feature, I guess

For me they looked very similar, here is the first one of them, I guess this is the desired result?

Yes, the result is as expected.

I was looking to update the reference as the Campos and Calderón paper was published in 2012, see http://doi.org/10.1155/2012/873078 The problem is that they reformulated some of their equations and I was not able to fit them easily to the ones in https://arxiv.org/pdf/1105.0957.pdf, which should be the paper you referenced. Could you have a look at it and propose if we should stay with the arxiv version or switch to the other one? I would then do the update to the comments.

Yes, I also had a look at the journal paper, which covers a more generalised problem and is harder to understand. I would stay with arxiv version since the formulas in the campos_zeros correspond directly to equations in the paper.

I'm also wondering if we should add a link to the scipy routine as you used it as a template.

Yes, that's a good idea.

[hagenw]

Ok, I updated the references accordingly.

For the test function, could you maybe send me a screen shot of the Matlab result, that I can try to get the same with Octave.

[hagenw]

I added a numerical test case as well for the order 1, 10, and 100 comparing to the results from the python code. This is not perfect as the outcome depends on the ordering of the zeros in the resulting vector, which can differ for different calculation methods, but I thought it is still better than nothing.

[hagenw]

I also did a NFC-HOA calculation in order to compare with Hahn and Spors (2017) Fig. 8:

conf = SFS_config;
conf.nfchoa.order = 30;
sound_field_imp_nfchoa([-2 2],[-2 2],0,[1.5 1.5 0],'ps',1/conf.c,conf)

conf.nfchoa.order = 120;
sound_field_imp_nfchoa([-2 2],[-2 2],0,[1.5 1.5 0],'ps',1/conf.c,conf)

It looks very similar, but we have the different of that prominent circular contribution which is not part of the plot in the paper. If I remember correctly we had that discussion already somewhere else, hadn't we? Is this due to the inclusion/exclusion of the lowest order?

[fietew]

For the test function, could you maybe send me a screen shot of the Matlab result, that I can try to get the same with Octave.

I am pretty sure, that the coloring is not relevant, but I will explicitly define the colors in the script and send you a screen shot.

Is this due to the inclusion/exclusion of the lowest order?

What are you referring to? What is the lowest order? What is included/excluded?

[fietew]

The plot under MATLAB R2015a:

[hagenw]

The colors in the plot work now also under Octave, but the legend was still far to big. I added a commit that replaces the legend with an updated title, indicating the shown orders for every plot.

[hagenw]

Coming back to the computed figure and the comparison to Hahn and Spors (2017) Fig. 8. Forget what I have set about the order. But do you have an idea why we have the obvious difference between those plots (of course besides the difference in amplitude)?

[fietew]

Are you using normalisation for the plot? Also the range of the colormap in Hahn and Spors (2017) Fig. 8. is [-2,2].

[hagenw]

I'm using normalisation, otherwise there is no energy in the plot. But you can do the following to get a very similar energy as in Hahn and Spors (2017) Fig. 8:

conf = SFS_config;
conf.nfchoa.order = 30;
conf.plot.caxis = [-0.3 0.3];
sound_field_imp_nfchoa([-2 2],[-2 2],0,[1.5 1.5 0],'ps',1/conf.c,conf)

will result in

This looks similar, but still we have a more pronounced visible circle as indicated by the arrow in the figure.

[fietew]

Was #56 the discussion you meant?

[fietew]

Is it possible, that you have chosen a different number of loudspeakers? In the paper, 60 has been chosen. For 64 loudspeakers, one loudspeaker is exactly at 45 degrees, which is the virtual point source direction in your plots. For 60 loudspeakers, 45 degrees lies between two loudspeakers.

[hagenw]

Ah, brilliant, that's it:

conf = SFS_config;
conf.nfchoa.order = 30;
conf.secondary_sources.number = 60;
conf.plot.caxis = [-0.3 0.3];
sound_field_imp_nfchoa([-2 2],[-2 2],0,[1.5 1.5 0],'ps',1/conf.c,conf)

and we will get the same result as in the paper.

I would say this is then nearly ready. There is only this discussion on the Bessel polynomial left, @narahahn should give his go, and after a small commit clean up we could merge.

[hagenw]

I made another test and tried to improve Fig. 3.14 from my thesis and it worked quite well.

Original NFC-HOA part of the figure:

New NFC-HOA part of the figure as calculated with this branch:

@fietew: From my side, I have no further change requests for this branch and would like to start cleaning the commit history up, if you agree.

[fietew]

I agree.