Scanning speed and velocity mapping of Cybrid

I'm trying to sort out how all the various filters described here compare with a simpler approach. However I see a number of options, have we settled on a "best one"? Perhaps the one @JayKominek is utilizing which is

https://github.com/jkominek/piano-conversion/blob/a8b33c9164264ce179180f5ef9f337700e9bbd1d/firmware/adc-board/Core/Src/main.c#L134-L142

#define FIR_LEN 13
// TODO want this in DTCM also
// scipy.signal.savgol_coeffs(13, 2, deriv=1, use='dot', delta=0.01)
// the delta doesn't matter, so it was chosen to bring the values to ~1.0
float filter[FIR_LEN] = { -3.2967032967033107, -2.7472527472527473,
		-2.197802197802188, -1.64835164835163, -1.0989010989010746,
		-0.5494505494505213, 2.9531932455029164e-14, 0.5494505494505779,
		1.098901098901124, 1.6483516483516676, 2.1978021978022086,
		2.7472527472527473, 3.296703296703284 };

and

https://github.com/jkominek/piano-conversion/blob/a8b33c9164264ce179180f5ef9f337700e9bbd1d/firmware/adc-board/Core/Src/main.c#L1035-L1044

static inline __attribute__((always_inline)) float computefilter(int chan, int idxs[])
{
	float v = 0.0;
	for(int i=0; i<FIR_LEN; i++) {
		// TODO check to see what the compiler-emitted assembly looks like.
		// if we're not using the fused multiply and add instruction, fix it.
		v += filter[i] * sensorbuffer[chan][idxs[i]];
	}
	return v;
}

Consider that:

• We need the velocity only at the time of strike (before it is slowed down by whatever the hammer strikes, obviously)
• The hammer is in free flight from let off to the time of strike
• The previous two points imply that the hammer velocity is approximately constant and more accurately slowed down by gravity

My question is: instead of doing this hyperlocal measurement at the very end (with all the noise implications of measuring only few points all nearby), why we could not measure a few points from let off to strike and make the measurement basically averaged along a much longer travel (possibly considering gravity slow down, for which we have analytical formula)?

And by "few points from let off to strike" I really mean a few points (again, to average noise) and not one point alone as I am currently making (if you looked at my code). My gut feeling tells me that this would be computationally simpler (not a big deal on the STM32, but somewhat on a pico, see below) and probably more accurate for the purpose we need.

Thoughts?

vagfilm No idea about implementing Kalman on pico (Python?)

Only now I noticed that I never responded to this. I guess you can use Python on the pico, but I am using C (like Jay is on the STM32). And, to reiterate, the pico is a testing platform which might (and hopefully will) turn out to be decent enough to be viable for those who don't want to order custom PCBs, but it is not our final goal.

RIP My gut feeling tells me that this would be computationally simpler (not a big deal on the STM32, but somewhat on a pico, see below) and probably more accurate for the purpose we need.

The "linear difference" mechanism you're suggest requires deciding beforehand(/algorithmically) where the start and stop points ought to be.

The "compute a velocity constantly when you're close enough to the peak to matter" trades away CPU cycles, which I have plenty of, for the very simple decision procedure. "Is this the peak velocity? Yes, generate note. No, keep waiting for peak velocity."

I've been nearly hopping with excitement ever since I realized how much simpler the logic is when I've got a constant stream of velocities to work with. I wouldn't return to what you're suggesting unless somebody else wrote the code for it, and then wrote automatic calibration procedures, and proved it all worked better.

RIP with all the noise implications of measuring only few points all nearby

That 13-point Savitzky-Golay looks smooth enough on my data that you could convince me it was generated by some purely analytical process, not subject to noise. I've got the processing power to extend that into 20-something points, easily.

RIP have we settled on a "best one"?

I don't see a need to decide on a best one. I (at least) just need to find a "good enough" one (for now). I'm doing engineering, not science. 😃

I consider the Kalman filter interesting as well. I looked at that sometime back in the fall, but my ADC readings were still too rough, and it couldn't overcome the noise. It might be feasible, now. Since it could also produce a running estimate of velocity, the same simple decision procedure would apply, and there'd be a minimum of calibration.

JayKominek the peak velocity? Yes, generate note. No, keep waiting for peak velocity."

You don't take into account any absolute position? How can you detect a silent strike without accounting for position?

vagfilm

So, I wrote a little quickly since I was addressing DelVento who'd been reading the relevant code. The slightly longer version is, I only run the velocity calculation when the position is "high enough". I consider that more of a performance optimization than something necessary for correctness.

