# DIY 3D Speaker Scanner - the Mathematics and Everything Else

OP

#### NTK

##### Addicted to Fun and Learning
Forum Donor
I got my attention diverted to checking out the AKTools. I don't have the signal processing or the statistics toolboxes, and was trying to port the functions from Octave as I go trying to get some of the demos to work. Hasn't gotten very far.

I am going to go back to work on the MATLAB fitting/reconstruction routines. I have the code for the spherical Bessel and Hankel function completed a little while ago, so there shouldn't be any show stoppers for me. I will also try using these same functions from AKTools as I may need some of its other functions in the future.

#### Kessito

##### Member

https://nl.mathworks.com/matlabcent...ic-transform-gaunt-coefficients-and-rotations

He has already done the least squares fitting for arbituary measurement points

Here is the link to his PhD thesis for which he did the work:

https://aaltodoc.aalto.fi/handle/123456789/22499

https://aaltodoc.aalto.fi/bitstream.../isbn9789526070377.pdf?sequence=1&isAllowed=y
good working code also:

function [F_N, Y_N] = leastSquaresSHT(N, F, dirs, basisType, weights)
%LEASTSQUARESSHT Spherical harmonic transform of F using least-squares
%
% N: maximum order of transform
% F: the spherical function evaluated at K directions 'dirs'
% if F is a KxM matrix then the transform is applied to each column
% of F, and returns the coefficients at each column of F_N
% respectively
% dirs: [azimuth1 inclination1; ...; azimuthK inclinationK] angles in
% rads for each evaluation point, where inclination is the polar
% angle from zenith: inclination = pi/2-elevation
% basisType: 'complex' or 'real' spherical harmonics
% weights: vector of K weights for each direction, for weighted
% least-squares
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Archontis Politis, 10/10/2013, updated 7/2/2015
% [email protected]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the harmonic coefficients
Y_N = getSH(N, dirs, basisType);
Npoints = size(dirs,1);

% non-weighted case for uniform arrangements
if nargin<5 || ~exist('weights','var') || isempty(weights)

% perform transform in the least squares sense
F_N = pinv(Y_N)*F;
else

% perform weighted least-squares transform
F_N = (Y_N'*diag(weights)*Y_N) \ (Y_N'*diag(weights) * F);
end

end

Last edited:

#### Kessito

##### Member
And last but not least:
a 1 deg 3D anechoic measurment set of set of a 3-way speaker:

https://depositonce.tu-berlin.de/bitstream/11303/7506.2/25/Documentation.pdf

coming from this treasure vault:

https://www.sofaconventions.org/mediawiki/index.php/Files

2.2 Source and receiver descriptions Three different sources and three receiver types were used for creating the BRAS database: A 3-way dodecahedron loudspeaker and a Genelec 8020c 2-way near-field monitor were used for measuring single channel impulse responses, and QSC K8 2-way PA monitors were used for measuring binaural impulse responses. A G.R.A.S. AF40 1/2” free field measurement microphone was used for the acquisition of single channel impulse responses in the free-field environments, while a 1/2” B&K 4134 diffuse-field measurement microphone was used in reverberant scenarios. The FABIAN head and torso simulator was used for measuring binaural impulse responses. Directivities by means of impulse responses and/or third octave spectra are stored in comma-separated value (CSV) and/or MATLAB files inside the folder 2 Directivities. They are provided on a Gauss-like spherical sampling grid with an angular resolution of 1◦ along azimuth, and elevation, with a total of 64,442 sampling points. Each line in the CSV-files holds the data for one sampling point of the grid as specified in Figure 1 (c, d). The structure of the provided MAT-files is shown in Figure 1 (e, f). Note that different coordinate conventions are used for loudspeaker directivities (Figure 1a), and head-related impulse responses (HRIRs) (Figure 1b). Moreover, the coordinate convention of the directivity files is independent from the absolute coordinate system of the scene geometries given in the SketchUp files. The directivities of the dodecahedron loudspeaker, and the two measurement microphones were assumed to be omnidirectional.

open science boys; [email protected]#k yeah!

Of course even better for testing is if we could get a raw dataset from Amir, to see if we can reeplicate the excact results!

Last edited:

#### mwmkravchenko

##### Member
Of course even better for testing is if we could get a raw dataset from Amir, to see if we can reeplicate the excact results!
That would require the same loudspeaker at the very least.

Now I am wondering how dependent on the mircrophone this will also become. I have 6mm mics up to 25mm mics. Obviously there will be a bandwidth difference between 6, 12 and 25mm microphones. But there is also a directional pickup pattern difference.

