# Reducing harmonics in the resonance testing signal

Currently, the velocity of the resonance test signal is a triangle wave, with square wave acceleration.
By incorporating brief constant velocity sections into the test signal, it is possible to achieve a better sine wave. All we have to do is reduce the maximum velocity:

``````diff --git a/klippy/extras/resonance_tester.py b/klippy/extras/resonance_tester.py
index db03e3207..64b9afcc7 100644
--- a/klippy/extras/resonance_tester.py
+++ b/klippy/extras/resonance_tester.py
@@ -91,11 +91,11 @@ class VibrationPulseTest:
gcmd.respond_info("Testing frequency %.0f Hz" % (freq,))
while freq <= self.freq_end + 0.000001:
t_seg = .25 / freq
-            accel = self.accel_per_hz * freq
-            max_v = accel * t_seg
+            accel = self.accel_per_hz * freq * 1.08813
+            max_v = accel * t_seg * 0.74202
"M204", "M204", {"S": accel}))
-            L = .5 * accel * t_seg**2
+            L = 0.46672 * accel * t_seg**2
dX, dY = axis.get_point(L)
nX = X + sign * dX
nY = Y + sign * dY

``````

I obtained this factor by maximizing the SNR of the fundamental harmonic over higher harmonics. The optimal duration of the acceleration is 0.742 time shorter than in a pure square wave. There is no closed form solution, but you can ask for the derivation details. This should yield a 40% reduction in the ratio of higher harmonic power to fundamental harmonic power.

Hereâs the square-wave acceleration signal currently used by the resonance tester:

SNR: 2.07

• blue: square wave signal,
• green: fundamental harmonic,
• orange: sum of the 10 following harmonics, acting like the differences of the blue and green curve with a low-pass filter.

Introducing cruised segments between accelerations reduces the higher harmonics:

SNR: 3.45
(note the segments with zero acceleration around x=0.25)

Velocity signal:

• orange: current velocity signal,
• blue: velocity with cruised sections,
• green: true sine wave (fundamental harmonic)

The effect on the position signal more subtle.

The question arises as to whether we should reduce these higher harmonics. I donât believe that their presence is negatively impacting the accuracy of the test. Currently, the accelerometerâs spectra are averaged throughout the entire sweep, rendering it inconsequential which vibration mode is stimulated at any given moment. When we stimulate at 50 Hz, a certain amount of energy is inevitably distributed to 150 Hz, 250 Hz, and so on, with the resonance amplitudes contributing to the final results within these frequency bins. This leads to an increase in gain at higher frequencies, but this outcome is actually desired. I believe the `accel_per_hz` parameter was added to intentionally increase the sensibility at higher frequencies.
On the other hand, having 1/3 of the energy in higher harmonics might exerts more strain on the printer when measuring high frequencies. Although, it is possible that the springiness of the steppers take them out.
`TEST_RESONANCES` is sometime used with a very slow `hz_per_sec`, targeting a single frequency in order to pinpoint a resonance or find rattling parts. This use case may benefit from a cleaner test signal.
This could also prove handy if a future method for measuring resonances correlates the test signal frequency with the accelerometer output. Increasing the number of movement segments would yield an even cleaner test signal, albeit not entirely free of harmonics due to inherent non-linearities in the stepper and mechanical components.

In the end, itâs a bit of a moot point, but the modification is straightforward.

Iâm currently testing this to see if I see positive changes in the spectrograms.

2 Likes

Very interesting work. Thank you.

Years ago, the company I worked for (IBM/Celestica) had a very well equipped failure analysis and reliability testing testing lab with a lot of vibration (along with thermal) testing capabilities. Vibration testing was always done with single frequency sinusoidal vibration with no harmonics.

Iâve just deleted a long story about a customer I had in which we got them off of an âomnidirectional broadband vibration systemâ (almost literally a jack hammer welded to a metal plate) that produced random motion with varying g-forces with a single direction, single frequency sinusoidal vibration table and decreased their run in failure rate from 30% to less than 1% and increased their MTBF from two years to over twenty.

What I learned from this is that subjecting something to motions at higher frequencies and g-forces than a single frequency sweep, determined for the object under test, on a single axis results in:

• Failures that cannot be reliably determined to be at a specific frequency. It is very possible that the problem is outside the normal operating envelope of the object under test
• Immediate damage to the object under test
• Decreased MTBF (reliability) of the object under test