(Hand waving, without reviewing my data: ) Quiet/partial/incomplete/whatever strikes produce wildly low note velocities which convert to MIDI notes with velocity <1. For debugging purposes I round those up to 1, but when the code is more finalized, they'd just be dropped, as even with high resolution MIDI, you can't go below 1.0.

But by only running the velocity code when the hammer is high enough, the code wastes less time, and processes fewer events which would just be dropped later.

JayKominek The slightly longer version is, I only run the velocity calculation when the position is "high enough"

Ok. Thanks for clarifying. I think that is a good approach...
If my logic is correct, when the hammer lefts the resting position and does not return to resting position but also does not reach the max position, then you can send a note-on == 1, and some VST's use that as silent strike with the damper lifted and let the string resonate. Position of the hammer (or position of the key) is important for some events.

On the subject of lightly played, near-silent and silent notes: The velocity peak is when the hammer escapes the key into free flight, but it's slowing down all the way from there to the string (if it even reaches the string). This may be a small difference for loud notes but the quieter a note gets the bigger the difference will be. Maybe better to take the velocity from a position calibrated to be just below the string?

One very basic newtonian question (told you before that I am no expert in physics, even basic things...).
Velocity of hammer at impact always equals sound, or momentum/acceleration plays a role? I mean, if the hammer hits the string at velocity N when played from the jack (ie, hammer accelerating) does it cause the same string displacement (loudness) as an hammer at velocity N that is in free flight from rest position because of a staccato hit (and is probably deaccelerating)? Maybe in practice this does not matter (off the jack may never reach the same impact speed as a strong free flight hit)...

vagfilm Velocity of hammer at impact always equals sound, or momentum/acceleration plays a role?

The difference between velocity and momentum is just the hammer mass, which all else being equal, would certainly matter. You can expect a 20g hammer striking at 6m/s to sound louder than a 10g hammer. However, this mass matters in the same way that string tension, string weight, felt damping, hardness and shape etc, all matter. These are all essentially the same from blow to blow on a given string, so all would get lumped into the string model.

vagfilm I mean, if the hammer hits the string at velocity N when played from the jack (ie, hammer accelerating) does it cause the same string displacement (loudness) as an hammer at velocity N that is in free flight from rest position because of a staccato hit (and is probably deaccelerating)?

This is a slightly different discussion, with perhaps some controversy. In the very simplest model, the hammer is in free flight when it hits the strings regardless of how you attack the key, so the only thing the string can see is the final velocity of the hammer. How you accelerate to get to that velocity, whether as an impulse, or an extended press is irrelevant. As long as you're still holding the key depressed after the strike, I think this is an adequate picture.
However, if you're considering different staccato hits, there can be differences in the timing of when the damper returns to the string relative to the hammer hitting the string. The quick strike can let the damper return to the string sooner and harder for a shorter note. I believe that for very short notes, the perception of loudness is also affected by the note length, so it's not quite apples to apples comparing identical strike velocities. However, if you want to compare notes of comparable loudness, you can have a short duration, high velocity key strike, or a longer duration, lower velocity strike, having the same 'loudness', but different tonality.

RIP And by "few points from let off to strike" I really mean a few points (again, to average noise) and not one point alone as I am currently making (if you looked at my code). My gut feeling tells me that this would be computationally simpler (not a big deal on the STM32, but somewhat on a pico, see below) and probably more accurate for the purpose we need.

The 'few points away from strike' boils down to a trade-off. Fewer points = more bandwidth = more noise, but also less latency.

My feeling is Savitzky-Galoy is overkill, but if the drawbacks aren't show-stoppers, it's perfectly reasonable to accept as being good enough and move on to other matters.

The drawbacks are the computational complexity and latency. In what I see of @JayKominek 's implementation, you have to keep 13 past samples around, and a floating point multiply and addition on each of them, for every cycle. The complexity is a moot point if the platform can provide it, but as you point out, can be an issue say, on the pico. The other matter is latency - the 13 point FIR, assuming it's linear phase, is a 6.5 sample group delay. I don't know what sampling rate this is at, and haven't played with his data, but this is probably neither a show-stopper nor immaterial.

With @gzpiano 's data, last I played with it, it looked like two single-pole FIR filters (each just one state, and a shift/add - this can all be fixed-point) is adequate for smoothing out any peakiness in the velocity for his original 10-bit data. The coefficients I last looked at, I think gave a group delay of 5 cycles, so not significantly better than Jay's, but the simple fixed-point calculations can be a win. OTOH, with the same data, a shorter Savitzky-Golay filter may be adequate too, bringing its delay down too, I haven't tried.

