Sine hard sync

[rigatoni_modular]

Coming across this thread years later... I'm looking at the blep2 line on the graph and wondering why I'm not getting the same result. I'm trying to get there by following this part of the linked Native Instruments paper:

Screenshot 2023-11-06 at 5.20.59 PM.png

However, in Desmos when I use a power series approximation of the integral of Sinc to get the blep and try to follow the equation in that part of the paper to get the next antiderivative blep I just get something that gets large pretty quickly on either side of the x axis and it's not obvious to me what I'd need to subtract from it to get that nice zero-centered blep from xhy3's graph. Note that I've excluded the sin/cos term in that equation because I'm not even in the right ballpark.

Screenshot 2023-11-06 at 5.17.57 PM.png

[2DaT]

Coming across this thread years later... I'm looking at the blep2 line on the graph and wondering why I'm not getting the same result. I'm trying to get there by following this part of the linked Native Instruments paper:

Screenshot 2023-11-06 at 5.20.59 PM.png

However, in Desmos when I use a power series approximation of the integral of Sinc to get the blep and try to follow the equation in that part of the paper to get the next antiderivative blep I just get something that gets large pretty quickly on either side of the x axis and it's not obvious to me what I'd need to subtract from it to get that nice zero-centered blep from xhy3's graph. Note that I've excluded the sin/cos term in that equation because I'm not even in the right ballpark.

Screenshot 2023-11-06 at 5.17.57 PM.png

The recursive formula works for "blep residuals", that means we need to subtract the "naive" step from it after we compute it.

Here is corrected version that computes analytic blep residuals.

https://www.desmos.com/calculator/j3jzbjc10m

Bleps 1, 2 and 3 look about right, but I am not 100% sure about oscillating nature of higher bleps.

Some things to notice: if you want to tabulate these, you need to use windowing and the higher the blep, the more smooth windowing it needs. I don't have the exact recipe how to window analytic bleps, but in theory it shouldn't be too hard to have something workable with some experimenting.

Windowing should counteract the infinite growth of higher order blep residuals too.

Also do notice that first blep residual is non-continuous at midpoint. That means that you can't linearly interpolate the table through 0. But it's possible to sign-invert the right half of the blep, interpolate through that and negate the sign in the blep insertion code.

[mystran]

Some things to notice: if you want to tabulate these, you need to use windowing and the higher the blep, the more smooth windowing it needs. I don't have the exact recipe how to window analytic bleps, but in theory it shouldn't be too hard to have something workable with some experimenting.

I have never managed to get windowing of analytic BLEPs to work acceptably even though I've tried... and I'd actually argue that it's "sort of wrong" from the theoretical point of view 'cos you should really integrate w(x)*sin(pi*x)/(pi*x) where w(x) is the window in order to reconstruct "what would have been produced had we filtered in continuous time." If you think of what the target spectrum for higher order BLEPs looks like, poor results from windowing those is not entirely unexpected either.

Now, perhaps windowing analytic BLEPs can work... but I have not been able to get this to work. What I have been able to do is window a sinc and numerically integrate the result up to four times (with good quality), which is sufficient to anti-alias piece-wise cubic waveforms. I'd imagine more is possible if the prototype filter is sufficiently flat around DC, but it becomes increasingly difficult to get that for higher orders with windowing. Cubics are doable basically just by bumping Kaiser alpha as long as the kernel is long enough... and they've been kinda sufficient for my purposes.

Beyond cubics, I've not managed to these to behave, except for kernels (eg. Lagrange) that are of no real practical value... but I think it's doable in theory.

[2DaT]

Now, perhaps windowing analytic BLEPs can work... but I have not been able to get this to work. What I have been able to do is window a sinc and numerically integrate the result up to four times (with good quality), which is sufficient to anti-alias piece-wise cubic waveforms. I'd imagine more is possible if the prototype filter is sufficiently flat around DC, but it becomes increasingly difficult to get that for higher orders with windowing. Cubics are doable basically just by bumping Kaiser alpha as long as the kernel is long enough... and they've been kinda sufficient for my purposes.

