Input shaping

[tombrazier]

Description

Use record and replay of Bresenham algorithm calculations to effect ZV input shaping for X and Y axes.

It's quite rough and ready at present. At least the following are on my mind:

• It requires record buffers which consume more SRAM than I like. Longer term there are ways of doing this with no buffers.
• The "echo" comes after the main signal. I think there should really be two "pulses" one before and one after the logic for non-shaped axes for better alignment between the shaped and non-shaped axes.

How it works

Input shaping involves scaling a input control signal and overlaying it with one or more scaled and time-shifted "echoes" of itself to produce an output signal that damps oscillations. For example if an acceleration of 1000mm/s^2 started at time t, input shaping might change that to 500mm/s^2 at time t and then a boost of another 500mm/s^2 at time t + x. The value of x would be related to a resonant frequency of the physical printer.

The input signal could be anything that scales linearly. So it could be a time sequence of accelerations, as in the example above, or of velocities or of axis positions or even of stepper pulses. Input shaping works with all of these.

This solution treats calls to Stepper::pulse_phase_isr() as the input signal. For the shaped axes, it scales the effect of calling Stepper::pulse_phase_isr() and records and replays these calls. This is a simple place to insert the input shaping logic in Marlin and none of the rest of Marlin needs to know about it at all.

Stepper::pulse_phase_isr() is where the Bresenham logic happens. Each time the function is called it increments position counters for all the axes. Whenever a counter for any particular axis passes a specific value, a stepper pulse is generated for that axis. By overlaying two or more sequences of calls to the Stepper::pulse_phase_isr() logic, the rate of pulse generation is dependent on both the original input signal and the echo.

One concern that might spring to mind is whether overlaying two sets of logic that generate pulses might cause irregular timings of stepper pulses. Such irregular pulse trains cause uneven movement and are also known to cause some TMC stepper drivers to go to sleep. However this approach does not directly overlay pulse trains. Instead it overlays Bresenham logic and this results in clean pulse trains.

Here's a picture of how stepper pulses are generated for a move only in one axis, say X. The green line represents calls to Stepper::pulse_phase_isr(), the yellow line is the delayed echo which results in calls to Stepper::inputshaping_isr_x() (which duplicates the Bresenham logic), blue is the Bresenham counter which resets every second time it is incremented and red is the output stepper pulses which are generated whenever the Bresenham counter resets.

In that case the delay works out such that calls to Stepper::inputshaping_isr_x() do not causes stepper pulses but just make them happen on the next call to Stepper::pulse_phase_isr(). Here is an example where it's the other way around, i.e. Stepper::inputshaping_isr_x() generates pulses and Stepper::pulse_phase_isr() just increments the counter.

Either way we do not get double pulses or irregular pulse distances.

Changing direction

When there is a direction change for an axis between two segments, this complicates matters. With the echo, this direction change could occur anywhere in the time overlap between the segment ending and its echoes ending. This has two consequences:

• The Bresenham counter now needs to be signed and able to count in either direction.
• The Bresenham logic needs to recognise when the net effects of the original signal and the echo have changed from increasing to decreasing or vice versa.

As a first step, the logic in Stepper::block_phase_isr() which normally handles direction changes still does so if the echo buffers are empty when a segment is started. This will happen for the first segment printed in a run of segments. It will also happen for the first segment after a call to Planner::synchronize(). The latter is critical for homing where the direction change needs to be noted by the stepper module as soon as a bump segment starts printing. Otherwise the endstop logic will think he stepper module is still moving the print head towards the endstop.

More generally, the Bresenham logic can know for sure that a direction change has occurred if the Bresenham counter becomes sufficiently positive when the direction is backwards or sufficiently negative when the direction is forwards. This can introduce a minor delay in issuing the first step after a direction change but it does not appear to be problematic.

Compatibility

Input shaping introduces sections at the start and end of segments where the speed is combined with the speed of the previous/next segment or, in the case of a first or last segment, is halved. Whether laser power needs adjusting in these segments and how do so needs attention.

Input shaping may or may not work with direct stepping. This needs looking at.

For now, this PR introduces a sanity check preventing input shaping from being used with either of these two features.

The future

A better implementation would eliminate the buffer(s) that hold the echo data.