Just thinking out loud. Most affordable mics are 6mm capsules.

#### Kessito

##### Member
That would require the same loudspeaker at the very least.

Now I am wondering how dependent on the mircrophone this will also become. I have 6mm mics up to 25mm mics. Obviously there will be a bandwidth difference between 6, 12 and 25mm microphones. But there is also a directional pickup pattern difference.

Just thinking out loud. Most affordable mics are 6mm capsules.
No, most affordable are the mems, which you know how tiny they are

#### mwmkravchenko

##### Member
No, most affordable are the mems, which you know how tiny they are
I agree on the affordability.

And if that is what gets used, Which again I think is a great concept and I am thoroughly behind. There has to be a difference check created if there is any between the single mike multi point test versus the multi mems mic test. I have no doubt there can be a reconciliation between the two testing methods. Just to make the have valid inter-relatable testing the has to be comparisons to determine validity.

Anyone that has spent anytime behind a microphone knows that getting repeatable results is the toughest part. Close to repeatable is easy. Exactly repeatable hardly ever happens unless you have a very precise setup. Both for the microphone and the loudspeaker under test.

#### Kessito

##### Member
I agree on the affordability.

And if that is what gets used, Which again I think is a great concept and I am thoroughly behind. There has to be a difference check created if there is any between the single mike multi point test versus the multi mems mic test. I have no doubt there can be a reconciliation between the two testing methods. Just to make the have valid inter-relatable testing the has to be comparisons to determine validity.

Anyone that has spent anytime behind a microphone knows that getting repeatable results is the toughest part. Close to repeatable is easy. Exactly repeatable hardly ever happens unless you have a very precise setup. Both for the microphone and the loudspeaker under test.
Very true, calibration is mandatory. And probably some form of quick check calibration before measurement.
1-calibrate all mics by comparison method
2-do measurement with mics in situ with a known source ie a 1inch dome ( safe this source and dont touch it for other project) this is the reference
3- before measurement, do a check measurement with the known source and tje mics in the same position to check if there is "no" deviation from the reference situation ( no being ie <0.25dB)

#### mwmkravchenko

##### Member
3- before measurement, do a check measurement with the known source and tje mics in the same position to check if there is "no" deviation from the reference situation ( no being ie <0.25dB)
Be dancing with 1 to 1.5 over a wide measurement range. There are so many variables. If you can consistently get 1db or better you are doing a stellar job. If you can repeat measurements at 0.25db accuracy you are a measurement master. All hail the master micer

Joking aside, there are many variables in a setup and a test run and if you to easily get the same results that means that there is some smoothing going on that is glossing over the differences. As an example. Changes in air temperature. Forced air heating or cooling. even the fortuitous environmental noise that is just right not to be band rejected by the software. Although if I use discrete sine wave testing with my woofer tester pro it is pretty immune to noise. That testing is high resolution and quite useful but takes a loooonnng time to do. Oh and is rather nasty to listen to. But gets you the most realistic low frequency testing. Most short signals do not exercise the woofers enough to really get a real life measurement of what your woofers are actually doing if you use a sweep or a chirp or any of the noise flavours.

OP

#### NTK

##### Addicted to Fun and Learning
Forum Donor

https://nl.mathworks.com/matlabcent...ic-transform-gaunt-coefficients-and-rotations

He has already done the least squares fitting for arbituary measurement points

Here is the link to his PhD thesis for which he did the work:

https://aaltodoc.aalto.fi/handle/123456789/22499

https://aaltodoc.aalto.fi/bitstream.../isbn9789526070377.pdf?sequence=1&isAllowed=y
good working code also:
View attachment 119663
function [F_N, Y_N] = leastSquaresSHT(N, F, dirs, basisType, weights)
%LEASTSQUARESSHT Spherical harmonic transform of F using least-squares
%
% N: maximum order of transform
% F: the spherical function evaluated at K directions 'dirs'
% if F is a KxM matrix then the transform is applied to each column
% of F, and returns the coefficients at each column of F_N
% respectively
% dirs: [azimuth1 inclination1; ...; azimuthK inclinationK] angles in
% rads for each evaluation point, where inclination is the polar
% angle from zenith: inclination = pi/2-elevation
% basisType: 'complex' or 'real' spherical harmonics
% weights: vector of K weights for each direction, for weighted
% least-squares
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Archontis Politis, 10/10/2013, updated 7/2/2015
% [email protected]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the harmonic coefficients
Y_N = getSH(N, dirs, basisType);
Npoints = size(dirs,1);