Beyond cubics, I've not managed to these to behave, except for kernels (eg. Lagrange) that are of no real practical value... but I think it's doable in theory.

Well, I didn't do analytic BLEPs either.

It should be possible to at least alleviate some issues by using multiprecision library (such as boost::multiprecision) and use a very low step for numerical integration. We don't need very high oversampling factor for the table itself, but it will improve the accuracy of integration.

[Z1202]

I have never managed to get windowing of analytic BLEPs to work acceptably even though I've tried... and I'd actually argue that it's "sort of wrong" from the theoretical point of view 'cos you should really integrate w(x)*sin(pi*x)/(pi*x) where w(x) is the window in order to reconstruct "what would have been produced had we filtered in continuous time." If you think of what the target spectrum for higher order BLEPs looks like, poor results from windowing those is not entirely unexpected either.

I'm not sure whether from a theoretical standpoint either of integrating post- or pre-windowing is more wrong than the other. Windowing gives you approximate spectrum anyway, so both are kinda wrong, and I'm not sure whether one option of approximating before or after the integration is more theoretically wrong than the other.
It's possible that standard window functions are not too suitable for windowing higher-order BLEPs (that is post-integration) merely because of their slow spectral rolloff speed. E.g. Kaiser window's spectrum rolls off at 6dB/Oct, which might kinda defeat the purpose of higher-order BLEPs which are supposed to adjust the spectrum to make sure that it rolls off much faster than that. Just a theory.

[mystran]

It should be possible to at least alleviate some issues by using multiprecision library (such as boost::multiprecision) and use a very low step for numerical integration. We don't need very high oversampling factor for the table itself, but it will improve the accuracy of integration.

It's not an integration issue. It is literally an issue with the prototype kernel. If it is not flat enough around DC (ie. enough derivatives are at least approximately zero), the BLEPs go wonky independent of any numerical integration issues.

Numerical integration precision is something you DO need to pay attention to, but .. if you just integrate the first half and use symmetries, it's not too bad. The real problem is the prototype kernel.

[mystran]

]I'm not sure whether from a theoretical standpoint either of integrating post- or pre-windowing is more wrong than the other. Windowing gives you approximate spectrum anyway, so both are kinda wrong, and I'm not sure whether one option of approximating before or after the integration is more theoretically wrong than the other.

Well.. it depends on what you take as the theoretically "correct" result, but the original idea on which BLEPs are based on is that if we filter a waveform that is a sum of a linear function and some Heaviside steps, then under the assumption that our filter kernel is identity with respect to the linear function, we can only compute the filter "residue" with respect to the steps.

So.. in a sense, we are directly reconstructing "what would have resulted" had we filtered the continuous-time waveform using our original sinc-filter. As the Heaviside step is the integration operator, our BLEP is the integral of our chosen kernel (ie. the windowed sinc, typically). Similarly for the convolution of two Heaviside steps, which gives us a discontinuity in 1st derivative, which we can reconstruct using the second integral of our chosen kernel.

So it is not that about whether or not the kernel itself is an approximation, rather it's about the fact that we want to reconstruct "what would have resulted" had we filtered in continuous-time with whatever approximate kernel we might have.

This view into the whole algorithm also directly explains why higher orders are increasingly difficult: the filter must be identity with respect to higher degree polynomials, or put another way, the filter response around DC must be flat to a higher order. You can easily verify this empirically, just construct Lagrange interpolators which by definition maximize the flatness around DC; the order will directly (and exactly) predict how many times you can integrate the result (and you can do this analytically, because they are piece-wise polynomial) before the construction falls apart.

Unfortunately said Lagrange interpolators are not practically terribly useful for the purpose of constructing anti-aliased oscillators as anyone who is willing to try will quickly discover.

ps. I would like to add that I do not have fully satisfactory, mathematically rigorous proof that having zero derivatives around DC in spectral domain translates into identify with respect to polynomials in time-domain.. but this seems "reasonable" in terms of treating polynomials as integrals of DC signal and empirically it appears to be the correct condition.

