Hi,
I would like to drag and drop a folder with one shot samples into my plugin, and split/filter the contents of the folder and subfolders into several menus/groups by name.
Here’s an example: samples with “Clap, CP, CL” in their name will be added to the appropriate menu, and for example “Snare, SN, SD” to another, and so on…
Can anyone help me with this?
Any math heads in the house?
Basically trying to find the intersection of two hyperbolas given a pair of quadratic equations (but in Max).
So with two equations:
0.000171(x)^2 - 0.00023(y - 101)^2 = 1
0.00054(y)^2 - 0.00012(x - 101)^2 = 1
I want to arrive at all the possible intersections. Which in this case are:
(-76.8718, 94.2476)
(101.819, 43.0349)
(208.908, -66.6292)
(-406.089, -242.887)
And looks something like this:
Or here’s an embedded one:
So I can plot things to find the answers, but I need to be able to compute this on the fly, and do so in Max.
///////////////////////////////////
context:
using distance time of arrival (TDOA) I’m trying to compute the absolute positino on a strike on a drum with 4 sensors mounted on it.
the results are looking very promising so far, but I’m struggling with moving the computation step into Max:
(the red x is the predicted point, as the mean of the intersections for all the equations in the top left)
a_1x^2 - b_1y^2 = 1
y = \pm \sqrt((a_1x^2 - 1)/b_1)
then take each branch of that solution, plug into the other quadratic equation. for example with the positive branch:
a_2x^2 - b_2y^2 = 1
a_2x^2 - b_2((a_1x^2 - 1)/b_1) = 1
and after some simplifying (please check my arithmetic here, i’m still having my coffee)
x = \pm\sqrt(1 - (b_2/b_1) - a_2 + (b_2a_1/b_1))
Oof.
Ok, pardon some of the (potentially) dumb-ish question(s) as I am definitely punching above my weight in terms of applying this kind of math.
So this much makes sense as this is just a restructuring of the equation in the form that I have it above. I guess in my case there will be a few more twists and turns because things are rotated/translated (e.g. (y - 101)^2 in the first equation).
So in my specific case there will be a bit more algebra moving bits around, but the idea is to solve for y, straight up, with a ± to give two branches.
For the next step here:
Am I literally plugging (a_1x^2 - 1)/b_1) into there, or do I plug in both versions of y from above? (i.e. two floats (±), or populating that bit of equation with the coefficients a_1/b_1).
And then this bit:
So once producing two branches for y as above, this would be a condensed restructuring to solve for both versions of x in one line?
As above, this would need a bit more algebra go jiggle things about for the rotation/translation.
But in effect I would be producing two functions (each with a ±) that would get 4 coefficients plopped in the top of(?).
Lastly, I guess I’ll find this out when I get there, but sometimes there will be no intersection between the hyperbolas. I don’t know if this will produce a nan/INF or does there need to be code checking for things like that so a loop doesn’t hang or something?
oop sorry i forgot about the constant coordinate shift, i don’t think it makes things any more complicated - well i guess i was wrong there.
maybe i can try and update and actually write a solver to plug your example into.
sorry i know this is a max thread but i don’t have any particular thoughts about wirting solvers in max, seems a potentially clunky tool for it.
but the short answer is: you end up with 4 pairs of solution coordinates. some solutions ask you to take square root of a negative number. those intersections don’t lie on the real number plane. you can avoid nan from the computing environment by checking the argument to sqrt.
That would be super amazing/useful.
Obviously implementing it in Max will be unpleasant but thankfully in my circumstance the rotation/translation stuff is constant, as in, I’m projecting the hyperbolas into this specific orientation always, as it’s based on the physical position of the sensors on the drum. So I end up with four rotations/translations where one strike on the drum produces four equations/hyperbolas:
0.000171(x)^2 - 0.00023(y - 101)^2 = 1
0.00054(y)^2 - 0.00012(x - 101)^2 = 1
0.000592(x)^2 - 0.000117(y + 101)^2 = 1
0.000163(y)^2 - 0.000247(x + 101)^2 = 1
So I’ll solve for each pair of hyperbolas, which would theoretically give 4 intersections per pair (barring ones where they don’t intersect), but the signage (i.e. which quadrant the initial time offset was measured) would reduce each set of intersections to a single (potentially) valid point. Yielding up to 6 (potentially) valid intersections.
(you can make it least clunky in max by using ‘gen’(without the tilde)… i might be able to help there when the solver is ready… tho i’m not sure i understand what i’m reading yet, still, looks amazing 
)
That’d be super useful. You popped into my mind as I know you’re quite on the maths stuff, but I didn’t have a concrete enough way to get to even formulating what the problem is…
A buddy of mine made a version that does it using scipy in Python, so trying to figure out what’s underneath fsolve (may be minpack ?) to try and porting that over to Max/gen/codebox.
I’ve updated it to take a list of coefficients as well as computing all 6 intersections (instead of 4).
solve-scipy.py (2.2 KB)
thank you! (unfortunately, i’m only good at finding some inefficient workaround, haha… ‘maths’ is something else entirely… 
)
the scipy stuff makes sense overall, but then i look at ‘fsolve’ and realize it finds the root of functions:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html
that made my head-explode 
however! i’ll keep looking at stuff, as just your mention of minpack helped me find the exact part of the source to look for:
i have no clue, though, which one is relevant, ‘hybrd’
https://www.math.utah.edu/software/minpack/minpack/hybrd.html
‘hybrd1’
https://www.math.utah.edu/software/minpack/minpack/hybrd1.html
‘hybrj’
https://www.math.utah.edu/software/minpack/minpack/hybrj.html
or ‘hybrj1’
https://www.math.utah.edu/software/minpack/minpack/hybrj1.html
(and then it still might take me years to understand them 
)
and not that i understand this, either, but they mentioned that they use modified ‘Powell’s dog leg method’:
i wish i understood all that 
(but hoping adding all these links might help others, too)
those things you’re referring to are numerical solvers. i’m suggesting that with 2 quadratics and 2 unknowns, an algebraic / non-iterative solution should be possible. and it would be more straightforward to implement (once the algebra is worked out for your constraints.)
without the offsets it’s quite simple to solve with direct substitution, as shown above.
with the offsets, i tried for 5 minutes this morning but that specific direct substitution approach doesn’t quite work (or at least, it became too hairy for me rather quickly.) you give the problem this way:
\textbf{(1)}. ax^2 - b(y-n)^2 = 1
\textbf{(2)}. cy^2 - d(x-n)^2 = 1
where a, b, c, d, n are constant.
bearing in mind that all coefficients can be negative, you can expand all the terms, group them on one side, and put this in a more standard form:
\textbf{(1)}. a_0 + a_1x + a_2y + a_3xy + a_4x^2 + a_5y^2 = 0
\textbf{(2)}. b_0 + b_1x + b_2y + b_3xy + b_4x^2 + b_5y^2 = 0
where some of the general coefficients (say: a_1, a_3, b_2, b_3) are zero.
[ed] more specifically, given your original constants a, b, c, d, n, in the more general form i ended up with:
a_0 = -bn^2 - 1
a_1 = 0
a_2 = 2bn
a_3 = 0
a_4 = a
a_5 = -b
b_0 = -dn^2 - 1
b_1 = 2dn
b_2 = 0
b_3 = 0
b_4 = -d
b_5 = c
there are a couple ways to tackle this general form in algebraic geometry, finding the roots of the two equations. again some of the roots will be imaginary. the constraints should make things easier. i’m not sure if my 2-year-old will afford me the time to continue to work it out today.
BTW, i’m not seeing right how it can have 6 solutions. i think according to Bézout’s Theorem it should always have 4 including those in the complex plane.
[ed] oh. i missed your update where it is clear that you are actually solving for four hyperbolas under some kind of rotation. i don’t quite follow it but my intuition is that this is maybe not the most tractable way to formulate the problem. perhaps e.g. spherical coordinates would be simpler.
i definitely don’t feel like formulating it this way and plugging it into an iterative nonlinear solver in realtime is the way to go. just an intuition though.
[ed] oh, here are two math stackexchange questions treating the same problem. i haven’t quite grokked the answers in my few seconds of reading, and still need to work it out for myself.
looks like the one answer just plugs the equation into wolfram alpha and “magically” reduces to one variable… again i don’t quite see how without doing a lot of steps on paper.
[ed] ok, “magically” means computing the resultant, presumably with a Sylvester matrix.
here’s what that same technique yields given your specific constraints. it is a single (nasty) quartic equation for x, which has again a closed formula yielding four solutions.
a_5^2 b_4^2 x^4 + a_4^2 b_5^2 x^4
- 2 a_4 a_5 b_4 b_5 x^4 + 2 a_5^2 b_1 b_4 x^3 - 2 a_4 a_5 b_1 b_5 x^3
+ a_5^2 b_1^2 x^2 + 2 a_0 a_4 b_5^2 x^2 + 2 a_5^2 b_0 b_4 x^2
- 2 a_4 a_5 b_0 b_5 x^2 + a_2^2 b_4 b_5 x^2 - 2 a_0 a_5 b_4 b_5 x^2
+ 2 a_5^2 b_0 b_1 x + a_2^2 b_1 b_5 x - 2 a_0 a_5 b_1 b_5 x
+ a_5^2 b_0^2 + a_0^2 b_5^2 + a_2^2 b_0 b_5
- 2 a_0 a_5 b_0 b_5
= 0
a bit nasty, but computers don’t care. you can take those solutions for x and plug them back into the corresponding quartic for y (i guess.) omitting that here but here is the walpha command to produce it:
the difficulty now is finding the roots of the quartic, which is extremely complicated even for a computer…
honestly i’m surprised that it is this tricky, i guess i’ve only encountered degenerate cases before. in this case i still suspect there is some clever way to simplify the direct substitution and avoid this quartic insanity.
[ed] oh, i suppose i should also chase down the procedure outlined in this SE answer:
that seems more practical
[ed] i put this to my more properly mathy friend (well, my mom) and got:
… but. i just gave her the uncontextualized problem with a pair of quadratics, which is also where i started. @Rodrigo , i think this warrants a step back: i have not really grasped your initial formulation of the geometry problem. and:
i think we are not capturing all the constraints in the actual geometry, which might make things way simpler
it still doesn’t feel to me like this is the most elegant or tractable approach to estimating position from sensor DOA/TOA with a known geometry in 3D. (that is, trying to solve it with conic sections in a projected cartesian plane, but having to brute-force the equations with general nonlinear solvers.) maybe we could go back to a basic statement of the geometry of the problem, and try again from there. (and maybe without totally glorming up the thread.)
Oh man, what an epic response!
A lot of the maths goes comfortably over my head, but I’ll try to more fully unpack the process and context as that may be useful here.
///////////////////////////////////////////////////////////////////////////////////////////////
So as mentioned above, the idea is to triangulate the position of a strike on a drum using TDOA from multiple sensors.
The jumping off point is this paper from NIME proceedings back in 2014. In their approach they use 6 optical sensors on a drum and then use TDOA to find the position of a strike.
The main/relevant bits are:
and this:
Where S1/S2 are the physical sensors, and the hyperbola P represents all the possible places where the strike could have happened given the delta values of time of arrival.
They then follow that through to solve for each pair of sensors (S1/S2, S2/S3, S3/S4, etc…) to arrive at:
This is the jumping off point.
///////////////////////////////////////////////////////////////////////////////////////////////
So in my case I’m using completely different input sensors (a DIY/3d-printed thing I’m developing for SP-Tools, and only four of them.
That gives me a physical setup that looked like this at the time of capturing all my initial test data:
I did a load of testing moving in cardinal directions by fixed amounts etc… and that gives me audio that looks like this:
For the sake of simplicity (and for solving the harder parts first) I did manual/visual peak picking from this to arrive at delta offsets. In the final version this will either by cross-correlation or amplitude-based peak picking, but for now I’ve manually tagged deltas.
This particular set of offsets gives me: 0 80 125 82 (amount of samples, recorded at 44.1k).
At this point the math slightly departs from their paper as in testing, and reading through some of the sources with my friend (this hefty book) the way they are computing the speed of sound through the medium (s) produces some wildly different results. From the paper they are computing s as:
Which in our case would yield a speed of sound through the drum head as being 135.4m/s (r = 0.1778m (i.e. 7") and f = 248Hz. In our case, because I was thorough/detailed on how I was hitting the drum, we can actually compute s based on real-world measurements:
So with this, we get an s of 84.32m/s, which is much closer to simply computing 2rf, which produces 88.19m/s. This becomes important later on when it comes to the intersections.
After computing s and having delta offsets (which are converted to seconds, to have consistent units), that produces a quadratic equation per pair of sensor points (S1/S2, as per their paper).
This takes us kind of to the orignal point in that when solving for two positions (S1 and S2 in my picture above, which I have labelled cardinally as (N)orth and (E)east) I end up with the equations:
0.000171(x)^2 - 0.00023(y - 101)^2 = 1
0.00054(y)^2 - 0.00012(x - 101)^2 = 1
(these are in mm/s rather than m/s as I think that will be a more useful final unit)
When these two are plotted that yields this:
This would produce (up to) four intersections in theory:
(-76.8718, 94.2476)
(101.819, 43.0349)
(208.908, -66.6292)
(-406.089, -242.887)
But given the initial delta values of 0 80 125 82, that means that the (valid) intersection of these two particular hyperbolas would be (-76.8718, 94.2476) (top left quadrant intersection) as the sound arrived there first (the 0 in the list of deltas).
So the math/logic of this is that given the TOA of the wavefront at those two sensors (N and E) that (-76.8718, 94.2476) is the point that the strike happened.
This is then repeated for each of the remaining pairs (ES, SW, WN):
Then the mean of each of those intersections produces (-86.5685, 90.12725):
Which is really accurate as this was a strike 21mm diagonally in from point N on this plot.
I did have a mistake when I posted about 6 possible intersections as even though there are two more geometric intersections in this projection (red+green and blue+purple) those are incidental to the projection and those hyperbolas don’t actually intersect in “reality”. As in, I can compute North/South and West/East as additional hyperbolas, which would produce two more quadratic equations, but I have not (yet). Mainly because it’s not as straightforward to rotate/translate the hyperbolas onto this projection.
Now the reason I went on and on about the constant s is if I follow the paper to the T and use 2PIrf/2.045 I instead get this projection:
In this case two of the hyperbolas don’t intersect at all (red+blue, green+purple) and the ones that do intersect would produce a strike that would be physically off the drum.
///////////////////////////////////////////////////////////////////////////////////////////////
Phew, that was more than I was planning on typing out!
Was useful for me to write the whole process/thinking out though.
So yeah, that’s what’s going on, and why/how these equations are being arrived at. And for the moment, each set of points I’ve put into it returns a good/plausible result. The problem is trying to do this in “realtime” where a strike on the drum triggers an onset detection algorithm which then does a windowed cross-correlation or waits for 3 more peaks to happen in amplitude-land to produce some time deltas which are then shoved into all the maths after that, hopefully to return (“instantly”) an XY coordinate of where the strike happend on the drum.
yes, i’ve been properly nerdsniped. it is a fun little problem.
firstly, maybe moderators would like to see this whole subthread split off as it has little to do with Max per se.
so…
thank you for the explanation. i’m sure the equations are correct, not questioning that. but i still have a hunch that they could be reformulated to make a more tractable set of quadratics, possibly by some coordinate transformation at the outset. (e.g. it’s enticing that the sensors all lie on a circle…)
i’ll read the paper just to be sure i’m clear.
wearing an engineering hat, i’m always inclined to explore closed-form solutions before resorting to an iterative solver when the solution needs to be produced in bounded real-time. (in other words, there are strong benefits to having known numerical accuracy and known compute time for the solution.)
but here, that made more obvious sense when i had not spotted that the coordinate shift here leads to fairly disgusting quartics that can’t be easily solved using radicals. and the same engineering hat knows that there comes a point where the right thing to do is give up and fall back on less-ideal solutions.
so maybe indeed the best answer for the problem is to embed a library in max that can perform numerical root-finding (like BLAS, LAPACK and its descendants.) that seems like a generally useful thing to have and i’m a little surprised if no-one has done it before actually. (btw, it still might make most sense to have it crunch the two quartics - i more than half suspect the solver is already converting to that form.)
but… that is actually not an interesting problem for me, it is plumbing, the kind of task i’d happily perform for a reasonable consulting fee, but not something i’m gonna do for kicks in my free time. (sorry!)
no promises, but maybe i can go so far as to put these things (the quartics or whatever) in the form of a least-squares problem with a jacobian, so that they can be solved iteratively (what fsolve is doing under the hood before calling minpack routines, i assume.) the result would be C/C++ code calling the low level libs, and you could take over embedding that in a max external. (actually… nah, i just don’t feel like doing this for free, i’m sorry, it is too much work for a known result when other tasks are waiting.)
additionally (for fun), i think i will continue to pick at both trying to slightly reformulate the initial geometry problem, and perhaps pick at some of these algebraic methods of tackling quartics efficiently (because seems like a good thing to know about.) but no promises on timeline, it’s probably a weekend activity.
by the way, and not very directly relevant: i’ve noticed some historical errors in the paper related to prior work by Buchla: the Buchla “Thunder” did not use IR sensing, but a clever arrangement of capacitve sensors. they might be (somehow) thinking of the unreleased “Thunder II” which did use optical sensors to detect strike location on a drumhead. (a product was designed and prototyped to completion for release by E-mu, but i suppose the deal was never finalized.) like “Lightning,” “Thunder II” lacked 21st-century DSP resources and so employed more tricks in the analog domain: in “Lightning” the X/Y is captured quite directly by projecting through an aperture onto a photocell grid; in “Thunder II” a ring of alternating sensors and emittors are mounted below a reflective drumhead, thus capturing deflections in the surface with high accuracy and speed. (it’s almost the same idea as here, but using bespoke sensor/emitter arrangement i think makes the geometry nicer - i’m not sure it even had to use TDOA. i should take another look at it.)
the “Lightning” sensor geometry might interest you if you are interested in tracking LED emitters - it is a pleasingly low-tech solution that works quite well if (say) you can use a transparent drumhead and mount the sensor grid some distance beneath it. (doesn’t even need a computer - the most complicated part of the circuit is if you want to do what Don did and use specific PWM frequencies to distinguish emitters and filter out interference - then they must be demodulated but this is also accomplished by standard analog methods.)
ok one more update. i’ve found this paper that is a very specific treatement of hyperbolic position location using TDOA. section 3.1 gets into specific solver algorithms, noting that the nonlinear equations are hard to solve and outlining 3 methods for linearizing them:
[ http://archives.njit.edu/vol01/etd/2000s/2002/njit-etd2002-037/njit-etd2002-037.pdf ]
these are in 3-space BTW so a bit more complicated-looking. but i think same concepts apply. also worth noting that the final design involves a VHDL circuit for solving the final plane equations.
working through this section a bit more i think it is promising for your application: they are able to simplify by means of adding a 4th (technically redunant) sensor, and working all 4 equations at once. fortuitously, your application also has four sensors. there are 6 pages of algebra to work through, and it would be interesting to start it in two dimensions and try to acheive something similar.
i’m not 100% sure, but i guess it would also work to just use their results. make one set of coordinates zero throughout (let’s arbitrarily make it x given that they give everything in terms of z.):
and you have exact solutions in (z, y).
though, i don’t yet understand how this yields four solutions - i only see one radical. maybe over-determining it like this eliminates some of them already. extra nice if that’s the case.
