• WANTED: Happy members who like to discuss audio and other topics related to our interest. Desire to learn and share knowledge of science required. There are many reviews of audio hardware and expert members to help answer your questions. Click here to have your audio equipment measured for free!

DIY 3D Speaker Scanner - the Mathematics and Everything Else

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
One thing I noticed; there is a tranformation in the coordinates between the output grid of the simulated.sources and the final output grid. I norticed this when I wanted to subtract both grids to view the error.
I havent found yet where it is caused in the code
 
OP
NTK

NTK

Major Contributor
Forum Donor
Joined
Aug 11, 2019
Messages
2,656
Likes
5,819
Location
US East
One thing I noticed; there is a tranformation in the coordinates between the output grid of the simulated.sources and the final output grid. I norticed this when I wanted to subtract both grids to view the error.
I havent found yet where it is caused in the code
I used, as I wasn't playing much attention, different parameters to generate the output grid in the 'XZ' plane. In "step_3.m" line 62, I specified the angle range from -pi to pi. The input grid was generated using the default range of 0 to 2pi. I wonder if that is the problem. I didn't compare the inputs to the outputs to calculate the errors, so I have missed that.

Try changing the angle range in "step_3.m" line 62 to [0 2*pi] and see if it gives your desired results.

Sorry.
 

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
I used, as I wasn't playing much attention, different parameters to generate the output grid in the 'XZ' plane. In "step_3.m" line 62, I specified the angle range from -pi to pi. The input grid was generated using the default range of 0 to 2pi. I wonder if that is the problem. I didn't compare the inputs to the outputs to calculate the errors, so I have missed that.

Try changing the angle range in "step_3.m" line 62 to [0 2*pi] and see if it gives your desired results.

Sorry.
yup that did the trick, thanks!

a picture of a soucre with offset 0.6 0.6 0.4 with and without the correction:
the last plot is the difference plot between input and reconstruction (error plot)
1620208837885.png
 
Last edited:

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
Hereby I send the updated code I made, in the readme.txt I kept track off all made changes.
 

Attachments

  • MATLAB_release_2021_05_05.zip
    111 KB · Views: 104

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
NTK, could you please explain how I can get the response at an arbitruary (x y z) position if the spherical coefficients are calculated?
I have been fiddling with the code and the AKtools for a couple of hours, but I am not managing. For some reason I only manage on pre-calculated grids, but is should be possible once you hace the coefficients to get the response at any point right?
 
OP
NTK

NTK

Major Contributor
Forum Donor
Joined
Aug 11, 2019
Messages
2,656
Likes
5,819
Location
US East
NTK, could you please explain how I can get the response at an arbitruary (x y z) position if the spherical coefficients are calculated?
I have been fiddling with the code and the AKtools for a couple of hours, but I am not managing. For some reason I only manage on pre-calculated grids, but is should be possible once you hace the coefficients to get the response at any point right?
A very quick rundown of the fitting and reconstruction process:
  1. Use the coordinates of the measurements points (r, θ, φ), the frequency (ω), and the expansion order N to generate the measurement PSI matrix. The function to generate the PSI matrix is sph_PSI_mix().
  2. Input the measurement PSI matrix and the measured pressure amplitudes into the function lstsq_solve(), and it will return the expansion coefficients (and mean square fitting error).
  3. Use the coordinates of the reconstruction points (r, θ, φ), the same frequency (ω), the same expansion order N, and use sph_PSI_mix() to generate the reconstruction PSI matrix.
  4. The pressure amplitudes at the reconstruction points are obtained by the matrix multiplication of the reconstruction PSI matrix obtained in step 3 and the expansion coefficients in step 2.
I'll try to create an example script either later today or tomorrow to make it clearer.
 
OP
NTK

NTK