So, if you have a way to minimize the harmonics during resonance testing, then I would say letâs implement that!

FWIW, the standard specification for environmental testing of Airborne Equipment, RTCA DO-160, also describes vibration testing.
Quite complicate subject but one of the test procedures basically boil down to:

1. Find resonance frequencies with a sinusoidal acceleration input sweep across the spectrum. Whereas resonance frequency is defined as a peak of at least twice the input acceleration amplitude
2. Apply a random acceleration input curve across the spectrum to torture the equipment
3. Repeat 1. to check if the resonances have shifted (indication that something unwanted has happened)
4. Inspect for damages of the specimen

But I guess our intention is not to provoke failure in our printers and test for robustness but only find the resonance frequencies.

1 Like

Kind ping @dmbutyugin

Thanks, that is an interesting investigation. Generally speaking, I do not mind changing the form of the test signal. And the spectral leakage is a bit unfortunate, since it means that potentially we may need to put more energy into the printer to obtain the same result. That said, there are a few considerations, which may make it more troublesome than worthwhile.

First, the spectral leakage, even though it occurs, isnât that bad. If we look at the ideal frequency response (so, if the system followed the generated test precisely), the spectrogram would look like this:

Note that the coloring here is logarithmic. So indeed we observe multiple harmonics being produced. That said, not all of them are lost, for example, at 50 Hz response we get the energy of 10th harmonic from 5 Hz oscillations, energy of 9th harmonic from 5.5(5) Hz oscillations and so forth until 1st harmonic of 50 Hz oscillations. So, some energy isnât wasted. Now, if we look at the actual power spectral density, it looks like this:

Indeed, the ideal PSD grows linearly until the frequency of the end of the test, at which point we observe a large drop in the energy. In fact, if we look at it at the log scale and wider frequency range, it looks like

So, in fact the energy at the frequencies where we observe just the higher harmonics is almost 100 times smaller. So, @Piezo, Iâm unsure about your estimations of the higher frequency harmonics, perhaps the amplitude is like you calculated, but their energy density appears to be substantially smaller.

Another point is that you considered velocity, but this isnât whatâs being measured by the test. Instead, the test (well, accelerometer) measures acceleration. And that one doesnât have a triangle shape, it has a rectangular shape. Cutting velocity doesnât make acceleration look more âsine-wavyâ, but still keeps it rectangular, just increases its frequency (and introduces different spacing between the acceleration pulses).

And the third point is arguably quite a bit more important. It is that the segments of constant velocity (zero acceleration) are actually quite undesirable for 3D printing resonance testing. The reason is that you may hit one of the natural frequencies of a stepper motor(s) or the belts. Basically, motion with velocity `v` on GT2 belt with a 1.8 degree stepper motor and a pulley with `n` teeth will produce vibrations at the frequencies

``````v * 200 / (2 * n),
v * 200 / (2 * 2 * n),
v * 200 / (2 * 4 * n),
v / 2
``````

this will affect the measurements, and if you happen to hit a resonance frequency, this may substantially skew the results. Even if the motor will not move the full 1, 2 or 4 steps during the testing, I suspect the vibrations will still be produced. And I think it would be particularly difficult to compensate for that in the post-processing.

FWIW, it might be possible to solve this problem using either:

• Varying the cruising velocity throughout the test. Alas, I donât have a good idea how it should look like, e.g. should the cruising velocity increase or decrease throughout the test? And it should ideally change such that the vibrations excited at different speeds have the same spectral density.
• Instead of âapproximatingâ a sine wave by cutting the top of a velocity triangle, it may be beneficial to adjust the acceleration in a part of a test move in a more elaborate way. For example, if you target velocity, you may approximate the sine wave with 4 lines, having 2 different inclinations (i.e. acceleration near 0 will be higher than near the velocity maximum). I suspect you could still produce improved approximation that way, and avoid constant velocity cruising. If instead you decide to make acceleration profile more âsine-wavyâ, then you could split the current rectangular acceleration pulse into 3, with the middle having higher acceleration and the two side pulses having smaller acceleration, hopefully approximating sine wave better.

BTW, @Piezo, did you test the proposed change with different `coeff` values in `max_v = accel * t_seg * coeff`, and how did it look like?

1 Like

First, thank you for the feedback

