Sensor Fusion Implementation - Drift Correction Problems

I am trying to implement Gyro and Accelerometer sensor fusion, as described in "Motion Sensors Introduction" page 6. The goal is to correct roll and pitch gyro drift with gravity. Infortunately, I end up with a result that either spins continuously or drifts to the wrong axis (and no mater how I change those it doesn't match up).

Step (1): Obtaining position q(t). I have this working correctly.
Step (2): Convert accelerometer to world coordinates. Aw(t) = q(t) * Ab(t) * q(t)^-1. I think I have this as well.
Step (3): Create feedback quaternion, qf(t) = [0 Awy(t) - Awx(t) 0] * gain. This is giving me problems.

First of all, I am now sure wher to put commas. I assume the intent is qf(t) = [0, Awy(t), -Awx(t), 0]. Second, this spins like crazy unless I set the real component to 1 and not 0. Finally, I had to apply gain to factors and not the whole quaternion. In the end, it doesn't drift into the right place, especially after a turn. Sometimes is starts spinning in alternating axes.

Here's the complete code I have:

    // This works; Q gives proper rotation. Here avgAngV is in rad/s.

Quatf qYaw(Vector3f(0,1,0), avgAngV.y * dt);

Quatf qRoll(Vector3f(0,0,1), avgAngV.z * dt);

Quatf qPitch(Vector3f(1,0,0), avgAngV.x * dt);

Quatf q = qYaw * qPitch * qRoll;

Q = q * Q;



A = VectorFromAccelUpdate(update);

// This doesn't work:

float gain = 0.05f;

Quatf Qw   = Q * Quatf(A.x * dt, A.y * dt, A.z * dt, 0) * Q.Inversed();

Quatf qfeedback(Qw.y * gain, -Qw.x * gain, 0, 1);

Q = qfeedback * Q;

Q.Normalize();

I wonder if there is any sample code that demonstrates this? Has anyone tried implementing this approach?

Part of the problems is that I don't truly understand the math involved in the part(3). Multiplying quaternions accumulates rotation. Wouldn't adding a gravity-based acceleration result in continuous spin - why should it stabilize at all?