Major Contributor
Forum Donor
Joined
Aug 11, 2019
Messages
2,656
Likes
5,819
Location
US East
Here are 2 scripts to show the process of fitting and reconstruct sound pressures. I also added a function to calculate the dB_SPL, with the input parameter being the complex sound pressure in Pascals. You can run them with @Kessito's code (attachment in post #304). I used his concentric cylindrical measurement grid generator that include points in the top and bottom ends.

The first script 'fit_and_reconstruct_1.m' is the standard method. The second script 'fit_and_reconstruct_2_sfs.m' uses soundfield separation (SFS). Because of some really nice mathematical magic, the real differences between SFS and non-SFS is just one line of code (compare lines 139 in both scripts). Octave seems to have some problems with its numerical computations for the SFS script at the higher expansion orders. Haven't figured out why yet.

The scripts assumed we already know the optimal expansion order, which is one BIG assumption. When you play with the SFS script, you can see that the results are very sensitive to the expansion order. In the example I used 2000 Hz frequency, which is probably at or beyond the limit of SFS in that situation. You can see the best case error is still quite high. Klippel's NFS transitions time gated measurements for the higher frequencies.

I haven't started experimenting with different schemes to find the optimal expansion order. It will be my next big thing to attempt. I am also going to start experimenting with Professor Per Christian Hansen's regtools (https://www.mathworks.com/matlabcentral/fileexchange/52-regtools) for more robust ways of calculating the fitting coefficients (i.e. less sensitive to errors in the inputs).
 

Attachments

  • fit_and_reconstruct.zip
    4.6 KB · Views: 106

No. 5

Active Member
Joined
Nov 1, 2019
Messages
144
Likes
121
A huge thank you to @NTK for your continued work on this (and @Kessito for testing it out)!

For those of us who's ability to use Octave doesn't extend beyond basic math, could there be a tutorial/directions at some point on how to use those scripts?
 
OP
NTK

NTK

Major Contributor
Forum Donor
Joined
Aug 11, 2019
Messages
2,656
Likes
5,819
Location
US East
A huge thank you to @NTK for your continued work on this (and @Kessito for testing it out)!

For those of us who's ability to use Octave doesn't extend beyond basic math, could there be a tutorial/directions at some point on how to use those scripts?
Certainly! Thanks for your continuous interest and support :) I will write up a little tutorial explaining the 2 scripts I posted in post #307. Should be able to post it in a few days.
 

Tokyo_John

Active Member
Forum Donor
Joined
Mar 6, 2021
Messages
214
Likes
289
This is a nice bit of work, thanks for sharing it.

It is similar in nature to what some of us do in geophysics, using seismic records to infer the directionality and rupture pattern of an earthquake over a fault zone. In that problem the measurements are not always ideally placed, as we only have a certain number and often restricted locations of seismographs. Additionally, unlike clean air, the propagation medium is strongly heterogeneous and the waves bounce and refract and scatter along the way from source to receiver. Finally, there are both sound waves and elastic shear waves, which offers more information but a more complex recording (3-axes of movement, rather than a scalar pressure field).

The basic problem of determining the magnitude and fault orientation for an earthquake, called the "moment tensor solution" (which has 6 independent components for a "double-couple" moment in 3D, 3 unknowns for location, and 1 unknown for event time) is carried out in a way that is very similar to the approach you describe here, involving an analytical expansion of the solution in terms of coupled spherical harmonics and orthogonal basis functions that account for radial variations in seismic wave velocities inside the Earth. Combining data from hundreds of seismograph stations results in a non-square sensitivity matrix...typically in seismology we make up for non-ideal coverage of seismic stations by adding as many recordings as possible. The inversion is performed using singular value decomposition (SVD), which returns the least squares solution in terms of generalized eigenvalues and eigenvectors, from which one can infer the dominant modes of the solution as well as the covariances. It is a beautiful technique and I highly recommend SVD. The general theory is described in a book titled "Theoretical Global Seismology" by Dahlen and Tromp.

