McAulay Quatieri, what am I doing wrong???

McAulay Quatieri, what am I doing wrong???

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 -
previousMagnitude)/(SAMPLE_FREQUENCY*FRAME_DURATION);

// 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 -
previousW)));
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
length++;

Quote:}

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

Quote:}

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);

sampleCount++;

return sample;

Quote:}

Hi,
Here, I never read about people struggling with Precise' MQX operating
system for the DSP. And nobody replied to my MQX related questions.

Is this the wrong newsgroup for such stuff?