What is the best way to reduce aliasing in this oscillator?

[alestsurko]

Hi everyone. First of all sorry for may be some dumbest in this question. I'm a beginner in DSP and my math skills is far from perfect. But I'm learning. I found here and in another places a lot about anti-aliasing methods. I'm not fully understand how to implement it in all cases, but I understand how it works. I'm trying to make something like a sine oscillator, but from two parabolas. And I also want to modulate the parabolas width. Here is Julia code for trivial implementation of this oscillator:

Code: Select all

function parWM(w, c::PhaseCounter)
  phase = phasor(c)

  if phase < w
    return -4 * (phase/w - 0.5)^2 + 1
  else
    return 4 * ((phase-w)/(1-w) - 0.5)^2 - 1
  end

end

phase is a phase counter and 0 ≤ phase ≤ 1;
w is width and 0 ≤ w ≤ 1.

width = 0


width = 0.3


width = 0.5


width = 0.7


width = 1


I'm looking for efficient method to reduce aliasing in this oscillator. Not only looking, but want to understand what is better and why. And I hope for your help... On this forum I found polyBLAMP implementation by Tale:

http://www.kvraudio.com/forum/viewtopic ... r#p6107378

And it's working well in my case when width = 1 or width = 0 (when it's just a parabola). It also suppress aliasing with other width value, but there is audible aliasing.

Also I tried to apply polyBLAMP from this paper:

http://www.ness-music.eu/wp-content/upl ... squeda.pdf

But I failed... (I almost sure, that I did something wrong).

I'd not write, but I'm in stuck. The only solution currently I see is using Tale's polyBLAMP with oversampling. But I'm not sure this IS a solution and I believe there is a much better way. I'm also planning to use width in range [0.05, 0.5] or so, but currently it's not matter, I think. Also I'm trying to suppress aliasing in a single-parabola version of this oscillator for start:

Code: Select all

function parWM_bi_osc_lf(w, c::PhaseCounter)
  phase = phasor(c)

  if phase < w
    return -4 * (phase/w - 0.5)^2 + 1
  else
    return 0
  end

end

[logicalhippo]

Since the oscillator is quite complicated, oversampling might be your best bet.