Thanks @luns for the reply. It was strictly a naive Newton' physics question, and I was maybe not clear, but I was referring to differences in free flight travel on the same hammer (so, no change in mass).

luns I believe that for very short notes, the perception of loudness is also affected by the note length

As I said above, my question did not relate to this, but you are correct. To an extent, perception of loudness depends on volume contrast, so very short notes are perceived as louder than slightly longer ones.

luns The other matter is latency - the 13 point FIR, assuming it's linear phase, is a 6.5 sample group delay.

By default scipy.signal.savgol_coeffs generates coefficients which evaluate the polynomial at the centre of the window, but you can get it to generate coefficients to evaluate any position you like including the endpoints:

scipy.signal.savgol_coeffs(13, 2, deriv=1, use='dot', delta=0.01)
generates @JayKominek 's original coefficients, FIR computes derivative at 6 sample delay.

scipy.signal.savgol_coeffs(13, 2, deriv=1, use='dot', delta=0.01, pos=12)
generates different coefficients, FIR computes derivative at "12th" position in window, i.e. current sample.

It's not a free lunch though, the filter becomes worse in other ways (noisier and can overshoot rapid changes) if you move right to the end of the window. But tuning it for things like having a slightly longer FIR yet with a slightly shorter delay on the filter is possible.

JayKominek I've been nearly hopping with excitement ever since I realized how much simpler the logic is when I've got a constant stream of velocities to work with. I wouldn't return to what you're suggesting unless somebody else wrote the code for it, and then wrote automatic calibration procedures, and proved it all worked better.

Well, writing the code is easy, but calibration is tedious and not fully automatic as you say. I have a basic idea but it requires to do a couple of hammer strikes on each key. Getting rid of that is an excellent point which I underappreciated until you stated it so clearly in this message! Thanks!

