Second-order section structure filters for NFC-HOA driving function

[narahahn]

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:

• transform the zeros in the s-domain to the z-domain using the bilinear function
• compute the SOS coefficients using zp2sos (SOS is a matrix)
• filter the input (pulse) using sosfilt

[hagenw]

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.

[narahahn]

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:

• The zeros (z) of the spherical Hankel function is computed. (s-domain)
• The poles and zeros of the driving function is obtained by scaling z. (s-domain)
• The poles and zeros are transformed into the z-domain either by bilinear transform (default driving function) or by matched z-transform (matchedz driving function). Now we are in z-domain!
• Second-order section coefficients (sos) are computed by zp2sos.
• Each numerator of the sos is normalized so that the gain at fs/2 is 0 dB.
• The impulse (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.

[hagenw]

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.

[narahahn]

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?

[narahahn]

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:

• There is something odd in computing the zeros using the multi-precision toolbox.
• The current NFC-HOA driving function is unstable for high orders. It seems that the limit of 86 has to be lowered.

[spors]

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

[spors]

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.

[spors]

I have added the Jupyter notebook to sfstoolbox/data repository https://github.com/sfstoolbox/data/blob/master/sphbesselh_zeros/sphHankel_zeros.ipynb

[hagenw]

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.

[hagenw]

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.

[fietew]

I think the problem is at least twofold:

• How to compute the coefficients of the denominator and numerator polynomials correctly, as factorials of large arguments are needed.
• How to compute the zeros of the polynomials correctly

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

[narahahn]

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.

[hagenw]

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.

[advanpix]

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.