952860 - DC offset on the square and sawtooth oscillators

Status: RESOLVED FIXED

(Core :: Web Audio, defect)

Description

[Paul Adenot (:padenot)]

STR:
- Go here: http://soledadpenades.com/files/t/cascadiajs-audio-tags/examples/04_chain/
- Set the waveform to square or sawtooth
- Play with the frequency slider

Expected:
- The signal is centered on 0

Actual:
- The signal has a DC offset

The easiest fix is to use the DC blocker we use in the triangle.

Comment 1

[Paul Adenot (:padenot)]

Attachment: Use a DC blocker on the square and triangle oscillators. r=

This changes the shape of the square a bit, but that's way better than having a
massive DC offset.

Comment 2

[Karl Tomlinson (:karlt)]

Comment on attachment 8357168 [details] [diff] [review]
Use a DC blocker on the square and triangle oscillators. r=

>Bug 952860 - Use a DC blocker on the square and triangle oscillators. r=

s/triangle/sawtooth/

I think it would make sense to also DC block mSaw, and mSquare so that they
don't drift.  This would be equivalent to leaky integration, and leaky integration is probably simpler than adding the DC blocker.  Does that make sense?

I suspect there's something else fundamental going on here.  The square wave
at 220Hz for example is relatively stable (when the frequency is not changed)
for about 30 s and then it jumps suddenly.  Perhaps inaccurate integration,
perhaps inaccurate sinc function - I haven't worked it out.  I experimented
with integration using 2 BPBLIT samples per audio sample because we are
essentially integrating one harmonic at the Nyquist - it may have settle
things a bit, but still jumped sometimes.

But we should have some leaky integration regardless.
I guess 0.999 would be time constant of about 20-25 ms at 44.1 kHz.
i.e. a high pass filter at about 40-50 Hz.

Comment 3

[Paul Adenot (:padenot)]