Iâm surprised that we donât get the same ratio of power between the 1-st and 3-rd harmonics. For an ideal square wave the powers should be 8/ĎÂ˛ and 8/(9ĎÂ˛) respectively. Thatâs a 1/9 ratio, which match the 1/3 amplitude ratio reported by wikipedia. I measured a 1/71 ratio on your plot.
Can you upload the simulated raw accelerometer data? (or is the code available?)

This shows the power of each harmonic (normalized to the power of the fundamental):

The reduction are mostly about the 3-rd and 5-th harmonics. The 7-th, 15-th, 23-th, etc harmonics have slightly more relative power than the pure square-wave because of the normalization.

For a +/- 1 square wave, the fundamental harmonic amplitude is 4/Ď. The signal from the higher harmonic the square wave with cos(tau*t)*4/pi subtracted from it.

Iâm unsure which energy density calculation you are referring to. I defined the power/energy density as the auto-correlation over a period divided by the period. The +/- 1 square-wave have has a energy density of 1, and its first harmonic has a energy density of 8/ĎÂ˛ â 0.81. The higher harmonics represent a power of 1 - 8/ĎÂ˛ â 0.19.
Actually, Iâve used the square root of that in the SNR calculation (for no good reason?). This was reported as SNR = sqrt(8/ĎÂ˛ / (1 - 8/ĎÂ˛)) = sqrt(4.28) = 2.07

The optimized pseudo square-wave has a power of lambda â 0.742, and its first harmonic has a power of 0.684. This leave 0.057 remaining power for the higher harmonics. The 40% leakage reduction was reported at equal first harmonic amplitude: 40.1% * 0.19 / 0.81 â 0.057 / 0.684

For reference the Fourier coefficients of the modified square wave are:
With lambda=1 this yields the coefficients of the usual square wave, while lambda â 0.742 is the optimized version.

Hereâs my notebook with the calculations: gist

I havenât optimized the velocity profile per se. The method I took was optimizing the harmonic content of the acceleration while keeping the same trapezoidal moves; the visualization of the velocity is just an afterthought.
The spectra of all position derivatives are related, so making the velocity look more âsine-wavyâ also reduces the spectral leakage of the acceleration.
In fact, the modified square have has a better correlation with a sine wave, relative to its power density.

Iâm no entirely sure that this will be an issue. The cruised segments are very short, potentially shorter than the belt pitch or than an electrical rotation (4-steps). For the resonance to establish, few cycles are likely needed. With `accel_per_hz=75`, the cruised distance is 1.78/freq millimeters. Lower frequencies will be more affected, so yes this have to be tested for.

Currently the velocity is equal to `accel_per_hz/4`. We could randomly dither between `accel_per_hz/4` and `accel_per_hz * 0.742/4`. Or fade from `accel_per_hz/4` to `accel_per_hz * 0.742/4` as the 3-rd and 4-th harmonic go out of the frequency band of interest.

Yes, the sine wave can always be better approximated with more segments. Having an even number of segments per half-wave will avoid a cruised section.
However, Iâm not sure if this is practical at high frequencies. Some printers already struggle with the current profiles. Instantiating and issuing that many M204 commands might become problematic. I have a branch that adds a acceleration parameter to `toolhead.move` that could help with that.
Otherwise, this could be improved (and simplified) if we bypass the move planer and target directly the trapezoidal queue.

Not yet, sadly. My printer is in a dire condition amidst a toolhead re-design. Iâll will likely find the time to test next week end, but that will be on a sliding bed only.

Thank you both for the interesting story and references. This reinforce the idea that the test signal should be as clean as possible. It is still not clear to me in which extent the steppers (and maybe step compression) already band-limit the signal.

Hereâs the raw data:
raw_data_ideal.zip (1.3 MB)
You can pass it directly to `scripts/graph_accelerometer.py` (or `scripts/calibrate_shaper.py` if you will).
I was referring to Power Spectral Density (PSD) [link], which in our case is calculated normalized to the segment (of a fixed duration), so essentially to time.

Sorry, I didnât go through your evaluations yet, but maybe my data will give you something to look into.

I did the analysis of a single window using code mostly copied from klipper, and I do get the expected PSD ratio for the harmonics of 1/9, 1/25, etc: notebook@gist
Iâll have to dig a bit deeper to find what `graph_accelerometer.py` does differently.

