Synchronize a software waveform generator

Question

I'm having some trouble with a software waveform generator I'm attempting to put together. Especially with synchronizing the waveform.

Goal

Generate a mains synchronized (50 Hz) sinusoidal waveform to output on the DAC of a STM32F103.

Current Method

I developed a NCO (numerical controller oscillator) based on the accumulator method from this wiki. This works, but my first attempt did not have precise enough calculations for the accumulator increment to get a reasonable close to the original frequency, so I had to move to q7.24 for this.
My samplerate is 10 kHz, with an interrupt every 2 ms.

Problem

I'm measuring the mains frequency with a 4 MHz free-running timer, and calculate the frequency to provide a q15.16 variable. This provides the frequency with reasonable precision. It is satisfactory for all other products that require frequency measurement.
It is however, not accurate enough for this purpose. There will be error due to the deviation of the 12 MHz crystal.

This is causing a problem, since the signal will need to synchronize to this 50 Hz reference. Which it will need some feedback for.

I developed phase measurement using the square wave subtraction method. This can calculate the phase within the limits of the samplerate.
With this result I feed a PI control loop that adds or subtracts from the accumulator increment in the NCO.

This keeps it in phase, but it still drifts around significantly on changes of mains. (and it changes a lot)

Is there a better way to synchronize a software numerical controlled oscillator? Or can I improve my current method?

Code

// Calculate increment for ac waveform generator accumulator
uint64_t f = (FREQ()->freq) << 8;  // Convert frequency q15.16 to q7.24       
uint64_t t = (f * 1678);           // Multiply by sampleperiod, f * (1/fs) = f * 0.0001
inc = t / (1<<24);                 // q24*q24=q48, need to shift back to q24
// 50 hz approx 83886 q24

// Generate AC Waveform
// Add to (f * (1/fs)) to accumulator
accumlator24 += inc;
accumlator24 += phase_correction;
// Simulate 24 bit unsigned integer
if(accumlator24 & 0xFF000000u){
  accumlator24 = accumlator24 & 0x00FFFFFFu;
}
// Convert accumulator from q24 to q16 for fixed point sine
uint32_t accumulator16 = accumlator24 / 256;
// Multiply q16 accumulator "0.0 to 1.0" with 2*PI
x = fix16_mul(accumulator, fix16_2pi);
x = fix16_cos(x);
// Fill phase detect buffer (bitbanded)
if( x > 0 ){
  WFG_phase_bb[WFG_phase_dac_i++] = 1;
}else{
  WFG_phase_bb[WFG_phase_dac_i++] = 0;
}
// Multiply process value
x = fix16_smul(x, amplitude);
// Move to DAC range (0.0-1.0)
x = fix16_sadd(x, F16C(1,0));
x = fix16_smul(x, F16C(0,5));    
x = fix16_clamp(x, 0, 0xFFFF);
// Copy lsb of fix16 (0.0-1.0) to dac buffer
dac1_data[i] = x;

---

Solution

In my implementation above the phase measuring of the generated vs reference waveform was lacking.
Because the resolution of this is equal to the samplerate due to the reference being a binary signal. I have improved the phase measurement from the binary comparison to the multiplication method. This allows for the mentioned topologys by Dan Boschen to be used.

I've added input normalization. This is actually remarkable easy to do if the signal is clean. Multiply by (1/peak).

if( input > peak ){
    peak = input;
}else{
    peak *= 0.995; // decay the peak, simulates a rectifier
}
factor = 1 / peak;
normalized = input * factor.

This normalized signal can be multiplied with the generated signal. This gives the phase as offset the multiplied signal. A low pass filter is used to get a stable value.

The phase detection with the multiplication method is proven to be accurate and provides a stable phase control loop.

Answer 1

It is not clear to me how you are determining your phase error so there may be issues there, and it also looks like you are adding the phase correction to the accumulator directly. Consider adding the corrections to the frequency control word instead which would be analogous to a traditional analog PLL such as that shown in the graphic below (with the NCO being the equivalent to a VCO); this would allow the NCO to be an additional integrator in the loop which is necessary to implement a desired second order type 2 PI loop with your use of a PI Loop filter. If you feed the corrections to the phase accumulator directly, you are treating the NCO as a phase control device instead of a frequency control device.

A type 2 implementation will track a frequency offset with zero error, while a type 1 implementation (which is what would result if the control was fed to the phase accumulator directly) would have an offset that is proportional to the inverse of the loop gain.

For the phase detector implementation, consider using a simple multiplier followed by a low pass filter to remove the double frequency output. (This multiplier can be analog or digital-- in your case since you are dealing with an analog signal an analog multiplier may be simplest. If you can be ensured the duty cycle is 50%, an XOR gate is perfect for this purpose. Otherwise a passive mixer is often a good choice). The average value will be proportional the phase error with a quadrature offset (so the Loop would therefore lock in quadrature to the input signal). Care should be taken that the low pass filter simply removes the double frequency component and is not too tight in bandwidth so as to compete with the loop bandwidth (it is preferred to set the loop bandwidth with the loop filter alone). Depending on implementation (see diagrams below) a "Phase/Frequency Detector" may be a preferred choice for the phase detector given that it won't lock to alias frequency offsets. If a multiplier type phase detector is used, the VCO/NCO output must be constrained to a narrow frequency range to avoid such false locks.

Three variants of a PLL approach are shown below. Depending on your need for digital and analog signal types in the rest of your system and programmability of loop dynamics one may be preferred over the other.

For the last implementation shown (which I think would be simplest) there are plenty of integrated chip solutions out there to make this even simpler. For example this chip: https://www.analog.com/en/products/hmc984.html#product-overview incorporates a phase/frequency detector together with a high range divider. The divider is 14 bits but 18 bits are required to get to 240,000, so the approach in this case would be to "pre-scale" the 12 MHz VCXO output with a divide by 16 to create a 750KHz clock, which would then be divided by 15,000 with this chip. There are many other integrated solutions to consider if you went down this path.