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

inline float Trajectory::A(int n)Quote:}

{

// Linearly interpolate magnitudes

return (previousMagnitude + (magnitudeSlope * n));

float Trajectory::tick()Quote:}

{

// 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:}