Edit: It appears that the harmonic ratios you observed are an artifact produced by the time average of the spectral densities. This shows clearly when averaging over a varying number of windows:

The first windows is taken at the middle of the sweep, then more and more windows are added to the average in both direction.

My handwavy explanation is that higher harmonics ramp up faster in frequency than the fundamental. Therefore their energy is spread over more bins in the average.

I really appreciate the detailed response but I feel like there is a disconnect here with a basic assumption; the purpose of this test is, or rather should be, measuring the 3D printerâs response to a vibration input thatâs in the form of a simple sine wave. The goal of doing this is finding resonances that can affect the quality of the print and/or lead to the damage/premature wearing out of parts down the road.

I believe a properly set up system will not have any natural frequencies which cause resonances during normal operation and the output will be a fairly flat line of low amplitude across the vibration frequencies being tested.

The ideal vibration input should have a spectral signature of an impulse at the specified frequency with no harmonics. Including harmonics as paort of the sine wave, at any power level, will subject the printer to vibrations outside the specified range and could result in false positive readings, in the best case.

Now, I donât think anybody has done any studies as to what the response graphs actually mean, what can be done to reduce the spikes shown on them and what are the issues that arise if the spikes arenât flattened. I am interested in pursing this to see what resonance testing means in terms of affecting the quality of prints and longevity of the printer.

To be more precise, the ultimate purpose of this test is to find an input shaper configuration such that a printer can print without ringing/echo defects at the maximum acceleration and velocity.

Alas, thatâs not really what we observe in 3D printing, otherwise we wouldnât need input shaping. Most 3D printers happen to have resonances exactly in the operating range, so both the users and the test would be expected to hit some resonances under normal circumstances. FWIW, the transfer function of an oscillator is far from being flat (and 3D printer kinematics can be approximated by some combination of them).

Separately, thereâs really no âright rangeâ of frequencies for 3D printer. Ideally, we should identify all the resonances and for each one of them estimate how much they affect the print quality and compensate for that accordingly (if needed).

Thereâs actually a fair amount of literature discussing this topic. For example, previously I received a recommendation for the book âSystem Identification. A Frequency Domain Approachâ by Rik Pintelon,
Johan Schoukens, which describes several options, swept sine including. Well, we donât really use âsineâ due to practical reasons. But in principle there are other alternatives, e.g. Pseudorandom Binary Sequence.

Iâm afraid you are not doing the right justice to a lot of people who worked on the systems motion control research and input shaping theory over the last several decades. Unfortunately, Iâm nowhere near being an expert myself. FWIW, the effects of input shapers on the frequency response in terms of transfer functions are fairly well studied. Although the literature on input shaping for 3D printing is indeed rather limited, other adjacent areas of system motion control are studied in more depth. If you want to explore the topic more, I can recommend a few articles for starters (that I could find):

• Frequency weighted H2 optimization of multi-mode input shaper, M. Goubej, et. al.,
• Input Shaper Design in Convex Optimization Framework with Frequency Domain Constraints, H. S. Bae, J. C. Gerdes,
• An Adaptive Algorithm for the Tuning of Two Input Shaping Methods, M. Bodson,
• Basic algorithms of input shaping autotuning, M. Gniadek, S. Brock,
• A limited-preview filtered B-spline approach to tracking control - With application to vibration-induced error compensation of a 3D printer, M. Duan, D. Yoon, C. E. Okwudire,
• Effects of input shaping on two-dimensional trajectory following, W. E. Singhose, N. Singer,
• Many excellent articles by W. E. Singhose et. al. on various input shapers that are now standard.

And of course any further research in this area will be highly encouraged and welcome!

OK, I see where the differences come from, but Iâm afraid I donât fully follow. The averaging is applied in the same manner to all frequencies, therefore it should show the âtrue average power spectral densityâ? And Welchâs algorithm to estimate PSD does exactly that. Basically, if we were performing vibrations at constant frequency, then the results will be like youâve estimated, but since weâre performing a sweep, we get what we get instead. Am I wrong here?

Yes, this is what I meant. I should have specified more clearly that I analyzed the spectral content of the base wave form and not the chirp. I canât do this kind of math and frankly I donât think it is needed to show that reducing the harmonic content of the non modulated waveform will reduce the leakage over the whole test.

