McAulay Quatieri, what am I doing wrong???

McAulay Quatieri, what am I doing wrong???

Post by FA » Sun, 20 Oct 2002 02:21:37

Hello All,
 I'm trying to implement the McAulay & Quatieri algorithm (in C++) but I'm
having a bit of a problem with the cubic polynomial phase interpolation bit
(implementing the phaseless model - integrating the frequency works fine).

The following is a code snippet. Would anyone care to shed some light on
what I'm doing wrong.

Trajectory is the class that maintains state in the synthesiser and is
passed the next magnitude/frequency/phase sinusoid (after peak matching
which isn't posted here).

For info playing a single note with several harmonics sounds like the
correct pitch (or close), but there is an unnatural "bell like" quality a
bit like a ring modulator effect. I'm sure it's something silly, but I can't
see what I've done differently to what is stated in their paper.

Regards and TIA, Fraser.

void Trajectory::setState(const Sinusoid& rhs)
 // Sets the frame parameters, called once per frame

 // Set magnitude
 previousMagnitude = magnitude;
 magnitude = (&rhs == this) ? 0 : rhs.magnitude; // If &rhs == this we
assume that the trajectory is dying and so we set magnitude to zero

 // Set new linear interpolation slope
 magnitudeSlope = (magnitude -

 // Set frequency (in hertz)
 previousFrequency = frequency;
 frequency = rhs.frequency;

 // Set phase (in radians and possibly wrapped)
 previousPhase = phase;
 phase = rhs.phase;

 // Set new phase interpolation values from McAulay & Quatieri paper
 const float T = FRAME_DURATION;
 const float T2 = T*T;
 const float T3 = T2*T;

 const float previousW = M_PIx2*previousFrequency;
 const float W = M_PIx2*frequency;

 M = rint((1/M_PIx2)*((previousPhase + previousW*T - phase) + (T/2)*(W -
 alpha = (3/T2)*(phase - previousPhase - previousW*T + M_PIx2*M) - (1/T)
*(W - previousW);
 beta = -(2/T3)*(phase - previousPhase - previousW*T + M_PIx2*M) +
(1/T2)*(W - previousW);

 // Reset sampleCount
 sampleCount = 0;

 // Increment trajectory length


inline float Trajectory::A(int n)
 // Linearly interpolate magnitudes
 return (previousMagnitude + (magnitudeSlope * n));


float Trajectory::tick()
 // Called for each sample

 const float t = (float)sampleCount/(float)(SAMPLE_FREQUENCY);
 const float t2 = t*t;
 const float t3 = t2*t;

 const float previousW = M_PIx2*previousFrequency;

 const float ph = previousPhase + previousW*t + alpha*t2 + beta*t3;

 const float sample = A(sampleCount) * cos(ph);


 return sample;