Another approach would be to represent your signal as a doubly-integrated asymmetric pulse wave, which you could implement with BLEP and then numerically integrate (you'd probably need a high-pass filter to prevent errors from accumulating). I'm not sure how well this would work.

[JCJR]

Hi alest

I do not know your solution. Merely offering a Captain Obvious observation of which you may already be aware.

If no aliasing occurs at w = 0 and w = 1, then your antialias code is working to smooth the nonlinearity of the wraparound of phase thru 1 back to 0 at the end of each wave cycle.

At other values of w, you get a less-intense nonlinearity somewhere in the middle of each wave cycle. In your images it is the sudden slope transition at the mid wave zero crossings, when the phase makes transition from phase < w to phase > w.

If you can apply your smoothing method to those midpoint nonlinearities, in addition to smoothing the phase wraparound nonlinearities, maybe it would be well anti aliased for all values of w. That is a very pretty looking waveshape modulation. And requires not much math from the cpu, at least in the non antialiased form!

Does it also alias if w = 0.5? If your remaining aliasing comes from zero crossing nonlinearities, then w = 0.5 shows a smooth zero crossing and might be expected to sound about as clean as w = 0 or w = 1?

[alestsurko]

Thanks guys! I already thought that nobody would answer...

Since the oscillator is quite complicated, oversampling might be your best bet.

Another approach would be to represent your signal as a doubly-integrated asymmetric pulse wave, which you could implement with BLEP and then numerically integrate (you'd probably need a high-pass filter to prevent errors from accumulating). I'm not sure how well this would work.

Well, after all my searches, I decide to do oversampling. Your words added me a determination

If no aliasing occurs at w = 0 and w = 1, then your antialias code is working to smooth the nonlinearity of the wraparound of phase thru 1 back to 0 at the end of each wave cycle.

At other values of w, you get a less-intense nonlinearity somewhere in the middle of each wave cycle. In your images it is the sudden slope transition at the mid wave zero crossings, when the phase makes transition from phase < w to phase > w.

If you can apply your smoothing method to those midpoint nonlinearities, in addition to smoothing the phase wraparound nonlinearities, maybe it would be well anti aliased for all values of w. That is a very pretty looking waveshape modulation. And requires not much math from the cpu, at least in the non antialiased form!

Does it also alias if w = 0.5? If your remaining aliasing comes from zero crossing nonlinearities, then w = 0.5 shows a smooth zero crossing and might be expected to sound about as clean as w = 0 or w = 1?

Actually the Tale's polyBLAMP is working quite well for a parabola (like my oscillator with w = 0 and w = 1). So I do polyBLAMP two times - for two parts of the oscillator. In Julia it looks like this:

Code: Select all

function polyBLAMP(t, dt)
  if t < dt
    x = t/dt - 1
    return -1/3 * x^3
  elseif t > (1 - dt)
    x = (t - 1)/dt + 1
    return 1/3 * x^3
  else
    return 0
  end
end

function parWM_bi_osc(w, c::PhaseCounter)

  phase = phasor(c)

  if phase < w
    x = -4 * (phase/w - 0.5)^2 + 1
  else
    x = 4 * ((phase-w)/(1-w) - 0.5)^2 - 1
  end

  scale_a = 2 * c.incr/w
  p1 = scale_a * polyBLAMP(phase, c.incr)
  p2 = scale_a * polyBLAMP(mod(phase + 1 - w, 1), c.incr)

  scale_b = 2 * c.incr/(1-w)
  p3 = scale_b * polyBLAMP(mod(phase, 1), c.incr)
  p4 = scale_b * polyBLAMP(mod(phase - w, 1), c.incr)

  return x + p1 + p2 - p3 - p4
end

But while it works for a parabola, it has much more aliasing with a width modulation version (but less when compared with a trivial implementation):

Code: Select all

function parWM_bi_osc_lf(w, c::PhaseCounter)
  phase = phasor(c)

  if phase < w
    return -4 * (phase/w - 0.5)^2 + 1
  else
    return 0
  end

end

And, as I understand it, the cause is a different ramp form in the last case. PolyBLAMP for such a ramp described in the "ALIASING REDUCTION IN SOFT-CLIPPING ALGORITHMS" (I gave a link above). The difference with the Tale's polyBLAMP is in the second polynomial. But likely because of my math skills I had no success with implementation of the polyBLAMP from that paper

Anyway today I did a simple "real-world" implementation and it's not so un-efficient as I thought it would be. Especially that the synth will be monophonic. So for now I decide to do oversampling. Maybe in the future, when I improve my math skills, I'll find a better solution.

[Mayae]

Can you make some example sounds with it? Would like to hear how the oscillator sounds

[JCJR]

Does the simple antialiased version even sound bad at w = 0.5? I'm not an expert, merely curious.

[alestsurko]

Hi guys! Sorry for the late reply...

Does the simple antialiased version even sound bad at w = 0.5? I'm not an expert, merely curious.

Sorry that I didn't answer in the first time. Yes, I remember that there was slightly more aliasing at w=0.5 comparing to trivial implementation. But I've done so much versions... As I remember that was when I did polyBLAMP only once, but I'm not sure.

Can you make some example sounds with it? Would like to hear how the oscillator sounds

Sure:

https://dl.dropboxusercontent.com/u/286 ... ample.aifc

Currently I'm using polyBLAMP and 4x oversampling. Performance is quite well on my mac, but on iPad 4th generation it's about 20% of the CPU usage for the only one oscillator! I've planned two... and filter... and reverb... and chorus... and some saturation...

[Mayae]

Interesting. The CPU usage sounds quite high, I've done way more in C# on older android phones? Are you using Julia still?

[alestsurko]

Interesting. The CPU usage sounds quite high, I've done way more in C# on older android phones? Are you using Julia still?

No, the DSP part is in C++. In iOS simulator (on mac) it's only 2-3%.

[alestsurko]

I've improved performance to 5-6% of the CPU usage. But comparing to usual mac's performance it's still too high. Or is it normal for the iOS/mobile world?

[mystran]

Assuming a "slow" modulation rate for the "pulse-width" (ie. slow enough that aliasing due to side-bands generated by the modulation itself can be neglected and every cycle can be treated as stationary; this is the usual assumption built into BLEP-based methods), your highest order discontinuity is in 2nd derivative, hence 2nd order BLEP ("BLORA"?) can give you the "exact" solution without oversampling.

More generally, given a piece-wise polynomial wave-form with Nth order segments, you will need to handle at most N derivatives (+ the zero-order step discontinuity), except for special cases where some of the derivatives happen to match (as would be the case for a saw-waves 1st derivative and your "width=0" parabola's 2nd derivative, etc).

Various methods for generating such "higher-order BLEPs" have been described elsewhere (by myself and others), but assuming you have a table of such kernels, the way you would anti-alias any transition from one polynomial segment to the next is fairly simple: first evaluate both polynomials at the point of discontinuity to get the step-discontinuity (as the difference of the two values) and insert appropriate relevant "0-order" BLEP.

Then take the derivatives of both polynomials and repeat the evaluation of difference, except this time you correct it using a "1st order" BLEP (or BLAMP as some people would call it). Repeat this process, taking progressively higher order derivatives and correcting them with higher order BLEPs. If the polynomials (or the one with higher order) are of order N, then after Nth derivative, all the remaining derivatives are necessarily zero (and hence the difference between the derivatives of the two polynomials is also zero) and you have correctly compensated for the discontinuity.

In theory, this process works for arbitrary orders. In practice, it appears to be the case that higher-order BLEPs place increasing demands on both the prototype filter and the numerical precision of the implementation (in addition to the increase in work as you have to solve and correct more derivatives), so whether this approach can be effective for short PolyBLEPs I honestly don't know. With "long" kernels, it does work reasonably well at least up to cubics though.

In any case, hopefully this will give you some additional insight into the results you have been able to obtain.

[alestsurko]

Assuming a "slow" modulation rate for the "pulse-width" (ie. slow enough that aliasing due to side-bands generated by the modulation itself can be neglected and every cycle can be treated as stationary; this is the usual assumption built into BLEP-based methods), your highest order discontinuity is in 2nd derivative, hence 2nd order BLEP ("BLORA"?) can give you the "exact" solution without oversampling.

More generally, given a piece-wise polynomial wave-form with Nth order segments, you will need to handle at most N derivatives (+ the zero-order step discontinuity), except for special cases where some of the derivatives happen to match (as would be the case for a saw-waves 1st derivative and your "width=0" parabola's 2nd derivative, etc).

Various methods for generating such "higher-order BLEPs" have been described elsewhere (by myself and others), but assuming you have a table of such kernels, the way you would anti-alias any transition from one polynomial segment to the next is fairly simple: first evaluate both polynomials at the point of discontinuity to get the step-discontinuity (as the difference of the two values) and insert appropriate relevant "0-order" BLEP.

Then take the derivatives of both polynomials and repeat the evaluation of difference, except this time you correct it using a "1st order" BLEP (or BLAMP as some people would call it). Repeat this process, taking progressively higher order derivatives and correcting them with higher order BLEPs. If the polynomials (or the one with higher order) are of order N, then after Nth derivative, all the remaining derivatives are necessarily zero (and hence the difference between the derivatives of the two polynomials is also zero) and you have correctly compensated for the discontinuity.

In theory, this process works for arbitrary orders. In practice, it appears to be the case that higher-order BLEPs place increasing demands on both the prototype filter and the numerical precision of the implementation (in addition to the increase in work as you have to solve and correct more derivatives), so whether this approach can be effective for short PolyBLEPs I honestly don't know. With "long" kernels, it does work reasonably well at least up to cubics though.

In any case, hopefully this will give you some additional insight into the results you have been able to obtain.

Thanks!!