The total or average energy dumped at a given high frequency might be not that high, but at some point in time this band will see the full instantaneous power of some harmonic which is as I described.
My proposal is reduce this with a simple change in the current framework. The reduction will be transcribed proportionally in the power density average over the whole test, which is also worthwhile.

As a side node, the Welchâs algorithm is based on the short time Fourier transforms. It is dependent on the time-frequency compromise of the specific Gabor frame (nfft, noverlap). The Fourier transform on the whole `raw_data_ideal` signal shows a different picture:

I was trying to generate a band-limited version of the ideal signal, as I think we might get better estimations without the aliased components (the quilted pattern on your spectrogram). But it turns out that doing BLIT synthesis of a FM signal is not as easy as I initially thought.
I should probably direct my energy into real world testing instead, and investigate the âcruising resonancesâ that you mentioned.

Thank you for the comprehensive reply. I have a couple of comments back:

Can you explain why velocity is part of the discussion? As I understand it, velocity would be a variable dependent on maximum acceleration rather than something we test for on its own.

I understand this, but we should be able to get quite accurate sine wave simulations, no? It is my understanding that was what @Piezo 's first post was about.

I meant to say when it come to 3D printers (I thought that would be obvious from the context of the discussion) - as far as I can tell nobody has done any research on analyzing and applying the results from resonance testing on 3D printers themselves and what it means in terms of print quality and printer longevity.

I looked at the abstracts of the articles you provided in your reply but I donât see any that discuss the effects on 3D printers and how to minimize them in hardware (ie mass redistribution, stiffening and dampening parts) as opposed to software.

This is the optimal version I got when I tried that:

The square root SNR (to be consistent with the values in the OP) is 4.79.

• blue: the modified square wave with reduced acceleration segments,
• green: the fundamental harmonic
• orange: the next 10 harmonics

This plot shows the energy of harmonics, relative to the energy of the fundamental.

• blue: pure square wave
• orange: adding cruised segments as proposed in the OP
• green: â3-steps square-waveâ with reduced acceleration/deceleration segment, as shown above.

That plot illustrate exactly what is optimized: I minimized the sum of energy of the harmonics relative to the fundamental. `argmin (power - fundamental_power) / fundamental_power`.

The optimal duration and acceleration values for the added segment are again non trivial. I couldnât find a closed form solution for the duration of the acceleration. It is possible that itâs related to the Dottie number in some way. The reduced acceleration value is not simply the mean of the cosine over the segment.
This is the value I found:

• duration: `0.3905 * T/4` (fraction of a fourth of a period),
• acceleration: `0.3480 * accel` for a +/- `accel` waveform.

Now in practice, this requires issuing 4 `M204`s per period. My concern is that the current code will not be performant enough for frequent acceleration changes at the end of the frequency range.

1 Like

I applied few corrections to the OP:

• To output the correct frequency, the move distance must be adjusted by a factor `Îť (2 - Îť) â 0.9334`. Big mistake.
• To have the same fundamental harmonic amplitude, `accel_per_hz` should be adjusted by a factor of `1 / sin(Ď Îť / 2) â 1.0881`. Minor correction for easier comparison of the results. It might be better to change the default from 75 to 82 mm/sÂ˛/Hz than lying on the applied acceleration.

I did some preliminary testing on a narrow frequency band with `TEST_RESONANCES AXIS=X OUTPUT=raw_data FREQ_START=40 FREQ_END=40.5 HZ_PER_SEC=0.01`

I spied on the plannerâs outputs to check the modifications. These are the two moves composing a half-period, which are repeated in a loop:
Unmodified Klipper:

``````accel:6.2ms cruise:0.0ms decel:0.0ms startv:0.0 cruisev:18.8 accel:3037
accel:0.0ms cruise:0.0ms decel:6.2ms startv:18.8 cruisev:18.8 accel:3037
``````

Modified for cruised segments:

``````accel:4.6ms cruise:1.6ms decel:0.0ms startv:0.0 cruisev:15.1 accel:3305
accel:0.0ms cruise:1.6ms decel:4.6ms startv:15.1 cruisev:15.1 accel:3305
``````

Comparing the real-life spectra:

Blue is the unmodified pure square wave, orange is the wave-form with cruising introduced.