For inferring rupture patterns and heterogeneous energy radiation along a fault zone (asperities, etc.), other techniques such as "beam-forming methods" are used. This is another approach you could look into, if you want to fall further down the rabbit hole.

In any case, I'd be happy to give you feedback or point you to more detailed references if you're interested.
 
OP
NTK

NTK

Major Contributor
Forum Donor
Joined
Aug 11, 2019
Messages
2,656
Likes
5,819
Location
US East
@Tokyo_John Thank you very much for the feedback and offer.

Finding more robust methods to compute the fitting coefficients is one of the bigger items in my to-do list (the other one is the determination of the optimal expansion order). The numerical problem can become unstable when measurements are noisy and are taken in reflective rooms. With reflections, the sound sources are all spread out and the "acoustic center" is no longer well defined. The "previously optimal" (i.e. if there is no reflection) location of the expansion center will be a lot less optimal. The sound radiation pattern is also going to be much more complicated. All these will require us to use much higher order expansion functions, and we will need to solve for many more fitting coefficients.

I am reading up on regularized least squares. Earl Williams, one of the primary inventors and leading expert on nearfield acoustical holography, has a paper on the topic. I still have a lot to learn before I can fully understand his paper. I do think I have a decent handle on SVD, and am planning to start experimenting with Professor Hansen's MATLAB regularization tools (regtools) package.
https://www.researchgate.net/public..._Journal_of_the_Acoustical_Society_of_America

However, the problems folks like Earl Williams are trying to solve are a lot worse than what we have here. They are trying to figure out the vibration modes of the vibrating surfaces using mic measurements of the sound waves they produce. So hopefully we can get by with using just the basic techniques.
 

No. 5

Active Member
Joined
Nov 1, 2019
Messages
144
Likes
121
Certainly! Thanks for your continuous interest and support :) I will write up a little tutorial explaining the 2 scripts I posted in post #307. Should be able to post it in a few days.
Thank you for your continued time and effort in development!

This is something that I've been interested in seeing for the past four years or so, and my own limited mathematics experience has certainly hampered the progress that I would have been able to make myself, so I find it very exciting to see someone else who's capable and interested in doing something like this.

Thanks again. :)
 
  • Like
Reactions: NTK

Tokyo_John

Active Member
Forum Donor
Joined
Mar 6, 2021
Messages
214
Likes
289
Finding more robust methods to compute the fitting coefficients is one of the bigger items in my to-do list (the other one is the determination of the optimal expansion order).

This kind of information is, in principle, obtainable from the SVD, which is unconditionally stable. One can perform the magic of replacing 1/0 with 0 to get rid of trouble.

The numerical problem can become unstable when measurements are noisy and are taken in reflective rooms. With reflections, the sound sources are all spread out and the "acoustic center" is no longer well defined. The "previously optimal" (i.e. if there is no reflection) location of the expansion center will be a lot less optimal. The sound radiation pattern is also going to be much more complicated. All these will require us to use much higher order expansion functions, and we will need to solve for many more fitting coefficients.

Yes, this is definitely a problem, measuring the room instead of the source, ideally you want to be able to separate those, but the only way to know for sure where the signal is coming from is to use phased array approaches. It is somewhat more complicated, but something you might consider. Array approaches such as beamforming are commonly used by anti-missile systems, sonar signal processing, acoustical science, seismology, etc..

This group shares their code for beamforming using microphone arrays...
https://github.com/MiguelBlancoGalindo/MicArrayBeamforming

I am reading up on regularized least squares. Earl Williams, one of the primary inventors and leading expert on nearfield acoustical holography, has a paper on the topic. I still have a lot to learn before I can fully understand his paper. I do think I have a decent handle on SVD, and am planning to start experimenting with Professor Hansen's MATLAB regularization tools (regtools) package.

Sounds good, SVD offers a way to use norm regularization (or no regularization if one recomposes the solution and tunes the expansion just so).
 
Last edited:

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
I think that beamforming might not be the best option for acoustic center finding in this case, because we know what the input stimulis and min/max time of flight is.
I think that the multilaterarion option I am testing with is promising, I only need to get unwrapped phase for the frequencies of interest.
 

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
source separation with offset destection working pretty well:
Cardioid with 2 disturbing sources outside measurment cylinder: max 0.5dB deviation @ 250 Hz with 9th order

Next step I will use aktools room to simulate a " real" reflective small room around the setup

biggest challenge now is getting high enough order fit without "exploding" least squares I think
Original:
1620648992649.png

1620649153981.png


with disturbances and reconstructed with sfs:
1620649101014.png


1620649053460.png
 

Attachments

  • 1620648842268.png
    1620648842268.png
    228.9 KB · Views: 89
  • 1620648924558.png
    1620648924558.png
    81.1 KB · Views: 84

Ulf

Member
Joined
May 11, 2021
Messages
9
Likes
12
(sorry if this was already discussed) You can do better integration by using a Gaussian type grid, i.e. Lebedev grids for the sphere. Using carefully placed points (and weights) these are designed to integrate all spherical harmonics up to a given order exactly. The weights correspond to a different weighting of your least square system, which can also give more robust fitting results.
 

Kessito

Member
Joined
May 18, 2020
Messages
54
Likes
99
(sorry if this was already discussed) You can do better integration by using a Gaussian type grid, i.e. Lebedev grids for the sphere. Using carefully placed points (and weights) these are designed to integrate all spherical harmonics up to a given order exactly. The weights correspond to a different weighting of your least square system, which can also give more robust fitting results.

Yes I am (we are) aware, but I want to have the flexibillty of a cylindrical grid; this allows you to make a much easier array with the possibillity of measuring large cabs (I want to be able to measure cabinets of 1.2x1.2x1.2m. See post 296
 

Ulf

Member
Joined
May 11, 2021
Messages
9
Likes
12
Yes I am (we are) aware, but I want to have the flexibillty of a cylindrical grid; this allows you to make a much easier array with the possibillity of measuring large cabs (I want to be able to measure cabinets of 1.2x1.2x1.2m. See post 296
You can still use Gauss-Legendre in the z direction (maybe you are?). Not so much to save measurement time, but to make the solution more accurate. When you regularize in your least square solution you want the vector norm to correspond to the physical energy of the wave field, and the grid weights (can) give you that.
 
Last edited:

AdamG

Proving your point makes it “Science”.
Moderator
Forum Donor
Joined
Jan 3, 2021
Messages
4,636
Likes
14,918
Location
Reality
You can still use Gauss-Legendre in the z direction (maybe you are?). Not so much to save measurement time, but to make the solution more accurate. When you regularize in your least square solution you want the vector norm to correspond to the physical energy of the wave field, and the grid weights (can) give you that.
Welcome Aboard @Ulf.
 
OP
NTK

NTK

Major Contributor
Forum Donor
Joined
Aug 11, 2019
Messages
2,656
Likes
5,819
Location
US East
(sorry if this was already discussed) You can do better integration by using a Gaussian type grid, i.e. Lebedev grids for the sphere. Using carefully placed points (and weights) these are designed to integrate all spherical harmonics up to a given order exactly. The weights correspond to a different weighting of your least square system, which can also give more robust fitting results.
Thank you very much for the info. In my research on the implementation of NAH (nearfield acoustical holography), I haven't come across any work on finding the optimal measurement grids. Doesn't necessarily mean that there aren't any. The complication is that this will require a very elaborate robotic measurement system.

I downloaded the MATLAB scripts from the link below, and should be able to run tests with them and compare results with equidistance grids. But I probably won't be able to start for another couple of weeks.
https://people.sc.fsu.edu/~jburkardt/m_src/sphere_lebedev_rule/sphere_lebedev_rule.html
 
Top Bottom