Closed narahahn closed 8 years ago
As I have no experience with the SOS filtering stuff, could you maybe create a branch and implement an example of your solution? The only thing we have to consider is to test it also in Octave. If I remember correctly the first implementation of the NFC-HOA stuff in the time-domain was only running under Matlab.
I just created a branch and updated the SOS implementation of 2.5D NFC-HOA driving functions.
Previously, the matlab function zp2sos
was used with the pole-zeros in the Laplace domain. @fs446 pointed out that zp2sos
is actually intended for z-domain representations (either pole-zeros or transfer function coefficients).(http://de.mathworks.com/help/signal/ref/zp2sos.html) Mathematically, the results are identical but there might be numerical problems.
Now the driving functions are computed as follows:
z
) of the spherical Hankel function is computed. (s-domain)z
. (s-domain)default
driving function) or by matched z-transform (matchedz
driving function). Now we are in z-domain!sos
) are computed by zp2sos
.sos
is normalized so that the gain at fs/2 is 0 dB.pulse
) is filtered with sos
SFS_start;
conf = SFS_config;
conf.dimension = '2.5D';
conf.secondary_sources.geometry = 'circle';
conf.secondary_sources.size = 3;
conf.nfchoa.order = 16;
conf.plot.usedb = true;
x0 = secondary_source_positions(conf);
X = [-2 2];
Y = [-2 2];
Z = 0;
t = 200;
% virtual plane waves
xs = [0,-1,0];
src = 'pw';
conf.driving_functions = 'default';
sound_field_imp_nfchoa(X,Y,Z,xs,src,t,conf);
title('bilinear transform')
% print('-dpng','d_pw_bilinear')
conf.driving_functions = 'matchedz';
sound_field_imp_nfchoa(X,Y,Z,xs,src,t,conf);
title('matched-z transform')
% print('-dpng','d_pw_matchedz')
The difference between two driving functions is almost negligible.
I got the following error, when I try to execute your above code in Octave:
>> sound_field_imp_nfchoa(X,Y,Z,xs,src,t,conf);
error: bilinear: must have at least as many poles as zeros in s-plane
error: called from
bilinear at line 91 column 5
driving_function_imp_nfchoa_pw at line 102 column 14
driving_function_imp_nfchoa at line 93 column 13
sound_field_imp_nfchoa at line 94 column 3
If I remember correctly there is a slight difference in the behavior of the bilinear
function in Octave compared to Matlab. I will have a look at this next week.
I realized that bilinear
in Octave is slightly different than Matlab. The last input in Octave is the sampling period T
while in Matlab it is the sampling frequency fs
.(http://octave.sourceforge.net/signal/function/bilinear.html)
Could you check if it works now?
I computed the zeros of the spherical Hankel functions WITH and WITHOUT using the multi-precision toolbox. The zeros are almost the same below 26
and the distance in the z-domain is in the order of 10^7
.
However, for orders equal or greater than 26
, the zeros computed WITH the multi-precision toolbox seem very strange. N-2
multiple zeros appear at (1,0)
and the other two zeros are in odd positions.
The following figures show the driving function for each modes in the time-domain.
As can be seen in the first figure, a discontinuity exists between n=25
and n=26
. The results WITHOUT the multi-precision seem more reasonable. In the latter, however, some poles of the driving function (=zeros of the spherical Hankel function) begin to lie outside the unit circle in the z-domain and the driving function becomes unstable.
The synthesized sound fields are compared here. NFC-HOA of order 186
was simulated for 64
secondary sources. Since the order is much higher than the anti-aliasing order (31
), the result is likely to resemble WFS. However, this is only observed when the zeros are computed WITHOUT the multi-precision toolbox.
Due to the instability, the sound field for t=1000
look like this:
Result WITH multi-precision looks like typical band-limited NFC-HOA.
In the latter, I modified sphbesselh_zeros
so that all zeros (except `n=1,2') are computed by the multi-precision toolbox.
Long story short:
86
has to be lowered.I have computed the zeros of the spherical Hankel function with Python
import mpmath as mp
mp.dps = 30; mp.pretty = True
order = 75
B = np.zeros(order)
for n in range(order):
B[n] = np.math.factorial(2*order-n)/(np.math.factorial(order-n)*np.math.factorial(n)*2**(order-n))
B = B[::-1]
z = mp.polyroots(B, extraprec=300, maxsteps=300)
The resulting pole locations (matched z-transform) seem to be plausible for order N=75
Another result from the Python simulation for order N=150. Now zeros (which will later become poles in NFC-HOA) are located outside of the unit circle.
I have added the Jupyter notebook to sfstoolbox/data repository https://github.com/sfstoolbox/data/blob/master/sphbesselh_zeros/sphHankel_zeros.ipynb
I created a pull request (#65) to discuss the changes proposed in the nfchoa_sosfilter
branch.
In this issue we continue the discussion on the problems with finding the zeros of the Bessel function and related numerical problems.
I could reproduce the results from @narahahn above, there are really some problems with the zeros calculated by the Multiprecission-Toolbox. As there are also numerical problems in python as well, we should maybe have a look at alternative ways to calculate the zeros. I had a first start on this, see https://github.com/sfstoolbox/sfs-deprecated-code/tree/master/hankel_zeros, but never really tried it.
I think the problem is at least twofold:
Regarding the first problem, there is a recursion formula for the numerator coefficients in http://iem.kug.ac.at/fileadmin/media/iem/altdaten/projekte/acoustics/awt/winkel/pomberger.pdf, page 43. A while ago, I implemented a function utilizing this recursion formula. Dunno, if it leads to more stable results.
function [hankel, derivative] = sphhankel_laplace(order)
hankel = struct(...
'zeros',[],...
'poles',[],...
'scale',[],...
'delay',[],...
'nominator',[],...
'denominator',[]);
hankel = repmat(hankel,order+1,1);
for n=0:order
% nominator
hankel(n+1).nominator = zeros(1, n+1);
for k=0:n-1
hankel(n+1).nominator(k+1) = ...
((2*n-k-1)*(2*n - k))/(2*(n-k))*hankel(n).nominator(k+1);
end
hankel(n+1).nominator(n+1) = 1;
end
if nargout >= 2
derivative = hankel;
for n=0:order
% nominator
derivative(n+1).nominator = zeros(1, n+2);
if n == 0
derivative(1).nominator = -hankel(2).nominator;
else
derivative(n+1).nominator(1:n+1) = (n+1)*hankel(n+1).nominator;
for k=2:n+1
derivative(n+1).nominator(k+1) = derivative(n+1).nominator(k+1) + ...
hankel(n).nominator(k-1);
end
end
end
for n=0:order
% flip nominator polynoms
derivative(n+1).nominator = derivative(n+1).nominator(end:-1:1);
% zeros
derivative(n+1).zeros = roots(derivative(n+1).nominator);
% denominator
derivative(n+1).denominator = zeros(1, n+3);
derivative(n+1).denominator(1) = 1;
% poles
derivative(n+1).poles = roots(derivative(n+1).denominator);
% additional scaling factor
derivative(n+1).scale = 1j.^(n+1);
% additional delay
derivative(n+1).delay = 1;
end
end
for n=0:order
% flip nominator polynoms
hankel(n+1).nominator = hankel(n+1).nominator(end:-1:1);
% zeros
hankel(n+1).zeros = roots(hankel(n+1).nominator);
% denominator
hankel(n+1).denominator = zeros(1, n+2);
hankel(n+1).denominator(1) = 1;
% poles
hankel(n+1).poles = roots(hankel(n+1).denominator);
% additional scaling factor
hankel(n+1).scale = -1j^n;
% additional delay
hankel(n+1).delay = 1;
end
end
By using the recurrence relation, the driving functions can be computed up to order 150 without NaN
or Inf
errors. Unfortunately, the stability is still a problem.
The discussion on the second order section filters is continued at #65. The problems with the Multiprecission Toolbox is fixed with #70. The discussion on how to find the zeros of the spherical Hankel function is continued at #71.
Hi @narahahn and @hagenw.
What precision did you use in computations with Multiprecision Toolbox? Finding zeros for high-order spherical Hankel function requires ever-growing level of precision (at least for the code above - factorials grow quickly).
Default precision in toolbox is quadruple (~34 decimal digits ) and obviously it was not enough for the computations. The good thing is that toolbox is able to work with arbitrary level of precision. User just need to call mp.Digits(XXX) before doing any computations with toolbox (XXX is number of decimal digits to use).
I am author of Multiprecision Toolbox.
The SOS coefficients in the Laplace domain is mapped to the SOS coefficients in the z-domain by using the bilinear transform. The latter are stored in two cell arrays, b and a. The filtering is then performed in a for-loop. I think we can do this in a simpler way: