sfstoolbox / sfs-matlab Goto Github PK

This is a follow up on #57.

We calculate the zeros of the spherical Hankel function with sphbesselh_zeros(), which works only up to an order of 85 due to numerical problems.

For higher orders we have to find a solution, which will be discussed here.

Rename these functions in a more meaningful way. For example, we could replace 'mono' with 'f' and 'imp' with 't', to indicate the different domains we are working in.

In addition to this it is a very confusing, that there exist a folder for the time domain functions, a folder for the ir stuff and a folder for the binaural synthesis. Maybe this can also be done in a better way.

I experienced that the tapering window for focused has to be shorter/start earlier (around 50%) then for non-focused sources (around 30%). In the context of local wfs using virtual secondary sources (vss), it might be necessary to use a different tapering window for the vss and for the real loudspeakers ( secondary sources ). Shall we put it in the code again or is it an overkill?

Find a way to calculate the aliasing frequency in a better way.

Maybe try to use the geometric approach with distance between individual wave fronts.

The current formulation of the line source driving function is used for 2D and 3D WFS, but it only works in 2D.
To make it work, I think x0 - xs should be replaced with a vector perpendicular to the line source, that is x0 - xs - <x0-xs,v> v where xs is a point on the line and v a unit vector along the line.

For the mono frequent wave fields the frequency f is the equivalent. There f is directly given as an parameter to the functions.

Hence, we should get rid of conf.frame completely by replacing it by the time t (in ms) as an input parameters to the function. Then we have to ensure that at time t=0 the wave field starts to come out of the secondary sources.

The direction of the focused source (and therefor the secondary source selection) is done via the reference point at the moment. Replace this method with an additional vector.

While introducing our new concept to include non-regular grids, we have not considered/tested the gnuplot functionalities. At the moment the gp_save_matrix needs vectors for x and y and yields an error if I execute

[P,x,y,z] = sound_field_mono(X,Y,Z,[xs,ns,1],'ps',1,f,conf);
gp_save_matrix('sound_field_line_source.dat',x,y,real(P));

At the moment the first loudspeaker of your array has to be on the x-axis
(which means phi=0). If you have another setup (like we have in Pinta) you
have to manually edit the secondary_source_positions.m function in order to
get the first loudspeaker at the desired location:

phi = linspace(pi/2,(5/2-2/nls)*pi,nls)';

At the moment we can only calculate a tapering window for 2D arrays like circle, line, box. For 3D the edge detection is not working, which is a prerequisite for calculating a tapering window. The challenge is to find a edge detection algorithm that works for 3D and 2D geometries. If this is not possible maybe we should detect first, if we have a 2D or 3D geometry and apply different algorithms.

This will not only made it more easy for real life scenarios, but enable us also to be consitent with a 3D spherical grid. Because in this case we have to set the number of sources anyhow to get the right grid.

At the moment only the mono-frequent driving function is available for 2D WFS

At the moment get_ir is interpolating always

We have started to work on this by returning the delay from wfs_preequalization() as well as from driving_function_imp_wfs(), see #62.