[rigatoni_modular]

Coming across this thread years later... I'm looking at the blep2 line on the graph and wondering why I'm not getting the same result. I'm trying to get there by following this part of the linked Native Instruments paper:

Screenshot 2023-11-06 at 5.20.59 PM.png

However, in Desmos when I use a power series approximation of the integral of Sinc to get the blep and try to follow the equation in that part of the paper to get the next antiderivative blep I just get something that gets large pretty quickly on either side of the x axis and it's not obvious to me what I'd need to subtract from it to get that nice zero-centered blep from xhy3's graph. Note that I've excluded the sin/cos term in that equation because I'm not even in the right ballpark.

Screenshot 2023-11-06 at 5.17.57 PM.png

The recursive formula works for "blep residuals", that means we need to subtract the "naive" step from it after we compute it.

Here is corrected version that computes analytic blep residuals.

https://www.desmos.com/calculator/j3jzbjc10m

Bleps 1, 2 and 3 look about right, but I am not 100% sure about oscillating nature of higher bleps.

Some things to notice: if you want to tabulate these, you need to use windowing and the higher the blep, the more smooth windowing it needs. I don't have the exact recipe how to window analytic bleps, but in theory it shouldn't be too hard to have something workable with some experimenting.

Windowing should counteract the infinite growth of higher order blep residuals too.

Also do notice that first blep residual is non-continuous at midpoint. That means that you can't linearly interpolate the table through 0. But it's possible to sign-invert the right half of the blep, interpolate through that and negate the sign in the blep insertion code.

Wow, thank you! I was barely even expecting a revival of this 3 year old topic

Those residuals look good. I did end up realizing last night that I was approaching it wrong by working off the band-limited step rather than the delta h_0 residual.

A couple things I'm still curious about:
- The left side of the equation in that part of the paper multiplies the residual by `n`, so it seems like each residual n >= 2 that we have in Desmos should have the whole right side divided by `n`. Is that reasonable? Seems like it also will end up getting us the scaling factor that xhy3 was looking for.
- To use utilities I already have to compute a minimum phase version of these higher-order bleps, would I need to compute the minimum phase response for (trivial nth degree step - nth degree residual) rather than the residuals themselves? I like the convenience of just being able to mix in the residual but it's not super viable when users of my oscillator could sync it at any arbitrary sample and I can't use lookahead to center the linear phase residual

[2DaT]

- To use utilities I already have to compute a minimum phase version of these higher-order bleps, would I need to compute the minimum phase response for (trivial nth degree step - nth degree residual) rather than the residuals themselves? I like the convenience of just being able to mix in the residual but it's not super viable when users of my oscillator could sync it at any arbitrary sample and I can't use lookahead to center the linear phase residual

Minimum phase doesn't really work, especially with hard sync and/or high order bleps. It's hard to center the transition, which for min-phased isn't really obvious where the naive step really happens.

Lookahead means latency and realistically you don't need more than 16 samples of latency to get waveforms with next to 0 perceptible aliasing. It will prevent low-latency feedback loops from functioning properly, but there isn't much that can be done, apart from tricks and hacks, but that's whole another discussion.

As for the scaling, it depends on what kind of scaling you use for phasors. I prefer [0,1) and it's simpler to calculate fractions and offsets in, but that requires the usage of sin(pi*x)/pi*x in blep calculations.

[rigatoni_modular]

So for an oscillator with hard sync and linear phase bleps you'd just start mixing in the bleps immediately but delay the actual sync of the generated waveform until the transition in the blep?

Also, is the reason min-phase bleps of order higher than 1 are problematic because you have to take the alternate approach of doing the min phase conversion on the integrated bandlimited step itself, which grows really quickly as it gets further away from zero unlike the 1st order blep?

[2DaT]

So for an oscillator with hard sync and linear phase bleps you'd just start mixing in the bleps immediately but delay the actual sync of the generated waveform until the transition in the blep?