(In reply to Karl Tomlinson (:karlt) from comment #2)
> Comment on attachment 8357168 [details] [diff] [review]
> Use a DC blocker on the square and triangle oscillators. r=
> 
> >Bug 952860 - Use a DC blocker on the square and triangle oscillators. r=
> 
> s/triangle/sawtooth/
> 
> I think it would make sense to also DC block mSaw, and mSquare so that they
> don't drift.  This would be equivalent to leaky integration, and leaky
> integration is probably simpler than adding the DC blocker.  Does that make
> sense?
> 
> I suspect there's something else fundamental going on here.  The square wave
> at 220Hz for example is relatively stable (when the frequency is not changed)
> for about 30 s and then it jumps suddenly.  Perhaps inaccurate integration,
> perhaps inaccurate sinc function - I haven't worked it out.  I experimented
> with integration using 2 BPBLIT samples per audio sample because we are
> essentially integrating one harmonic at the Nyquist - it may have settle
> things a bit, but still jumped sometimes.

The paper from which I implemented this clearly stated that leaky integrator / dc blockers would be needed for this technique to work.

> But we should have some leaky integration regardless.
> I guess 0.999 would be time constant of about 20-25 ms at 44.1 kHz.
> i.e. a high pass filter at about 40-50 Hz.

I now remember why I chose 0.995: its cut frequency is lower than 20Hz, which is kind of the limit under which it does not really matter anymore (considering audio equipment and ears).

fc = -(ln(x) / 2PI) * nyquist

where fc is the cut frequency in Hz (-3db with a one pole filter)
      x is the pole value
      nyquist is half the sampling frequency

-ln((0.995) / 2PI) * (44100 / 2)= 17.590846331216195

Anyways, leaky integrating should work, but I switch to merely using `sLeak` in this patch, and we still get a bad DC offset. I think the problem we have is that we don't do anything when changing frequency, whereas I think we should be having new integration initial condition for the new frequency / phase combination. I'm also pretty sure my maths for the initial integration condition at t = 0 were wrong.

In the meantime, we can probably take a patch like this. 0.995 should be alright, in fact. Hopefully we will have mathematical oscillator soon so using the oscillators as LFO into an AudioParam will be okay (because using a DC blocker with very low freq oscillator will not be great).

Comment 4

[Karl Tomlinson (:karlt)]

Comment on attachment 8357168 [details] [diff] [review]
Use a DC blocker on the square and triangle oscillators. r=

(In reply to Paul Adenot (:padenot) from comment #3)
> The paper from which I implemented this clearly stated that leaky integrator
> / dc blockers would be needed for this technique to work.

Yes, I understood leaky integrators were necessary for floating point rounding
and at frequency changes.  Perhaps they are also enough smooth out the
integration discretization error.  The paper suggests varying the leak
coefficient when the frequency changes.  Perhaps that's an option for future
improvement.

> > But we should have some leaky integration regardless.
> > I guess 0.999 would be time constant of about 20-25 ms at 44.1 kHz.
> > i.e. a high pass filter at about 40-50 Hz.
> 
> I now remember why I chose 0.995: its cut frequency is lower than 20Hz,
> which is kind of the limit under which it does not really matter anymore
> (considering audio equipment and ears).

I guess my time constant to filter cut-off conversion must be wrong.  0.999
should give a lower cut-off and less distortion, so keep that if it works well
enough.  The choice now is a bit arbitrary and we can tune later when other
issues are addressed.

> Anyways, leaky integrating should work, but I switch to merely using `sLeak`
> in this patch, and we still get a bad DC offset.

Are you are suggesting that leaky integration is different from using
DCBlocker, or just that the parameter is different?  I would have expected
0.995 to kill offsets more rapidly, so I don't understand.

-      mSquare += BipolarBLIT();
+      mSquare = 0.999 * mSquare + BipolarBLIT();

-      mSaw += UnipolarBLIT() - dcoffset;
+      mSaw = 0.999 * mSaw + UnipolarBLIT() - dcoffset;

should be equivalent to this patch, and seems to be in my testing.  You can go
ahead and assume my r+ if you are happy to land a patch like that, with a
coefficient of sLeaky or closer to 1, whatever you think best.  (I won't be
around next week.)

> I think the problem we have
> is that we don't do anything when changing frequency, whereas I think we
> should be having new integration initial condition for the new frequency /
> phase combination. I'm also pretty sure my maths for the initial integration
> condition at t = 0 were wrong.

One problem is that the integration sample of the BLIT() should be in the
middle of a phase increment at the current frequency.  Currently the BLIT is
performed at t, to integrate from t - 1/(2*samplerate) to
t + 1/(2*samplerate), assuming the frequency is constant over that period.
The phase increment is applied to increment the phase from t to
t + 1/samplerate assuming the frequency is constant over that period.
When the frequency changes, one of these assumptions is incorrect.
I think either the calculation of the BLIT or the phase increment should be
offset by 1/(2*samplerate).

This also helps consider the initial conditions, which should be for the point 
1/(2*samplerate) before the BLIT sample.

The triangle is more complex because the double integration means that DC
offsets add up before they are smoothed by the leaking.  The mSquare
integration may need to be leakier for the triangle than for the plain square.

I also expect some jumps on frequency changes due to harmonic count changes
changing the positions of the bumps.  Resetting initial conditions may halve
this error if done well but it is tricky to set the initial conditions well
near the step in the function.  Perhaps blending in the harmonic count changes
may be an option.

> In the meantime, we can probably take a patch like this.

Yes.  I'm marking this version r- because the leaky integrator is much
more concise, and easier to see the equation, than using DCBlocker.

Comment 5

[Karl Tomlinson (:karlt)]

Attachment: blit integration experiments

I'm just putting up this patch to record what I experimented with.
These things at least we should do at some point:

* Leaky integrate mSquare for the triangle wave.

* Remove C6 because mSquare has no DC offset.

* Move the mSignalPeriod adjustment on mTriangle into the integrand so it is
  used while it is based on the current frequency.

These were other experiments that I'm not sure about or not sure I have right
yet:

* Integrate UniBLIT at half sample interval offset.

* Get another harmonic for BPBLIT so that frequencies > nyquist/2
  generate output.

* Integrate at BPBLIT at twice the sample rate.

Comment 6

[Paul Adenot (:padenot)]

(In reply to :Karl Tomlinson (back Jan 28) from comment #4)
> Comment on attachment 8357168 [details] [diff] [review]
> Use a DC blocker on the square and triangle oscillators. r=
> 
> (In reply to Paul Adenot (:padenot) from comment #3)
> > The paper from which I implemented this clearly stated that leaky integrator
> > / dc blockers would be needed for this technique to work.
> 
> Yes, I understood leaky integrators were necessary for floating point
> rounding
> and at frequency changes.  Perhaps they are also enough smooth out the
> integration discretization error.  The paper suggests varying the leak
> coefficient when the frequency changes.  Perhaps that's an option for future
> improvement.
> 
> > > But we should have some leaky integration regardless.
> > > I guess 0.999 would be time constant of about 20-25 ms at 44.1 kHz.
> > > i.e. a high pass filter at about 40-50 Hz.
> > 
> > I now remember why I chose 0.995: its cut frequency is lower than 20Hz,
> > which is kind of the limit under which it does not really matter anymore
> > (considering audio equipment and ears).
> 
> I guess my time constant to filter cut-off conversion must be wrong.  0.999
> should give a lower cut-off and less distortion, so keep that if it works
> well
> enough.  The choice now is a bit arbitrary and we can tune later when other
> issues are addressed.
> 
> > Anyways, leaky integrating should work, but I switch to merely using `sLeak`
> > in this patch, and we still get a bad DC offset.
> 
> Are you are suggesting that leaky integration is different from using
> DCBlocker, or just that the parameter is different?  I would have expected
> 0.995 to kill offsets more rapidly, so I don't understand.

Looking back at my patch where I thought I implemented the leaky integration variant, I now see I wrote something like:

> mSaw += (UnipolarBLIT() - dcoffset) * sLeak;

Which is not a leaky integrator and does not work. Sorry about that.

> 
> -      mSquare += BipolarBLIT();
> +      mSquare = 0.999 * mSquare + BipolarBLIT();
> 
> -      mSaw += UnipolarBLIT() - dcoffset;
> +      mSaw = 0.999 * mSaw + UnipolarBLIT() - dcoffset;
> 
> should be equivalent to this patch, and seems to be in my testing.  You can
> go
> ahead and assume my r+ if you are happy to land a patch like that, with a
> coefficient of sLeaky or closer to 1, whatever you think best.  (I won't be
> around next week.)

Yeah, I might tweak the coeff, we'll see.

> 
> > I think the problem we have
> > is that we don't do anything when changing frequency, whereas I think we
> > should be having new integration initial condition for the new frequency /
> > phase combination. I'm also pretty sure my maths for the initial integration
> > condition at t = 0 were wrong.
> 
> One problem is that the integration sample of the BLIT() should be in the
> middle of a phase increment at the current frequency.  Currently the BLIT is
> performed at t, to integrate from t - 1/(2*samplerate) to
> t + 1/(2*samplerate), assuming the frequency is constant over that period.
> The phase increment is applied to increment the phase from t to
> t + 1/samplerate assuming the frequency is constant over that period.
> When the frequency changes, one of these assumptions is incorrect.
> I think either the calculation of the BLIT or the phase increment should be
> offset by 1/(2*samplerate).
> 
> This also helps consider the initial conditions, which should be for the
> point 
> 1/(2*samplerate) before the BLIT sample.
> 
> The triangle is more complex because the double integration means that DC
> offsets add up before they are smoothed by the leaking.  The mSquare
> integration may need to be leakier for the triangle than for the plain
> square.
> 
> I also expect some jumps on frequency changes due to harmonic count changes
> changing the positions of the bumps.  Resetting initial conditions may halve
> this error if done well but it is tricky to set the initial conditions well
> near the step in the function.  Perhaps blending in the harmonic count
> changes
> may be an option.
> 
> > In the meantime, we can probably take a patch like this.
> 
> Yes.  I'm marking this version r- because the leaky integrator is much
> more concise, and easier to see the equation, than using DCBlocker.

I'll land the variant with the leaky integrator, and open a new bug quoting this last paragraph and comment 5.

Comment 7

[Paul Adenot (:padenot)]

Attachment: Leaky integrate sawtooth and square waves to prevent bad DC offset when changing frequency.

0.999 _looks_ better, and I can't really hear any difference, so I'll land that.

Comment 8

[Paul Adenot (:padenot)]

https://hg.mozilla.org/integration/mozilla-inbound/rev/58355a630c4c

Comment 9

[Ryan VanderMeulen [:RyanVM]]

https://hg.mozilla.org/mozilla-central/rev/58355a630c4c