I think we could insert a step between the planner and the stepper modules which plans individual motions. One segment would produce up to three of these, for the accel, cruise and decel. So the stepper module would stop knowing about the three phases and only concern itself with single motions each of which is governed by one equation of motion. For an accel/decel phase the equation would be d = u.t + a/2 . t^2. (d is distance, u is initial speed and a is acceleration). For cruise phase it would be the same with a = 0. For SCURVE it would be the same but a higher order polynomial. The beauty of this is the IS can then be implemented simply by modifying this stream of equations of motion because overlapping signals can be summed just by summing the coefficients of the overlapping polynomials.

Testing / tuning

Tuning follows a similar process to that for Klipper.. We still need to document this for Marlin. Klipper's tuning tower can be used for tuning with M593. Print as a single or double perimeter with little or no infill. You may need to set the fan on full to reduce warping.

When slicing the tower, watch out that the slicer doesn't slow or smooth the corners. Things to watch out for:

• Set a high speed (80-100 mm/s).
• Set high accelerations (>= 1000mm/s^2) for X and Y. Watch out for the many slicer settings that can affect this. e.g. acceleration limits for print and travel moves and generate M204 codes as well as per axis limits that generate M201 codes.
• For both speed and acceleration, also look out for perimeter specific settings.
• Make sure the cooling settings do not specify slowing down to increase per layer printing time.
• Set the junction deviation as high as possible (with custom gcode M205 J0.3). Or, if using classic jerk, set the jerk settings as high as you can.
• Do not use the Arachne engine - it smooths out corners.

M593's parameters are

 *  D<factor>    Set the zeta/damping factor. If axes (X, Y, etc.) are not specified, set for all axes.
 *  F<frequency> Set the frequency. If axes (X, Y, etc.) are not specified, set for all axes.
 *  T[map]       Input Shaping type, 0:ZV, 1:EI, 2:2H EI (not implemented yet)
 *  X<1>         Set the given parameters only for the X axis.
 *  Y<1>         Set the given parameters only for the Y axis.

If using PrusaSlicer the following will print a frequency tuning tower, ranging acceleration from 1500 to 7000:

; Input Shaping Frequency Testing
{if layer_z==layer_height}M204 P1500; Set Acceleration at Layer 1{endif}
{if layer_z==5}M204 P2000; Set Acceleration at Layer 5{endif}
{if layer_z==10}M204 P2500; Set Acceleration at Layer 10{endif}
{if layer_z==15}M204 P3000; Set Acceleration at Layer 15{endif}
{if layer_z==20}M204 P3500; Set Acceleration at Layer 20{endif}
{if layer_z==25}M204 P4000; Set Acceleration at Layer 25{endif}
{if layer_z==30}M204 P4500; Set Acceleration at Layer 30{endif}
{if layer_z==35}M204 P5000; Set Acceleration at Layer 35{endif}
{if layer_z==40}M204 P5500; Set Acceleration at Layer 40{endif}
{if layer_z==45}M204 P6000; Set Acceleration at Layer 45{endif}
{if layer_z==50}M204 P6500; Set Acceleration at Layer 50{endif}
{if layer_z==55}M204 P7000; Set Acceleration at Layer 55{endif}

Paste the code block above into "Printer Settings" -> "Custom G-code" -> "Before layer change G-code" under any existing code you might have there.

Also set your X & Y max acceleration in your start g-code with M201 X7000 Y7000 and disable Linear Advance (if enabled) with M900 K0.

And the following will print a tuning tower for damping ratio (or zeta) with values 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9:

; Input Shaping Zeta Testing
{if layer_z==layer_height}M593 D0.1; Set Input Shaping Zeta at Layer 1{endif}
{if layer_z==5}M593 D0.15; Set Input Shaping Zeta at Layer 5{endif}
{if layer_z==10}M593 D0.2; Set Input Shaping Zeta at Layer 10{endif}
{if layer_z==15}M593 D0.25; Set Input Shaping Zeta at Layer 15{endif}
{if layer_z==20}M593 D0.3; Set Input Shaping Zeta at Layer 20{endif}
{if layer_z==25}M593 D0.35; Set Input Shaping Zeta at Layer 25{endif}
{if layer_z==30}M593 D0.4; Set Input Shaping Zeta at Layer 30{endif}
{if layer_z==35}M593 D0.5; Set Input Shaping Zeta at Layer 35{endif}
{if layer_z==40}M593 D0.6; Set Input Shaping Zeta at Layer 40{endif}
{if layer_z==45}M593 D0.7; Set Input Shaping Zeta at Layer 45{endif}
{if layer_z==50}M593 D0.8; Set Input Shaping Zeta at Layer 50{endif}
{if layer_z==55}M593 D0.9; Set Input Shaping Zeta at Layer 55{endif}