Naive waveform (without any kind of antialiasing) should be delayed by half of the blep kernel length, so that when summed, the naive transition aligns with the exact middle of the blep residual (with sample and subsample precision). You do need a some kind of ring buffer for blep residuals and naive delay. It's certainly possible to do the whole thing with one ring buffer with 2 read positions, but that's a little bit more complicated.

[rigatoni_modular]

Oof, feels like kind of a pain to delay everything since this is for modular synthesis and there's plenty of feedback stuff going on there. Maybe I'll just try numerically integrating the minBLEP a few times and call it a day :/

[Z1202]

So it is not that about whether or not the kernel itself is an approximation, rather it's about the fact that we want to reconstruct "what would have resulted" had we filtered in continuous-time with whatever approximate kernel we might have.

I'm not sure I'm following the argument. Once there is some approximation step artificially inserted into your process, your result is "wrong", isn't it? As for the errors at DC, they are a bit braindamaging to reason about, IMHO, but in principle it seems to me that they should be inaudible anyway, as, worst case, they should result in a single fully continuous polynomial component (which is somehow offset by the remaining part of the function). But I'm even not 100% sure they are not 100% cancelling out automatically anyway if you use BLEP residuals, rather than explicit integration with windowed sinc. It might be even the other way around - direct integration with sinc might be more prone to polynomial offsets than using BLEP residuals. Nothing too conclusive here in my ideas though, it's long time since I properly dived into this topic.

[mystran]

So it is not that about whether or not the kernel itself is an approximation, rather it's about the fact that we want to reconstruct "what would have resulted" had we filtered in continuous-time with whatever approximate kernel we might have.

I'm not sure I'm following the argument. Once there is some approximation step artificially inserted into your process, your result is "wrong", isn't it?

Let me try to put it in slightly more formal terms. We are trying to sample the convolution of a continuous-time waveform with a low-pass filter. Convolution and integration (and we can express integration as a convolution) linear and associative, so we can reassociate without changing the results.

Now, under the assumptions that (1) the waveform we are trying to convolve is piece-wise polynomial of degree up to N and (2) our chosen low-pass filter is the identity operator with respect to polynomials of degree up to M>=N, the BLEP algorithm implements a linear convolution in continuous-time and then samples the result. We deconstruct the waveform into various integrals of Heaviside, convolve each part separately, swap the integrals with the convolutions (which is fine, they are linear) and drop the identity terms... then add it all back together and sample the result.

Note that up to this point I have only said "low-pass filter" and this is fundamental, because the way you design the low-pass filter and what it's response looks like is not important as long as it satisfies assumption (2) above to a good approximation. Strictly speaking our filter does not even need to have a low-pass response as far as the BLEP-convolution itself goes, that's only needed in the final sampling step to prevent aliasing.

Now.. the way you actually design the low-pass filter is up to you, it's effectively a separate thing from the BLEP-convolution algorithm that just needs some suitable kernel... but if you do choose w(x)*sinc(x) for some window w as your low-pass filter, then the multiplication with a window is not linear, it does not reassociate directly with integration, so if you distribute it after the reassociation steps of the BLEP algorithm then your algorithm is no longer computing a linear convolution! As far as I'm concerned, it is no longer entirely clear exactly what are we trying to compute.

The fact that numerically our (linear) BLEP algorithm will only reconstruct "approximately" correct linear convolution means that there's a noise floor, but with a correct implementation you can mostly have that noise floor below the aliasing attenuation of a reasonable kernel, so not a big deal.

[mystran]

Now.. the case where the waveform we are trying to sync is not a polynomial: what we are effectively doing is substituting a Taylor-series expansion on both sides of the hard-sync point and then anti-aliasing the transition between the Taylor polynomials. We can still do linear convolution here, it's just that the error-term of the Taylor-expansion will be left there as aliasing.

You can actually get rid of the aliasing "completely" if rather than doing on-demand Taylor expansions, you approximate the sine itself with a piece-wise polynomial (eg. piece-wise cubic) waveform. The downside is that you'll need to do more BLEPs now, because it's a BLEP waveform now, but now whatever errors are left are just harmonics and BLEPs will filter them out once they reach the kernel band-limit.