More stable algorithm to calculate zeros of spherical hankel function #152

fietew · 2017-06-23T12:06:20Z

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 · 2017-06-23T12:07:01Z

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

hagenw · 2017-06-26T11:11:28Z

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 · 2017-07-06T09:59:36Z

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 · 2017-07-06T11:18:40Z

It is ready. I will add a test function.

hagenw · 2017-07-06T14:55:32Z

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 · 2017-07-06T16:02:55Z

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 · 2017-07-07T06:50:20Z

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 · 2017-07-07T08:19:57Z

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 · 2017-07-07T11:54:28Z

SFS_general/sphbesselh_zeros.m

+% Zeros of nth-order ordinary Bessel polynomial are the same as for
+% exp(1/x)*K_{n+0.5}(1/x), where K_v is the modified Bessel function of second
+% kind. We, hence, define the target function and its first derivative as
+f = @(x) besselk(order+0.5, 1./x, 1);  % modified Bessel function scaled by exp


I would remove the comment % modified Bessel function scaled by exp as it my lead a user to expect the exp() function to be part of that line. In addition, I think the comment provided before is sufficient.

It took me a while, to figure out that besselk(nu, x, 1) is equal to besselk(nu,x)*exp(x). This is why I added this comment. But if it misleading, we can remove it.

I was a little bit confused by it.
Maybe it helps if we state explicitly the equation in the first comment. Something like

% Zeros of nth-order ordinary Bessel polynomial ... are the same as for modified % Bessel function of second kind K_v: % exp(1/x) * K_{n+0.5}(1/x) = ... % Hence, we define the target function and its first derivative as

Where you have to fill the ... part. Is it also second kind and then Y?

Could you help to fill the ... please.

According to https://en.wikipedia.org/wiki/Bessel_polynomials#Definition_in_terms_of_Bessel_functions
the nth-order ordinary bessel polynomial is denoted by y_n, which is also the spherical Neumann function in our nomenclature. I would suggest:

% Zeros of nth-order ordinary Bessel polynomial y_n are the same as for % the exponentially scaled modified Bessel function of second kind K_v: % \sqrt{2pi/x} * exp(1/x) * K_{n+0.5}(1/x) = y_n % Hence, we define the target function and its first derivative as

hagenw · 2017-07-07T12:30:45Z

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 · 2017-07-07T12:32:29Z

SFS_general/sphbesselh_zeros.m

-    B(n+1) = factorial(2*order-n)/(factorial(order-n)*factorial(n)*2^(order-n));
+
+function x = newton(f,fp,x0,TOL,MAXITER)
+    % Newton-Rapheson method for approximation of a single root


Where comes this implementation from?
When I'm looking at scipy's newton implementation it looks a lot more complicated.

The scipy implementation does handle many special cases and has some extensions. Here, the most simple version of the Newton-Rapheson method is implemented.

Now, I found also in the scipy version. Then maybe the easiest thing is to add not a link to the function header, but simply go with the current help text for newton().

I can't follow you. What do you want me to do?

Nothing :)
Just a confirmation that it is fine to not add a reference link to the scipy code for the Newton method.

hagenw · 2017-07-07T12:36:49Z

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 · 2017-07-08T15:42:33Z

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 · 2017-07-10T08:38:27Z

The plot under MATLAB R2015a:

hagenw · 2017-07-20T18:32:51Z

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 · 2017-07-20T18:34:58Z

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 · 2017-07-20T18:44:37Z

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

hagenw · 2017-07-20T19:03:12Z

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 · 2017-07-20T20:11:51Z

Was #56 the discussion you meant?

fietew · 2017-07-20T20:19:52Z

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 · 2017-07-20T21:47:32Z

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 · 2017-07-28T09:17:05Z

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 · 2017-07-29T09:31:51Z

I agree.

fietew assigned narahahn Jun 23, 2017

hagenw mentioned this pull request Jul 7, 2017

References in function headers #156

Closed

hagenw reviewed Jul 7, 2017

View reviewed changes

Fiete Winter and others added 2 commits July 31, 2017 16:08

Improve sphbesselh_zeros() for high orders

a338282

Add validation script for sphbesselh_zeros()

b52a9b0

hagenw force-pushed the sphbesselh_zeros branch from e2222a0 to b52a9b0 Compare July 31, 2017 15:10

hagenw merged commit cb3a7c1 into master Aug 14, 2017

hagenw deleted the sphbesselh_zeros branch August 14, 2017 08:23

hagenw mentioned this pull request Aug 22, 2017

Calculating zeros of spherical Hankel function #71

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

More stable algorithm to calculate zeros of spherical hankel function #152

More stable algorithm to calculate zeros of spherical hankel function #152

fietew commented Jun 23, 2017

fietew commented Jun 23, 2017

hagenw commented Jun 26, 2017

hagenw commented Jul 6, 2017

fietew commented Jul 6, 2017

hagenw commented Jul 6, 2017

hagenw commented Jul 6, 2017

fietew commented Jul 7, 2017 •

edited

Loading

hagenw commented Jul 7, 2017

hagenw Jul 7, 2017

fietew Jul 7, 2017

hagenw Jul 7, 2017

hagenw Jul 20, 2017

fietew Jul 21, 2017

hagenw commented Jul 7, 2017

hagenw Jul 7, 2017

fietew Jul 8, 2017

hagenw Jul 20, 2017

fietew Jul 20, 2017

hagenw Jul 20, 2017

hagenw commented Jul 7, 2017

fietew commented Jul 8, 2017 •

edited

Loading

fietew commented Jul 10, 2017

hagenw commented Jul 20, 2017

hagenw commented Jul 20, 2017

fietew commented Jul 20, 2017

hagenw commented Jul 20, 2017 •

edited

Loading

fietew commented Jul 20, 2017

fietew commented Jul 20, 2017 •

edited

Loading

hagenw commented Jul 20, 2017

hagenw commented Jul 28, 2017

fietew commented Jul 29, 2017

More stable algorithm to calculate zeros of spherical hankel function #152

More stable algorithm to calculate zeros of spherical hankel function #152

Conversation

fietew commented Jun 23, 2017

fietew commented Jun 23, 2017

hagenw commented Jun 26, 2017

hagenw commented Jul 6, 2017

fietew commented Jul 6, 2017

hagenw commented Jul 6, 2017

hagenw commented Jul 6, 2017

fietew commented Jul 7, 2017 • edited Loading

hagenw commented Jul 7, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hagenw commented Jul 7, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hagenw commented Jul 7, 2017

fietew commented Jul 8, 2017 • edited Loading

fietew commented Jul 10, 2017

hagenw commented Jul 20, 2017

hagenw commented Jul 20, 2017

fietew commented Jul 20, 2017

hagenw commented Jul 20, 2017 • edited Loading

fietew commented Jul 20, 2017

fietew commented Jul 20, 2017 • edited Loading

hagenw commented Jul 20, 2017

hagenw commented Jul 28, 2017

fietew commented Jul 29, 2017

fietew commented Jul 7, 2017 •

edited

Loading

fietew commented Jul 8, 2017 •

edited

Loading

hagenw commented Jul 20, 2017 •

edited

Loading

fietew commented Jul 20, 2017 •

edited

Loading