The same behavior should be added to delayline(). As this function is updated in general at the moment (see #50), this feature is postponed after finishing #50.

SOFA is a new upcoming common format for impulse responses:
http://sourceforge.net/projects/sofacoustics/

The number of taps for the WFS FIR prefilter is currently hard-coded in wfs_fir_prefilter() to Nfilt = 128.
If one is interested e.g. in the delay introduced by WFS reproduction, the filter length appears as a magic number.

I'd suggest to store the filter length in conf.wfs, maybe as hpreFIRnTaps (N) or hpreFIROrder (N-1).

(Feel free to assign the issue to me.)

Test the performance of all the demanding Toolbox functions and optimize the code.

conv is slow and needs a loop for convolving a multichannel impulse response.
Have alook at fftfilt and convolve from the LTFAT Toolbox for the implementation

Could you add a description for all the input and output parameters in get_shelve_lagrange().

How to handle the switching between different implementations of driving functions?
For example, if we like to switch to the ones from TU Delft?

This would be nice for aliasing_frequency() and is needed in tapering_window()

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

I have added several TODO messages in the text of some of the new localsfs functions and the SFS_config_example.m file. Could have a look at those and add the corresponding stuff to it.

It should use hanning(2_onset+1) and not hanning(2_onset)

In addition we should remove the leading and ending zero of the window. Note, that this is at the moment explicitly done in the tapering_window() function.
Further we could think of switching to chebwin() instead of hanning() and maybe renaming this function.

#77 added the ability of modal weighting of NFC-HOA orders to the time domain driving function of NFC-HOA.

This cannot be easily added to the current monochromatic implementation of NFC-HOA, but should be able to be added after implementing the new NFC-HOA structure as proposed in #63.

By executing the test function for the IIR pre-filter from current master (bdd0e4a) in Matlab R2013a

>> test_wfs_iir_prefilter

I got the following error:

Error using fdesign.arbmag/set
There is no enumerated value named 'Nb,Na,F,A'.

Error in /usr/local/MATLAB/R2013a/toolbox/signal/signal/@fdesign/@abstracttypewspecs/schema.p>set_specificationtype (line 65)

I suggest to enable two-channel headphone compensation. I think this could be done easily by allowing convolution() to accept two matrices as input that have one dimension of equal length.
I could do this myself if there are no problems likely to arise from this change?

As you can see in the above image the sound of the plane wave has some strange wobbles if the plane wave goes to the direction [1 -1 0] instead of [0 -1 0]. Have you an idea if this is correct?

By the way I renamed test_virtualsecondarysources.m to test_localwfs.m and added more array setups.

In driving_function_imp_nfchoa.m, there is a if-request to solve the problem of wrong propagation direction. I found that this comes from the error in computing the input for spatial IFFT (here). According to the Eq. (5) and (6) of Spors et al, 2011 the exponential term should have negative exponent.
Should I handle this issue myself?

Check if we can set conf.N to just 2048, because the default value of 2^15 causes a lot of computational time.

Try the following

conf = SFS_config;

conf.plot.realloudspeakers = true;

conf.tapwinlen = 1.0;
conf.secondary_sources.number = 36;
conf.secondary_sources.geometry = 'circular';

xs = [0, 3, 0];
x0 = secondary_source_positions(conf);
[x0, idx] = secondary_source_selection(x0, xs, 'ps');
x0 = secondary_source_tapering(x0, conf);
figure,
subplot(2,2,1);
draw_loudspeakers(x0, [1 1 0], conf);
xlim([-1.6,1.6]);
ylim([-1.6,1.6]);
title(sprintf('number of secondary sources: %d', size(x0,1))); 
subplot(2,2,2);
plot(x0(:,7), 'b');

xs = [0, -3, 0];
x0 = secondary_source_positions(conf);
x0 = secondary_source_selection(x0, xs, 'ps');
x0 = secondary_source_tapering(x0, conf);
subplot(2,2,3);
draw_loudspeakers(x0, [1 1 0], conf);
xlim([-1.6,1.6]);
ylim([-1.6,1.6]);
title(sprintf('number of secondary sources: %d', size(x0,1))); 
subplot(2,2,4);
plot(x0(:,7), 'r');

Although everything is symmetric, the number of selected sources is different. I think this is related to the usage of epsin https://github.com/sfstoolbox/sfs/blob/master/SFS_general/secondary_source_selection.m#L124 . While replacing it with 0 does make a difference, replacing it with -eps leads to the same result for both setups. Note, that the range for azimuth of the the selected secondary sources lies between phi_s +- 60 degrees, where phi_s is the azimuth of the virtual point source. Due to the azimuthal grid of 10 degrees resolution two secondary sources are exactly on this bound.

I only changed the config settings in that script to the new structure and it stoped working with the following message:

>> test_adaptivesecondarysources
Progress: [..........] 100%
Error using *
Inner matrix dimensions must agree.

Error in virtual_secondary_source_positions (line 267)
    xv(:,1:3) = repmat(xl, nls, 1) + repmat(Rd*nd, nls, 1) + x'*ndorth;

Error in test_adaptivesecondarysources (line 53)
  x0 = virtual_secondary_source_positions([],xs,src,conf);

Could you have a look at what goes wrong.

The dummy_irs() function leads to other results as the old irs function.
Be careful if you are using it.

Look at comments
Clean up code

Please add a description in the comment at the beginnning, where it is indicated. And maybe add some remarks to the calculation. for me it is not clear what is happening.

At the moment the aliasing is not working correctly for SDM, because of the lag of spectral repetitions in the kx domain.

e39a0ae introduces a new handling of delayline() which can now handle the default SOFA impulse response matrices in the form of [M C N], where

M ... number of measurements
C ... number of channels
N ... Number of samples

This has some benefits as for example the sofa_get_data() function returns the data in that format which is much easier to handle for smooth SOFA integration. Before we always used [N C] as the standard input format and several measurements were simply put into channels if needed.
At the moment both ways are supported by the current delayline() implementation, but we could consider to remove the [N C] completely in order to simplify and fasten the code, which is critical for this particular function.

I'm not sure if the simple usage of nearest neighbor points for performing the interpolation is the best approach for impulse responses. Here is an easy example:

x0 = [45 46 46 50; 0 -2 2 0];
xs = [47; 0];
x0_int = findnearestneighbour(x0,xs,3);
figure; plot(x0(1,:),x0(2,:),'xb')
axis square; axis([43 52 -3 3]); hold on
plot(xs(1),xs(2),'xr')
plot(x0_int(1,:),x0_int(2,:),'or','MarkerSize',10)
xlabel('Azimuth / deg')
ylabel('Elevation / deg')

As you can see we want to get the point shown as red cross. The red circles are the three nearest neighbors that are selected for the interpolation. Afterwards it is checked if the nearest and second nearest neighbor have the same azimuth or elevation. If so only these points are used. In this example all three points will be used for interpolation.
The problem is that in this case it would in my opinion make more sense to use just the two points that have both an elevation of 0°. The problem is I don't know if we are able to write an algorithm that does this sort of clever selection that works with every grid and point spacing.
Has anyone some experience on this?

Also add additional comments and references to the driving functions.

Handle also the possibility to switch between different implementations, for example Verheijen/Spors/Völk

As @fietew wrote in #45:

As far as I understood the current implementation, the SOFA file is accessed for each head orientation and for each loudspeaker position separately. This means, that an impulse response (having e.g. the index "42" ) is loaded from the hard disk drive into the RAM several times. If we really want to speed up things, then this architecture has to be reassessed. An alternative workflow could be:

• iterate over the head orientations and the secondary source position in order to get the indices of the impulse responses, which are needed from the SOFA file. This can be done by evaluating the header of the SOFA file, i.e. using nearest_neighbor, etc. (as you have done it in the current implementation)
• the indices are likely to be non-unique, as some head-orientations-ssd-position pairs might need the same impulse responses. However, functions like unique (http://de.mathworks.com/help/matlab/ref/unique.html) help us to remove those duplicates.
• sorting the indices might also help to identify segments of connected indices, i.e. groups of directly subsequent indices, e.g. [2,3,4,5] [7,8,9] [11,12]. As I experienced, that loading the impulse responses index by index is quite slow in SOFA, one could to something like https://github.com/TWOEARS/binaural-simulator/blob/master/src/%2Bsimulator/DirectionalIR.m#L189-L202 .
• after loading the impulse responses, they have to be processed according to the head-orientations-ssd-position pairs they are need for, i.e. distance adjustment, angle interpolation, adding up all loudspeakers.

The major drawback of this approach is, that there might be a huge number of impulse responses which are loaded into the RAM at once. One could define a maximum number of irs loaded in one step and then iterate as long as there are impulse responses left to be processed. However, one would have ensure, that there is always at least one head-orientations-ssd-position pair which can be processed with the current impulse responses in the RAM.

If someone has an idea how to do this or some time left we can work on this, but the priority is not very high at the moment.

conf.t0 works not for all cases. Especially the creating of binaural stimuli with ir_wfs_25d has problems with focused sources and plane waves. This issue includes also the fix_ir_length function.

Hi,

Thank you first for this great project!

I encountered the following error when I was trying to generate an impulse response according to 'QU_KEMAR_anechoic_3m.sofa'.

Commands:

conf = SFS_config_example;
hrtf = SOFAload('QU_KEMAR_anechoic_3m.sofa');
ir = get_ir(hrtf,[0 0 0],[0 0],[rad(30) 0 3],'spherical',conf);
nsig = randn(44100,1);
sig = auralize_ir(ir,nsig,1,conf);
sound(sig,conf.fs);

Error:

Undefined function 'fieldnames' for input arguments of type 'double'.

Error in SOFAload (line 153)
    for field = fieldnames(ObjTemplate)'

NB, we found that SOFA API was not provided directly with SFS, therefore we then dug it out from http://sourceforge.net/projects/sofacoustics/ . Also 'QU_KEMAR_anechoic_3m.sofa' was missing, we found several from the database, they were all trapped in this error.

We tried to debug but couldn't manage to make it work.

Looking forward to your reply.
Cheers!

Hi,

We followed the example to generate impulse responses for SSR.
Commands:

conf = SFS_config_example;
SFS_start;
SOFAstart;
irs = SOFAload('QU_KEMAR_anechoic_3m.sofa');
brs = ssr_brs_wfs([0,0,0],pi/2,[0,3,0],'ps',irs,conf);
wavwrite(brs,44100,16,'brs_set_for_SSR.wav');

However, wave file looks wrong (all 720 channels)

We also experimented with ssr_brs_point_source(), which actually suits our need the most.

brs = ssr_brs_point_source([0,0,0],pi/2,[0,3,0],irs,conf);

We got very similar results.

And the correct one we got from the example impulse response in SSR is like

Would like to hear your comments. Many thanks!

At the moment a lot of code is duplicated to handle the 2D and 2.5D case. This will become worse if also consider the 3D case.

Hence, we should consider to implement the switching of the dimensionality via a configuration option. For example: conf.dimension

Try to use the bsxfun Matlab function.

The sound field is at the moment scaled at conf.xref. This might be a little bit misleading, because one would expect conf.xref to have only a contribution for 2.5D synthesis.
Maybe we should always scale to the center of the plot.

Add explanation how to use a BRIR data set of multiple loudspeakers measured in a room.
The best would be to use the 64 loudspeaker data set from Rostock.

Unlike for monochromatic cases, it seems that the Ambisonic order cannot be freely chosen for time-domain NFC-HOA driving function. It only works if the order is an integer multiple of the number of secondary sources, of vice versa. I started to implement for arbitrary integer orders. Would it be worth doing this?