% non-weighted case for uniform arrangements
if nargin<5 || ~exist('weights','var') || isempty(weights)

% perform transform in the least squares sense
F_N = pinv(Y_N)*F;
else

% perform weighted least-squares transform
F_N = (Y_N'*diag(weights)*Y_N) \ (Y_N'*diag(weights) * F);
end

end
You guys are catching on. It is really quite simple. The gist is to solve this equation (from the Part 4 PDF in post #1).

The left hand column in the matrix form is the measurement data at the measurement locations, given by (r, θ, ϕ). The majority of work is to assemble the matrix in the right hand side of the equation. Each of its elements is calculated by multiplying the spherical harmonics Y_mn and the spherical Hankel h_n(1) or spherical Bessel j_n functions, with the coordinates (r, θ, ϕ) as the parameters. There are 2(N + 1)^2 columns in the matrix, where N is the highest expansion order. k is the acoustic wavenumber (k = ω/c).

Then, use the MATLAB backslash operator (\) to solve for the coefficients C_mn and D_mn. The C_mn terms describe the outgoing wave. So, discard the D_mn terms and we have field separation. Repeat for all frequencies of interest, and that's basically it!

To evaluate when the expansion order N is sufficient, I am thinking of using convergence. Increase N until reconstructions at, for example, the 70 spinorama locations (or a subset of them if we are not running full sphere measurements) no longer change significantly.

So that really is all that is to it (well, at least to my understanding)

#### mwmkravchenko

##### Member
The last time I worked with any Calculus was 34 years ago. And I was good enough to get it done barely. This is beyond what I have taken for mathematics. But. I can learn. The equations are there. All I need is a understanding of what each abbreviation stands for. Oh and to look up Hankel and Sperical Bessel function. That I vaguely remember a bit about. I think I understand the underlying concepts. But I have never taken the time to understand exactly what you are doing. But I will read your 5 papers and get my little brain sore from thinking

#### napilopez

##### Major Contributor
Forum Donor
Hello all of you smart people. I've skimmed through the thread, so pardon me if I've missed this. Nor do I mean to derail the conversation much.

But as someone who has been creating Spinoramas the 'old' DIY way for a year or two, I was wondering if there was a way to extrapolate and simplify some of the findings in this thread to just a single on-axis measurement?

My thought is, I've already had much success approximating anechoic results by splicing time-windowed measurements with nearfield bass measurements. In terms of directivity, the results are nearly identical. For example from my recent JBL HDI-1600 measurements vs the NFS:

After a couple of tweaks, the on-axis is even closer:

I'm mostly happy with these results, as they are within ±1dB for most of the range. However the biggest problem with this method is the low resolution from roughly 200Hz-1kHz, where I often miss the narrow resonances and cannot tell if bumps are from my measurement rig or from the speaker.

Ground plane is not a viable method where I live, so I thought I'd maybe I could make the Klippel method work for just the on-axis sweep and keep doing things the old way for directivity. This I think would be a very useful compromise for my measurement goals.

So what would I need to do to remove reflections from just one angle as opposed to a whole sphere? Would I be able to measure a speaker from two distances and use that information to 'correct' the on-axis sweep?

If this derails the thread too much, I'm happy to start another thread. But I figured since the smart people were already gathered here... The math here is also way over my head, but I'm a good at teaching myself stuff.

Last edited:
OP

Forum Donor

#### No. 5

##### Member
That's fantastic! Thank you for highlighting it. I wish I had seen that sooner.

I completely agree that a full scan is the gold standard, and what I'd ultimately like to see a DIY solution to, but if you're already measuring outside, being able to take out the ground bounce is a huge advantage.

OP

#### NTK

##### Addicted to Fun and Learning
Forum Donor
That's fantastic! Thank you for highlighting it. I wish I had seen that sooner.

I completely agree that a full scan is the gold standard, and what I'd ultimately like to see a DIY solution to, but if you're already measuring outside, being able to take out the ground bounce is a huge advantage.
I saw your post #217. I thought one answer to your question is beamforming microphone arrays. At the time I was thinking about the more typical planar microphone arrays, which is probably not very practical for speaker measurement DIYers.

It wasn't until the discussion of this field separation method came up that this may be a more viable beamforming method. This configuration of beamformer is called endfire. Here is an link:
https://www.vocal.com/beamforming-2/endfire-microphone-array-beamforming-with-small-beamwidth/

Figure 2 shows, with 16 mics, we can get a highly directional response (20 degs beamwidth). The setup as shown in the picture in post #272 can simulate an endfire beamforming array. If we can take a series of measurements, say 50, exactly the same way, it is equivalent to having a 50 mic endfire array. Theoretically, as long as the front wall (wall behind the speaker) is far away, the other bounces can be filtered out and can give us a near anechoic response.

#### Few

##### Member
Do those spatial filtering diagrams specify a wavelength? I couldn’t find it. I would expect an array of the type in the photo would be super-directional at high frequencies but not directional when the wavelength is a few meters long. Maybe there’s a straightforward way of dealing with that, but it seems like a significant complication for a measurement system. Or am I missing something? Is the intent just to use multiple fixed microphones rather than moving a single one?

Few

OP

#### NTK

##### Addicted to Fun and Learning
Forum Donor
Do those spatial filtering diagrams specify a wavelength? I couldn’t find it. I would expect an array of the type in the photo would be super-directional at high frequencies but not directional when the wavelength is a few meters long. Maybe there’s a straightforward way of dealing with that, but it seems like a significant complication for a measurement system. Or am I missing something? Is the intent just to use multiple fixed microphones rather than moving a single one?

Few
I am not at all familiar with the beamforming math, so I'll have to dig into it. But, yes, it is frequency dependent.

This is an unrelated discussion to the reverse engineering of the Klippel NFS 3D scanner. I think I should take this discussion to a separate thread.

#### mwmkravchenko

##### Member
Do those spatial filtering diagrams specify a wavelength? I couldn’t find it. I would expect an array of the type in the photo would be super-directional at high frequencies but not directional when the wavelength is a few meters long. Maybe there’s a straightforward way of dealing with that, but it seems like a significant complication for a measurement system. Or am I missing something? Is the intent just to use multiple fixed microphones rather than moving a single one?

Few
I have to agree with the interesting part of the idea. But it is not the Klippel method as NTK posted. You are as bad at derailing discussions as I am!

#### Few

##### Member
Sorry, I really didn’t mean to derail. I was trying to judge the direct applicability of the microphone arrays to the main thread. If multiple mics reduce the number of “positions” that have to be measured, that’s clearly advantageous. I wasn’t sure if there was also an opportunity to apply beamforming ideas. I didn’t want to miss something!

OP

#### NTK

##### Addicted to Fun and Learning
Forum Donor
Here is the preliminary version of my MATLAB code. It isn't complete, and below is a brief description of what it does:
• I haven't finished the MATLAB code to simulate the measurements, and I used my Python code to generate the simulated measurements. There are 3 sets of simulated measurements, in the directors 'data\conc', 'data\cyl', and 'data\planar'. They are simulated double layer measurements but for one sided (i.e. only in the front hemisphere instead of full sphere). This is just to test/demonstrate how well/poorly the method works for one-sided measurements.
• Run the 'main.m' script for the full simulation.
• Each set of the simulated data is organized into 3 comma separated value (CSV) files, 'coordinates.csv', 'frequencies.csv', and 'measurements.csv'.
• 'coordinates.csv' are the coordinates of the measurement points: one point for each line, and each line with the x-, y-, and z- coordinates.
• 'frequencies.csv': one measurement frequency per line.
• 'measurements.cvs': These are the simulated pressure measurements (complex amplitudes) in complex numbers. The lines correspond to the measurement points, and columns correspond to the frequencies. Thus, the number of lines equals to the number of measurement points, and the number of columns equals to the number of frequencies.
• Each provided data set has about 3300 measurement points at 9 frequencies. A plot generated with the Python code for 2000 Hz is included to help evaluating the reconstruction accuracy.
• At the end of the simulation run, the variable 'p2034' will contain the reconstruction results at the CTA2034 grid points (at 2 m distance, not corrected to 1 m). The columns are the results for each expansion order, i.e. from 0 to Nmax. The results are not saved.
There are still much missing.
• Only mixed field reconstruction is provided. Field separation is not yet implemented but it is straightforward to do.
• Algorithm to decide the optimal reconstruction expansion order has not been developed. The combination of least squares residuals and convergence metric seems like sufficient information to make the decision. Still need to decide on a method to decide the optimal point.
• For more thorough testing, will need MATLAB code to simulate measurements, and more realistic simulated measurements.
• Near term pllan is to port the code for simulating measurements to MATLAB, implement field separation, run more numerical experiments/tests, and decide on the method to find the optimal expansion order.
Please let me know if you have any question or comments or need explanation/clarification. Sorry for taking so long.

#### Attachments

• 698.1 KB Views: 13

#### mwmkravchenko

##### Member
You have been a little busy! I'll have to get my little brain around using Matlab to take a look.