Tune zeta after tuning the frequencies and enabling IS. Use a fairly high constant value for acceleration so that you still see ringing even with IS.

Requirements

Any printer.

Benefits

Cleaner prints with less ringing.

Configurations

Related Issues

#16531

[tombrazier]

Oh I see we are both working on the tree Thinky. I have not included your changes and have done a forced push. Off to bed now, though.

[Riconec]

sorry for noob questions, is it something like klipper input shaping? How resonant frequencies will be measured?

[tombrazier]

sorry for noob questions, is it something like klipper input shaping? How resonant frequencies will be measured?

Yes, it's the same idea. This is a first implementation and I expect it will get more comprehensive as time goes on.

Resonant frequencies can, in principle, be measured with an accelerometer. But the easiest is to do what the Klipper documentation suggests and use a test print which demonstrates ringing.

[Riconec]

cool! I am up to testing when you have something working :D

[tombrazier]

Thanks @Riconec. I think it's close now. I'm getting a bit of drift that I need to track down.

[tombrazier]

Guess which one was printed with the code in 88223ba.

[tombrazier]

@Riconec It's ready for you to test. A bit rough and ready at present and you'll have to specify buffer sizes manually in Configuration_adv.h. But it seems to work.

[tombrazier]

@descipher Do you fancy having a look at how input shaping (this PR) affects lasers?

[Riconec]

damn, those are already improvements. I need to figure out how to connect adxl to measure it. and best would be to wrap my head around this algorithm because I still haven't figured out how it works. have you implemented it yourself or it is a port of some sort?

[descipher]

@descipher Do you fancy having a look at how input shaping (this PR) affects lasers?

I can do, I can tell you that AVR's will suck .. lol, its because AVR's already have timing issues even now since the initial introduction of laser capability. AVR's can do laser operations with only the basic on/off control and apps that generates gcode that compensates for accel/decel power linearity.

[tombrazier]

damn, those are already improvements. I need to figure out how to connect adxl to measure it. and best would be to wrap my head around this algorithm because I still haven't figured out how it works. have you implemented it yourself or it is a port of some sort?

Input shaping as an idea has been around for a long time. There are a bunch of papers out there which are all pretty mathematical. But it boils down to whenever anything changes, don't just apply the change all in one go. Apply part of the change and then a little later apply more, etc. The simplest of these is the "zero vibration" method which applies half the change now and half when the oscillation has reached its first peak. If you take damping into account, which this PR does not at present, then it's not half/half but some other ratio.

So, for example, if the print head is supposed to start accelerating at 1000mm/s^2, start accelerating at 500mm/s^2 and then a little later add another 500mm/s^2 to get to the full acceleration.

Anyway this is not a port in the sense of copying code from somewhere else. But it is not original either. If you want to tune it, see how it is done for Klipper.

[thinkyhead]

Theoretically these should be able to coexist, with either the direct stepping blocks being shaped or by having the direct stepping blocks skip the shaper altogether.

[tombrazier]

Agreed. They may even work okay as is. I think direct stepping would probably just not be affected by input shaping. But I don't have the capacity to check that out right now.

[thinkyhead]

I assume that Klipper's direct stepping protocol sends basic segments but the shaping happens on the printer's mainboard, so it makes sense to also do that for Direct Stepping. But maybe Klipper does some pre-calculation for the benefit of shaping before it sends those segments over.

[tombrazier]

Do we use the same block protocol as Klipper? Or is there a front end out there somewhere for feeding Marlin specifically?

[thinkyhead]

Direct Stepping G6 specifically requires https://github.com/colinrgodsey/step-daemon

[thinkyhead]

If shaping is optimized at all by using constants instead of configurable FREQ values, then it makes sense to make the editing of the FREQ values optional. If there are values we can pre-calculate from those edited frequencies, that will also help.

[tombrazier]

It's just shaping_delay_x and shaping_delay_y which are presently* currently constexpr. But I doubt there would be much or any performance difference if they were variables.

*"Presently" can mean both "soon" and "at present" in the UK - is it different in the USA?

[tombrazier]

I am also pondering adding an option to specify damping. The papers on input shaping take that into account by changing the ratio of the pulse trains. e.g. they could be something like 60/40 rather than 50/50. There's a formula related to damping ratio ζ but since I don't think users are going to know how to measure ζ I am thinking we just print a tuning tower with different ratios.

[tombrazier]

