sound_field_mono with non regular grids

[fietew]

Recently, I tried to evaluate a sound field on circles of different radii, which leads to non-regular grids in cartesian coordinates. At the moment sound_field_mono does only accept single values for X,Y,Z or a 2-element range vector for each dimension in order to build up a regular grid in cartesian coordinates. The functionality could be extended: If one of the the inputs, e.g. X, is a matrix of a certain dimensionality, then Y, Z have to either of the same size as X or scalars. Is it a good idea and/or do you see any problems for the implementation?

[hagenw]

Yes, the idea is not bad. I haven't thought of this before as I never needed it. Even the effort to implement this should be not to big. Maybe adjusting xyz_grid() can to the trick already.

[fietew]

Unfortunately, it's more complicated than expected. Implementing it in xyz_grid() is not that straightforward, since many cases have to handled and some functionalities of sound_field_mono have to be disabled when using non-regular grids. So far, I added the functionality to sound_field_mono directly. See non_regular_grid branch (6168ff3e470b5efae9d3ec7fc25cf0a4f80a09e9).

[hagenw]

I was playing around a little bit and am asking myself if it makes sense to allow for a complete arbritrary grid. For example, what I did is the following (you have to pull before):

conf = SFS_config_example;
x = sort(randi([-2000 2000],300,1)/1000);
y = sort(randi([-2000 2000],300,1)/1000);
X = repmat(x',[conf.resolution 1]);
Y = repmat(y',[conf.resolution 1]);
[P,x,y,z,x0] = sound_field_mono_wfs(X,Y',0,[0 -1 0],'pw',1000,conf);
P = norm_sound_field_at_xref(P,x(1,:),y(:,1),zeros(300,1),conf);
plot_sound_field(P,x(1,:)',y(:,1),zeros(300,1),x0,conf)

I created a random grid from vectors, using always the same x-axis for all y-values. If you don't do this, how you are able to plot the result? Maybe you could add a small example of what you are want to do with this.

[fietew]

It is up to the user, how he wants to plot the data afterwards. However, the scatter/scatter3 plotting function of matlab plots a marker for each [x,y,z] and is able to encode the soundfield (either real part or absolute value ) in the size and/or the color of the marker (I prefer the color). If you think, that a "default" plotting method for irregular grids should be present in plot_sound_field, then this could be an alternative. A short example:

conf = SFS_config_example;
x = sort(randi([-2000 2000],300,1)/1000);
y = sort(randi([-2000 2000],300,1)/1000);
X = repmat(x',[conf.resolution 1]);
Y = repmat(y',[conf.resolution 1]);
[P,x,y,z,x0] = sound_field_mono_wfs(X,Y',0,[0 -1 0],'pw',1000,conf);
colormap( 'jet' )
scatter(x(:),y(:),[],min(20, max(-20,real(P(:)))))

[hagenw]

That's a nice idea and works quite well. I will review the non regular grids code again the next days, add scatter as the default non-regular grid plotting method to plot_sound_field and merge it into master.

[hagenw]

I did some changes and moved most of the non-regular grid handling to xyz_grid(). With this it was easy to integrate it directly into the plot_sound_field() by using scatter() and also to integrate it into the time domain functions as you can see in this figure.

Some of the comments have still to be updated and I would like to tackle the problem of sound field normalization again before merging this into master.

[hagenw]

We should also add testing function for the non-regular grid in the validation folder.

[fietew]

I will do that, since I already have some scripts to non regular grids

[fietew]

Is it purpose, that the markers are non-solid, i.e. just the circle edges are colored. Maybe completely colored markers would be better?

[hagenw]

I haven't looked at how scatter() works, just copied the code from your example ;) So, if solid marker gave better results, just change it.

[fietew]

Regarding the normalisation I would suggest the following: The sound_field_imp and sound_field_mono function should return an additional normalization factor which can then be taken into account when a plotting routine is called. The computation of this normalization factor is different for time and frequency domain: For time you might just take the maximum abs value amongst the computed sound field values, for frequency you might explicitly compute the sound field at xref.

The main problem will be, that the factor must be provided to the plotting routine as an optional argument in the best case, since we might not want to modify all the calls of the plotting function in the whole toolbox. Any ideas?

[hagenw]

For mono and imp we have already two different normalization functions: norm_sound_field_at_xref() and norm_sound_field(). I guess the first one will not work for all non-regular grids and also returns an error if xref is outside of your listening area. In principal we could remove those normalization functions easily and move them to the plot function. There is only one problem for the mono case. Some time ago you proposed that we should allow the usage of xref positions outside of the listening area. This has the problem that we need to calculate what would be the amplitude at the position of xref, which is not possible inside the plotting function.

[hagenw]

OK, now I see the problem. The plotting function cannot know if we have a time-domain sound field or a monochromatic one. Mh, the idea with the additional parameter to the plotting function is not very nice in my opinion. On the other side the only alternative would be to have to different plotting routines, plot_sound_field_imp() and plot_sound_field_mono() or is there another alternative?

[hagenw]

One alternative could be to normalize the sound fields in both cases to 1 (as currently done in the time domain case). The problem, why this is not working at the moment in the monochromatic case is that the amplitude near the secondary sources is normally way to large. A solution might be, to exclude a predefined area around the secondary sources (let's say 10cm) and look only in the remaining field for the maximum value and scale that one to 1. I will try, if I'm able to implement something like this.

[fietew]

We should somehow split the configuration parameter conf.usenormalisation . I think, normalising the sound field to ensure constant sound pressure at a distinct point in space is a different purpose than normalising the sound field for visibility in plots.

[hagenw]

Yes, we could/should move it to conf.plot.usenormalisation. After our changes this will also be the only normalization applied at the moment. If we want to have also the "normalising the sound field to ensure constant sound pressure at a distinct point" we have to implement it in addition.

[fietew]

I agree with that option.

[hagenw]

OK, I'm working on it at the moment.

[hagenw]

Basically, everything discussed is now implemented. The only thing missing is the exclusion of the small area around the secondary sources for scaling, which means that monochromatic plots scale wrongly at the moment.

[hagenw]

I implemented another trick. norm_sound_field() now first tries to normalizes the center of the sound field to 1. If the absolute value found at its center is smaller than a given threshold (0.3 at the moment) it looked for the absolute maximum in the whole sound field and does the normalization accordingly, see https://github.com/sfstoolbox/sfs/blob/non_regular_grid/SFS_general/norm_sound_field.m#L58

The trick with the threshold detects if we have a monochromatic or a time-domain plot. If the value of 0.3 is wisely chosen and if it works as expected in all cases, has to be tested.

[fietew]

There are still some issues with plot_sound_field(), if the custom grid is a 1D-array. I guess the function can't distinguish between a regular sampled 1D vector (eg. x) and a non-regular 1D vector.

conf = SFS_config_example;
conf.plot.useplot = false;

alpha = 0:0.5:359;
r = 0.5;
[alpha, r] = ndgrid(alpha, r);
X = cosd(alpha).*r;
Y = sind(alpha).*r;
Z = 0;

sound_field_mono(X,Y,Z,[1,0,0,1,0,0,1], 'ps', 1, 1000,conf);

[P,x,y,z] = sound_field_mono(X,Y,Z,[1,0,0,1,0,0,1], 'ps', 1, 1000,conf);
plot_sound_field(P,x,y,z,[],conf);

A solution might to use the same input convention as for sound_field_mono:

• if you want a regular grid, then just define its extrema, e.g. [xmin, xmax], or a scalar, e.g. z, and the rest is done by the algorithm (imagesc supports the [xmin, xmax] input)

conf = SFS_config_example;
conf.plot.useplot = false;

X = [-2,2];
Y = [-2,2];
Z = 0;

[P,x,y,z] = sound_field_mono(X,Y,Z,[1,0,0,1,0,0,1], 'ps', 1, 1000,conf);

% compare
figure;
imagesc(x,y,real(P))
figure;
imagesc(X,Y,real(P))

• any dimension containing more then 2 elements will be interpreted as non-regular grids

[hagenw]

At the moment the non-regular grids (better custom grids) are handled in the following way:

• 2D sound field: you have to provide two of X,Y,Z as a matrix with the dimensions of the desired sound field P and the other one as a scalar
• 1D sound field: you have to provide one of X,Y,Z as a vector and the others as a scalar

This means your first example is not supported at the moment. The question is, should we support it? Otherwise we should of course add an error message.

For non-custom grids we always assume that you provide only [xmin, xmax] etc.

[fietew]

At the moment, your last sentence is not true imho. According to https://github.com/sfstoolbox/sfs/blob/de2aaf48655ea8bae214ad8f169f74976f5b711b/SFS_monochromatic/sound_field_mono.m#L163 , plot_sound_field is not called with X,Y,Z but with x,y,z (which are linear interpolated between the respective extrema, using xyz_grid and xyz_axes). Hence, x,y,z might contain 1D-arrays with more than two elements, which cannot be distuished from 1D-custom arrays.

[hagenw]

What I mean is the input to the sound field functions. The idea for the plotting function is different as here we assume you never provide your own x,y,z, but always the one calculated by a sound_field function. Is this to confusing and we should use X,Y,Z?

[fietew]

Independent of confusion, it hinders us to use 1D custom arrays as an input for the plotting function, since the plotting function interprets it differently. Maybe, I should come up with some drawing or scheme, to show which cases (i.e. regular or custom / 1D, 2D, 3D,... ) I would like to cover with the sound_field_mono() and plot_sound_field().

[hagenw]

Yes, this would be a good idea. Maybe, you could add it directly to the README. Later on, we then just clean it up a little bit and have the documentation ready as well.

[fietew]

After some brainstorming, I (re-)implemented some stuff. The following things have been done now:

• 3D scatter-plot ( scatter3) for regular and custom grids
• remove or simplify some of the helper functions for computing the grids, checking dimensions, etc.
• plot_sound_field() does no longer accept the axis x,y,z but the same inputs as for sound_field_mono()
• sound_field_mono() now returns the grid data (instead of the axes ) as the 2nd,3rd and 4th output value

TODO:

• sound_field_imp()
• documentation
• polishing of comments in *.m files

[hagenw]

Very nice, the code looks even cleaner now and we have a 3D plotting function :)

Octave has some problems with the scatter functions, I will have a look at it.

[hagenw]

One small problem remains, where we have to think about. In the old version you could do the following:

>> [P,x,y,z] = sound_field_mono_wfs([-2 2],[-2 2],0,[0 -1 0],'pw',1000,conf);
>> plot_sound_field(P,x,y,z,x0,conf)

This of course still works, but as x,y,z are now the actual grid values, plot_sound_field() always thinks that we have a custom grid and is using the scatter function.

[fietew]

Yeah that's a little drawback with this implementation. However, the inputs for plot_sound_field() and sound_field_mono() are now consistent. We could add a warning message to the plotting function in the next release (maybe something MATLAB-like "functionality has changed to blablabla. If you want to use the old functionality, then do this") and then drop the warning message in the release after that.

[hagenw]

I fixed the problem with the slow scatter plot in Octave by limiting the number of different P-values to 64 (the number of different colors that are used in the plot), see http://savannah.gnu.org/bugs/?40663

One funny problem remains, where I'm not sure what might be the reason. Run the following commands in Octave and Matlab and compare the result:

>> [P,x,y,z,x0] = sound_field_mono_wfs([-2 2],[-2 2],0,[0 -1 0],'pw',1000,conf);
>> plot_sound_field(P,x,y,z,x0,conf)

I got the following in Matlab (as it should look like)

And the following in Octave, where I don't know why red dominates

[hagenw]

I did further comparisons by running test_plot(). The results from Octave are always displayed to the left side, the ones from Matlab to the right side:

2D, custom grid

And now for the 3D case, where it seems that the algorithm in Octave works better. or say in the way we would prefer it:

3D

3D, custom grid

[fietew]

A vague guess about the behaviour of the 3D plots could be the following: Matlab plots the points in order to appearance as given in the inputs to scatter3. This means, that points which are farther from the view point (which should be occluded) might be plotted over nearer points. Gnuplot utilizes some kind of z-buffer (see https://en.wikipedia.org/wiki/Z-buffering ) as it uses OpenGL besides gnuplot ( see https://www.gnu.org/software/octave/doc/interpreter/Introduction-to-Plotting.html )

[fietew]

I'm getting an error, if I execute the current test_plot in Octave (v 3.8.1):

error: cellfun: dimensions mismatch
error: called from:
error:   /usr/share/octave/3.8.1/m/plot/appearance/axis.m at line 348, column 14
error:   /usr/share/octave/3.8.1/m/plot/appearance/axis.m at line 374, column 8
error:   /usr/share/octave/3.8.1/m/plot/appearance/axis.m at line 189, column 7
error:   /usr/share/octave/3.8.1/m/plot/appearance/axis.m at line 149, column 7
error:   /home/****/*****/sfstoolbox/SFS_plotting/plot_sound_field.m at line 213, column 5
error:   /home/****/*****/sfstoolbox/SFS_monochromatic/sound_field_mono.m at line 163, column 5
error:   /home/****/*****/sfstoolbox/SFS_monochromatic/sound_field_mono_wfs.m at line 99, column 29
error:   /home/****/*****/sfstoolbox/validation/test_plot.m at line 57, column 1

[hagenw]

Mh, strange. It happens during the call of axis image;, I don't know if we can do much to change this behavior. I'm running Octave 4.0.0 and get no error.

[fietew]

I updated to Octave 4.0.0 and that fixed the problem.

[hagenw]

As I like your new code, I would say we stay with the new behavior of plot_sound_field() regarding the handling of X,Y,Z and I added the following warning message which will appear whenever you handover a custom grid to plot_sound_field():

>> plot_sound_field(P,x,y,z,x0,conf)
Warning: You have provided a custom grid for the X,Y,Z input parameters of PLOT_SOUND_FIELD.
 If this was desired ignore this message, otherwise note the following:
 The behavior of PLOT_SOUND_FIELD has changed. Before, it worked in the following way:
   [P,x,y,z,x0] = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf);
   plot_sound_field(P,x,y,z,x0,conf)
 Now you should use instead:
   plot_sound_field(P,X,Y,Z,x0,conf)
 where X,Y,Z are the same input arguments as for the sound field function. 
> In plot_sound_field at 111

[hagenw]

I fixed also the 3D plotting problem in Matlab by adding

set(gcf,'renderer','opengl');

after the scatter3 command. It still works in Octave. Are there still computers around whithout OpenGL? I don't know what will happen then.

[hagenw]

I ran test_plot() under Octave 3.8.2 and it is working well. So, maybe we don't have too worry to much.

[hagenw]

So, in my opinion the code is now ready to be merged into master. The biggest change for the user is that there exist now the possibility to choose the normalisation method for potting:

% Normalisation method. Available methods are:
%   'auto'      - 'center' if center of sound field > 0.3, otherwise 'max'
%   'center'    - center of sound field == 1
%   'max'       - max of sound field == 1
conf.plot.normalisation = 'auto'; % string

I still need to clean up a little bit the function that generates the plots for the README and one bug remains, which I explain in the next comment.

[hagenw]

If you run test_plot() and compare the 1D y-axis plots regular grid vs. custom grid for the first round of time domain plots it looks like this (its figure 20 and 27):

Something seems to go wrong for the custom grid. Funnily this happens only for the y-axis plot, x- and z-axis are fine. By the way the same happens for the second round of time domain plots (they are the same, but shown in dB).

[fietew]

The virtual sound source is located at [0 -1 0]. As you are using randi to generate the custom grids, it is quite likely, that you compute the sound field at x2=-1, which yields Inf. The normalization of the sound field yields zeros for all the other x2. The regular grid from -2 to 2 with a resolution of 300 does not include x2=-1 and is therefore displayed correctly.

[hagenw]

That's true, I fixed it, now everything works :)