• The correction factors for the fundamental energy and the frequency work as expected.
• The even harmonics (2-nd peak @ 80 Hz, 4-th @ 160 Hz, etc) are not supposed to exist in a square-wave. I blame non linearities in the steppers and possibly the motion system (belts). Interestingly, most of the even harmonics are amplified with the modified square wave. This could be possibly be explained by resonances accumulated during the cruised segments, but I continue to doubt that this can happen in 1.6ms.
• The odd harmonics are reduced, but not in the proportions of the harmonic series.

Things Iâll experiment with:

• Perform the test centered at different locations in the stepperâs electrical cycle. I expect less even harmonics if the deviation is odd-symmetric on both sides of the test point. The test above was executed at a central position of 1.898 step (modulo 4 steps).
• Align the cruised velocity with the velocity at which my steppers resonance resonates the most and the least. The maximum is arround ~60 mm/s on my printer, which is reached when setting `accel_per_hz = 404`.
• Implement the 3-steps 2 accelerations wave-form from my previous post, as suggested by @dmbutyugin.
• Make n balance sheet of the energy spent in and out the frequency band of interest over the whole frequency sweep.
1 Like

Interesting. Thank you for sharing this.

Well, if someone prints at 20-30 mm/sec, they could afford setting much higher acceleration compared to printing at 100-150 mm/sec, thatâs all I meant.

We could generate a sine wave of better shape, or even a perfect sine wave, just thereâs certain limit at which it may no longer be practical. Separately, I meant that there are other ways to generate a good test signal that do not cause substantial spectral leakage and that do not involve sine waves generate.

Ah, yes, fair enough, itâs a misunderstanding on my part then. Yes, if you want to look into that particular aspect - it would be very interesting to get some results.

1 Like

Thanks, the results are very interesting.

I didnât check your math, but I think both changes are conceptually OK if youâve confirmed that they do whatâs expected of them (and you did).

I have a slight suspicion that there may be another reason for that. Basically, with some constant cruising velocity the test should generate the acceleration signal something like: (+a, 0, +a, -a, 0, -a, âŚ). I suspect this could cause even harmonics to actually amplify. If you want to, you could post your patch to the mainline Klipper code here as a diff (I suspect itâs very small), then I could generate an âidealâ acceleration signal for you for further debugging.

Your proposal here makes a lot of sense. I would further suggest you to not test just ~60mm/sec cruising velocity, but also some other speeds (e.g. where you may have significant VFA, if you know the corresponding range of speeds on your printer).

Analytically, I didnât find any even harmonic. Both waveforms have odd symmetries at +/- T/4, where the acceleration sign flips, or at the middle of the cruise segment. This cancel out the inner product with sine bases at `2*n*frequency`.

You can find the patch in the first post of this thread. This should be enough to replicate the wave-form. A simulated acceleration signal would be very welcome, thanks.

To experiment with the effect of the stepper position, I modified the `resonance_tester` code to sweep positions instead of frequencies. At each iteration a very small offset (83.33 nm) is added to the backstroke. Over the course of the 6 min long test, the toolhead travels 3 electrical cycles, 12 steps, or 2.4 mm.

(I tweaked the visualization with `vmin, vmax = np.quantile(pdata, [0.25, 0.999])` and `noverlap=M - M//4`.)

• The second and fourth harmonics go in and out twice per step.
• The 3rd harmonic, which is present in the ideal signal, also shows some modulation, but once per step.
• Overall the positions with the least amount of harmonics seems to be half-steps.
• The electrical cycle is marked by these wide band noise patterns at 150-400 Hz. Otherwise, there is no other pattern that is different on each of the 4 steps, as suggested here. But these results might be specific to this test (and 40Hz frequency).
• The 6-10th harmonics seems to show faster paced patterns, but itâs quite faint and low resolution. I could try to stretch the test even more in time.

The spectrogram above used the modified wave form. This one is using the original square-wave:

Very similar beside the stronger 3rd harmonic.

Edit: I repeated the test above but over 4 steps for a duration of 10 min. I plotted the power of the first 5 harmonics over the central position of the oscillation:

It appears that the pattern on the 2nd harmonic has a period of 2 steps.

Iâm not sure where Iâm going with this. I guess this shows that harmonics are generated by the stepper themself, so itâs less relevant to try to remove them from the test signal. On the other hand, there is less effect on even harmonics, likely because they are already present in the test signal, showing room for improvement.

@dmbutyugin Hi. Did you get a chance do this? (The patch is in the first post). Otherwise can you share the code/method to generate the waveform?