I have kind of run out of energy to write the runtime editing code. Anyone else fancy the job? Also cribbing Klipper's tuning guide.

[descipher]

Initial observations indicate it works with basic laser operations. I will look further into it and update you when I can. Not likely 2b this week.

[truglodite]

The above mentioned suggestion of reading the Klipper tuning guide and applying it to marlin IS I think leaves a lot of room for differing interpretations. Overall for Klipper the procedure goes 1) find resonant frequencies, 2) pick an IS model, then 3) find the max acceleration for that model based on ringing and smoothing. Currently Marlin does not allow selecting an IS model, and it uses zeta as a tuning parameter.

Using that mix of info, I tried something like this: 1) measure x and y resonance, 2) using those frequencies, print a variable zeta tower ~1000mm/s^2 below where skipped steps occur in step 1, 3) print another variable accel tower using the best zeta from step 2 to find max accel, 4) repeat 2 using different accel and 3 using different zeta to fill out the desired scope of testing.

The problem is the variable zeta tower didn't seem to have any sweet spot. A large part of that is probably the increasing height of different zeta values, and the resulting interaction with my printer frame (more ringing amplitude as you raise higher above bed). I recall reading something in the Klipper docs that they keep zeta fixed at 0.15. As is, it seems that the most productive method will involve tons of iterations using different acc/zeta to fill out a test matrix. I was hoping to learn a more efficient method, if it exists.

Should we be testing zeta 0-1, or maybe closer to 0.15 like 0.1-0.2?

Also, is my above iterative flow for tuning marlin OK? ...or should I add/modify any steps?

[tombrazier]

@truglodite Your flow is what I would recommend. I found that varying zeta between 0 and about 0.3 there is little difference but above that print quality degrades. Zeta of 0.15 does seem pretty common.

[tombrazier]

@VanessaE Interesting. Your reasoning makes sense to me. If the difference between X and Y is not too great you should be able to use the mean value for both motors and still get significant benefit. A more aggressive shaper, such as EI, would help even more.

Better, however, would be what you suggest, i.e. shaping per axis rather than per motor. This is not possible with the present implementation. I do keep thinking about an implementation that does the shaping in the planner rather than the stepper module (which would also allow us to ditch the buffers!) but it would mean we could no longer use bresenham because on corners straight lines turn into curves with IS. Anyway, shaping in the planner could work per axis or per stepper.

It might make more sense just to wait and see what FBS brings.

Delta also has symmetry across its axes (all three - for which IS is not currently enabled) and any variation in frequency will be to do with something else. So there is likely to be parallels when we do implement IS for three axes on delta.

[VanessaE]

@tombrazier

If the difference between X and Y is not too great you should be able to use the mean value for both motors and still get significant benefit

Unfortunately, the difference is pretty significant on my smaller machine, but if I'm being honest, that printer can do without IS anyway (it's just more convenient to test stuff on than my big machine).

My big machine is the one that really could use IS, but there's no chance of using your idea, because while its heavy toolhead makes the machine's X and Y moving masses more similar than on the small printer, that toolhead is top-heavy, which causes it to ring only in the Y direction, due to the X gantry having some flexibility, which allows the toolhead to twist, if it gets enough of a Y impulse (and that's after having done everything I can to stiffen the X gantry and improve the toolhead's center of mass).

Guess I'll have to wait.

[j0hn1k]

@VanessaE Interesting. Your reasoning makes sense to me. If the difference between X and Y is not too great you should be able to use the mean value for both motors and still get significant benefit. A more aggressive shaper, such as EI, would help even more.

Better, however, would be what you suggest, i.e. shaping per axis rather than per motor. This is not possible with the present implementation. I do keep thinking about an implementation that does the shaping in the planner rather than the stepper module (which would also allow us to ditch the buffers!) but it would mean we could no longer use bresenham because on corners straight lines turn into curves with IS. Anyway, shaping in the planner could work per axis or per stepper.

It might make more sense just to wait and see what FBS brings.

Delta also has symmetry across its axes (all three - for which IS is not currently enabled) and any variation in frequency will be to do with something else. So there is likely to be parallels when we do implement IS for three axes on delta.

problem to make IS on the axes related to the performance of the controller?

[BEprints]

A few noob questions, Do I need to enable acceleration and jerk in Cura when slicing the test tower I have the post process set and what does it mean by test Zeta

[tombrazier]

@BEprints I can't answer the Cura question but zeta is the "damping ratio" which specifies how quickly vibrations die down on your printer. It is a value between 0.0 and 1.0. On most printers it seems to be around 0.15. A zeta test tower will vary zeta with height and you pick the zeta value that gives the least ringing.