For the sake of discussion and understanding, I still would like to get an answer to my question, which was slightly different (and I'll get to it in a moment).

luns The 'few points away from strike' boils down to a trade-off. Fewer points = more bandwidth = more noise, but also less latency.

What do you mean less latency? (*) Assuming you have enough compute horsepower to do the math (which I think the STM32 does, perhaps not so on the RPiPico), you can do more points, say 50? 100? Again, these are just "a few".

I mean, we are measuring the hammer from let off (and even before) to strike, that is at least 1mm, perhaps 2mm or more. Ballpark that should be 1000 samples, perhaps more, and we are throwing 90% or more away and utilize only the last few dozen(s) to compute the final velocity. Well, the final velocity is the only thing we need because it is what it causes the string to vibrate in an acoustic, so I am fine with that. What I am asking (and this reconnects with what I was saying above to Jay and in my "stirring" post of this morning) is if this is the best way to proceed (considering also measurement noise) or if utilizing earlier samples could provide better estimates of this end velocity -- regardless (for now) of other complexities, e.g. calibration as Jay mentioned in his answer.

In other words: knowing what we know, i.e.:

• the hammer is in free flight, being slowed down only by gravity (and perhaps friction, which we can also model relatively easily)
• frequency, resolution and noise of the ADC (ok, we don't know much of noise characterization, but still)

Is the best estimate of the final velocity the one which throws away the vast majority of the samples? Perhaps the answer is yes, but I don't know and I don't see it being an obvious yes. Or perhaps the answer is "we don't know" or "it's a close call" or something like that. Or even "no, it's not the best, but it's good enough" and hence we still do it that way, because (as Jay said above) it is so much better from a continuously-self-calibration point of view and that trumps everything else.

Thanks guys (even ones who I am not responding directly, I read you all and appreciate the conversation). This is such a nice discussion both from a pure inquiry point of view, as well as the practical implications of building these instruments!

(*) I assumed that we were already doing what @xooorx just wrote as position in position 12th -- I'm gonna reply to that message separately

xooorx luns The other matter is latency - the 13 point FIR, assuming it's linear phase, is a 6.5 sample group delay.

By default scipy.signal.savgol_coeffs generates coefficients which evaluate the polynomial at the centre of the window, but you can get it to generate coefficients to evaluate any position you like including the endpoints:

scipy.signal.savgol_coeffs(13, 2, deriv=1, use='dot', delta=0.01)
generates @JayKominek 's original coefficients, FIR computes derivative at 6 sample delay.

scipy.signal.savgol_coeffs(13, 2, deriv=1, use='dot', delta=0.01, pos=12)
generates different coefficients, FIR computes derivative at "12th" position in window, i.e. current sample.

It's not a free lunch though, the filter becomes worse in other ways (noisier and can overshoot rapid changes) if you move right to the end of the window. But tuning it for things like having a slightly longer FIR yet with a slightly shorter delay on the filter is possible.

Depending on the sampling frequency, 6 sample delay might be irrelevant. For example at a relatively slow 10kSample/s 6 samples is 0.6ms (I believe @JayKominek is planning to do 20kSample/s per channel, i.e. 0.3ms).

Can you elaborate on the tradeoffs for using the last position (or one before last or the like if much better)? Not in the abstract, but considering the situation? In a dozen sample, I believe the hammer can approximated by being at constant velocity (or at most constant acceleration) so no rapid change expected, and the only noise is the one from the ADC itself (since the hammer will be close enough to the sensor that all the other sources will be screened away).

Thanks!

RIP luns The 'few points away from strike' boils down to a trade-off. Fewer points = more bandwidth = more noise, but also less latency.

What do you mean less latency? (*) Assuming you have enough compute horsepower to do the math (which I think the STM32 does, perhaps not so on the RPiPico), you can do more points, say 50? 100? Again, these are just "a few".

By latency, I mean the group delay of the filter.

When you're smoothing out a curve, to determine what a given point of the smoothed curve should be, you'd want to look at samples near that point. Not only do the samples before the point inform you of where the point ought to be, but samples afterwards have information too, so it's nice to be able to look into the future. However, if you're processing samples as they come in, you know everything up to and including the present time, but nothing of the future.

Let's say you want to estimate a smoothed value of sample n based on 6 previous samples, and 6 'future' samples, how do you do that? You wait until sample n+6 is available and retroactively figure out based on samples n-6 through n+6 what the smoothed value at sample n ought to be. However, it's now 6 cycles after when sample n was current. Or less hypothetically, you just finished estimating what the hammer position was 6 cycles ago.

This is what FIR filters (finite impulse response) typically do, the default Savitzky-Golay filter included - they look at a finite number of samples symmetrically before and after, together with a delay to make things causal, and this is normally done specifically to allow this symmetry. While it's possible to make asymmetric FIR filters, if you're going to take into account more past samples than future samples, there's no need to limit yourself to a finite window; you can look at samples going infinitely far back into the past. This is known as IIR.
It's not possible to re-evaluate all the (infinitely many) past samples one-by-one at each time step the way you do for FIR. Instead, you maintain a state that captures what you need of all that history, and let that state decay exponentially while you replenish it with information from new samples. The weight with which you consider the past samples falls off exponentially instead of abruptly going to zero at the edge of an FIR window. The time constant of the decay is in some ways comparable to the length of an FIR window, also contributing a group delay (latency) in a similar manner. Also, while the weighting of far past samples necessarily falls off exponentially, that exponent doesn't have to carry all the way up to the present: you can also have what's essentially an FIR filter to shape things closer in.

RIP Can you elaborate on the tradeoffs for using the last position (or one before last or the like if much better)? Not in the abstract, but considering the situation?

The tradeoffs will be against whatever the exact situation is:

How many bits are the ADCs? More bits = less quantisation noise = less smoothing required.

What's the sample rate? The higher it is the more aggressively you can smooth away the noise without damaging the signal.

What accuracy is required? Getting a number that's a good proxy for 127 velocity levels will be easier than, say, getting a 12-bit measurement of the actual objective velocity.

W.R.T. the position variable in the filter, here's the original centered filter (red) vs the pos=12 version (blue) showing the faster response, but also the increased noise and the overshoot:

luns Let's say you want to estimate a smoothed value of sample n based on 6 previous samples, and 6 'future' samples, how do you do that? You wait until sample n+6 is available and retroactively figure out based on samples n-6 through n+6 what the smoothed value at sample n ought to be

Does one only consider samples of a given sensor over time or also those of other keys? That is, taking into account the overall velocity trends of a musical piece.

n-player Does one only consider samples of a given sensor over time or also those of other keys? That is, taking into account the overall velocity trends of a musical piece.

It is up to us to decide. AFAIK (and certainly for what I am doing), currently I consider only the sensor itself and not the others. Using others would be difficult for me, and with my minimalist approach I would do it only if strong benefits are shown.

xooorx W.R.T. the position variable in the filter, here's the original centered filter (red) vs the pos=12 version (blue) showing the faster response, but also the increased noise and the overshoot:

If my understanding is correct, the black is the hammer position, right? If so, even the red curve may be more than adequate, because in an acoustic piano the sound is not heard immediately as the hammer starts to slow down upon touching the string, but only after it leaves it free to vibrate, so when it starts going down.

Regarding the overshoot, if that is an additive or multiplicative factor more or less constant for all the MIDI velocity levels, it's irrelevant, no? Of course in general the overshoot will depend on the characteristics of the thing being differentiated, but here we know well: a slowing down hammer, so we can use that knowledge to make assumptions?

xooorx How many bits are the ADCs?

The least accurate ones (mine) have 9 bits effective (12 nominal). The most accurate ones, 16 nominal (I dunno effective).

xooorx What's the sample rate?

The slowest I've seen discussed is 1kSample/s (per channel). The faster (mine) 167kSample/s (per channel).

xooorx What accuracy is required?

At the very least 0.5 MIDI velocities, i.e. 8 bits. Ideally some more, let's say 14 bit (but that is likely overkill).

With this info can you provide more specifics? 🤣

I think one thing that we have not said is what is the slowest and fastest velocity which we would like to detect. We have that paper quoted above in this conversation, but I also captured an ADC session with pp, mf, f and ff strikes. I have posted that data into here and it can be plotted or extracted to a text file with the ./plot_midi.py tool placed here -- I will add READMEs shortly after hitting "post reply" here.

The pro is that my data contains also timing information of when the sampling occurred (with us accuracy). The con is that the USB bus of the Pico (and perhaps also of the old receiving computer I am using) chocked if I tried to send it at 10kSample/s or faster. Here is the picture of one of the three data files

RIP I have posted that data into here and it can be plotted or extracted to a text file with the ./plot_midi.py tool placed here -- I will add READMEs shortly after hitting "post reply" here.

READMEs added. If anyone has troubles with the python and would rather have just the text files, let me know with a private message and I will post them somewhere too.

RIP With this info can you provide more specifics? 🤣

Probably, eventually, yes.

So far I've coded up a physical model of a key/hammer/string system (pretty basic but it's high resolution and high sample rate) which is then "sampled" by a model of a sensor & ADC e.g. 12-bit, 16KHz, 2% noise or whatever.

The next step is to try out the filtering we've been discussing on the modelled ADC data and see what can be achieved. The point of doing this with modelled sensor data instead of real sensor data is that we also get the ground truth from the underlying model of what the position and velocity actually were at any instant, so can objectively measure how well the filtering is recovering these from the sensor data.

I'll post code as soon as I've got some filtering going. Meanwhile here's a picture of some entirely artificial sensor data:

xooorx So far I've coded up a physical model of a key/hammer/string system

Fantastic! Now we are getting scientific 🙂

If you care, you can introduce in your model the sensor response too. In fact the ADC number is not linear with distance for three reasons:

• current from the sensor (ADC measures voltage but V=RI is quite linear)
• ADC conversion (both in value and in time)
• time response for the sensor

For the first reason: see Pag.4, top figure for the 'RCA6 sensor (not sure what "negative distance" means in this plot). With this sensor (the one I am utilizing) at the distance of most interest, i.e. the last few mm, it's linear enough you can probably ignore this issue (assuming the datasheet is correct, which might or might not be). For the CNY70, the difference is big, see fig. 9 (notice that unlike the other one, the vertical scale is logarithm!) and also notice how it flattens in the most important last half mm. The same problem is shown in fig.2 for the QRE1113.

The second problem is described in general at https://microchipdeveloper.com/adc:adc-inl and the specifics of the ADC I am used are on fig. 116 and fig. 117 (in pag 588-9) and in the Errata RP2040-E11 on pag 647 of the datasheet (you would not think a datasheet would have 600+ pages eh?) -- I said repeatedly, I am interested in characterizing this platform, but I expect that eventually it would show its limits and not be adequate. So you might want to consider the other MCU which @JayKominek is already utilizing and which I plan to switch to when he has settled on a hw design. Here is the datasheet and the relevant data is in fig. 39 at pag 171 and fig 91 at pag 284. @JayKominek may chime in with more details, since I did not dig too deeply into these 356 pages (less than the other datasheet, but denser…)

Then, there is the time factor. The sensors do not respond immediately to a changing light, but they have a bit of delay. The sensors' datasheets report information about this either in a figure or in a table. Links above, let me know if the way they report it is too cryptic for you

The data is posted in https://github.com/davidedelvento/Mybrid/tree/main/data/RPiPico%2BEAITRCA6 in binary format.

There was a file uploaded here as plain text file (UNIX end line) rather than MIDI, but it got inadvertently deleted as described at https://pianoclack.com/forum/d/1222-i-have-inadvertently-deleted-some-files-on-the-forum/17 and I am not reuploading unless asked