[BEprints]

Thanks, I have figured out the Cura side for the acceleration I used but not sure on damping side I seen there was a screenshot a bit further up with a post processing script in cura for input shaping so I will look into that

[BEprints]

I am the same nothing will show I have the ringing tower stl from klipper I just need the post processing script

[TheSunshined]

Man, first of all.
Thank you for making marlin better for all of the users, saying as ender 3 user with Profesional Firmware fork.

[ahmetcemturan]

Thank You so much for this. I was seriously considering moving to klipper ( Especially with the current Raspberry Prices ) but that idea has been cancelled now. Just surprised that there is hardly any wuzz about it. I stumbled on the update by chance and there is hardly any videos or tutorials about it.

[tombrazier]

@ahmetcemturan You're welcome. Just implementing it was a bunch of work. Effort now needed to popularise. But even before that we need to get the docs in order. The present recommended tuning guide doesn't work as well as a modified klipper tuning procedure.

[ahmetcemturan]

Yeah documentation is very important.   I consider myself experienced but did not notice that input shaping needs to be enabled in the adv.config...  I printed the test and wondered why there was no difference.  :))-- Ahmet Cem Turan+90 555 455 544417:58, 23 January 2023, tombrazier ***@***.***>: @ahmetcemturan You're welcome. Just implementing it was a bunch of work. Effort now needed to popularise. But even before that we need to get the docs in order. The present recommended tuning guide doesn't work as well as a modified klipper tuning procedure. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[ahmetcemturan]

One thing though: The Dampening factor, is it not implemented yet? I tried to run this "M593 D{(layer_num < 2 ? 0 : 0 + 1.00 * (layer_num - 2) / 297)} ; Hz Input Shaping Damping Test 0-100" with the ringing test file after I had found my frequencies.
Unfortunately I could not notice any difference throughout the print.
Atatced a photo of the print:

Also: In the Config file there is an f after the damping value. What is that for? Do I need to enter it if I want to change D via Gcode?

[tombrazier]

It is implemented. Its effect is fairly subtle. The ringing in that test print is definitely more pronounced at the top. I think it is lowest at about 10mm.

[ahmetcemturan]

That's what I thought too but that is D 0.01
This is a CoreXY machine btw. I am now trying to adjust without turning the test 45degrees and using same F value for x and y.

[Mjankor]

I noticed a difference on the dampening factor. IIRC it had larger gaps in the test print and smoothed corners at higher damping rates.
Here's a couple of test prints on a CoreXY. Mine has a fairly light gantry and carriage, so I had pretty similar frequencies on both X and Y axes, and on the per-motor "axes".
Speeds are (I think) 150mm/s and 1500 to 7000mm/s accel.

[tombrazier]

I think we are coming to the conclusion that turning the model 45° for CoreXY is a red herring. Resonant oscillations happen characteristically on X and Y cartesian axes on CoreXY, as they do on bed slingers.

In the present implementation of IS, the axes are referred to as X and Y but in reality they are actually for stepper motors A and B which, of course, are directly correlated with X and Y on a bed slinger but not on CoreXY (or delta). We should really have referred to A and B (mia culpa). (And for deltas we need to implement IS for C.)

I think @ahmetcemturan you are correct in using the same F value for X and Y (really A and B) as CoreXY systems are symmetrical for A and B steppers. You will need to pick a frequency around half way between and X and Y resonant frequencies. This should still give some benefit but not as much as an IS implementation that specifically targets X and Y in cartesian space.

[ahmetcemturan]

I found that the best way to find the resonant frequency is to measure the distance of the "waves" then calculate the frequency. Then I can run the ringing test again with a better resolution.. E.g measured/calculated 32Hz. 28-34Hz ringing tower.

[ahmetcemturan]

I use this Gcode: M593 F{ ( 25+(30 * int(layer_z)/60))} it is better defined as it switches every mm and is easier to calculate afterwards. Just measure and put in the height.

[Chris-Jayden]

Would the implementation be any different for systems with Core-XZ kinematics?
We have X-Axis/Motor Input Shaping.
We have Y-Axis/Motor Input Shaping.
Would we then need Z-Axis/Motor Input Shaping?

[tombrazier]

We do need to replicate the X and Y axis code for Z. Delta machines need this too. It's not hard to do but really needs a bit of a code refactor using templates to remove all the existing code duplication for X and Y.