improve PLL fractional divider algorithm

Software/VNA_embedded/Application/Drivers/Si5351C.cpp

@@ -2,6 +2,7 @@
 
 #include <cmath>
 #include <cstdlib>
+#include "algorithm.hpp"
 
 #define LOG_LEVEL	LOG_LEVEL_INFO
 #define LOG_MODULE	"SI5351"
@@ -50,7 +51,7 @@ bool Si5351C::ConfigureCLKIn(uint32_t clkin_freq) {
 	return success;
 }
 
-bool Si5351C::SetPLL(PLL pll, uint32_t frequency, PLLSource src) {
+bool Si5351C::SetPLL(PLL pll, uint32_t frequency, PLLSource src, bool exactFrequency) {
 	if (frequency < 600000000 || frequency > 900000000) {
 		LOG_ERR("Requested PLL frequency out of range (600-900MHz): %lu", frequency);
 		return false;
@@ -67,14 +68,14 @@ bool Si5351C::SetPLL(PLL pll, uint32_t frequency, PLLSource src) {
 		LOG_ERR("Calculated divider out of range (15-90)");
 		return false;
 	}
-	FindOptimalDivider(frequency, srcFreq, c.P1, c.P2, c.P3);
+	FindOptimalDivider(frequency, srcFreq, c.P1, c.P2, c.P3, exactFrequency);
 
 	FreqPLL[(int) pll] = frequency;
 	LOG_DEBUG("Setting PLL %c to %luHz", pll==PLL::A ? 'A' : 'B', frequency);
 	return WritePLLConfig(c, pll);
 }
 
-bool Si5351C::SetCLK(uint8_t clknum, uint32_t frequency, PLL source, DriveStrength strength, uint32_t PLLFreqOverride) {
+bool Si5351C::SetCLK(uint8_t clknum, uint32_t frequency, PLL source, DriveStrength strength, uint32_t PLLFreqOverride, bool exactFrequency) {
 	ClkConfig c;
 	c.DivideBy4 = false;
 	c.IntegerMode = false;
@@ -113,7 +114,7 @@ bool Si5351C::SetCLK(uint8_t clknum, uint32_t frequency, PLL source, DriveStreng
 				return false;
 			}
 		}
-		FindOptimalDivider(pllFreq, frequency * c.RDiv, c.P1, c.P2, c.P3);
+		FindOptimalDivider(pllFreq, frequency * c.RDiv, c.P1, c.P2, c.P3, exactFrequency);
 	}
 	LOG_DEBUG("Setting CLK%d to %luHz", clknum, frequency);
 	return WriteClkConfig(c, clknum);
@@ -346,13 +347,26 @@ bool Si5351C::ResetPLL(PLL pll) {
 }
 
 void Si5351C::FindOptimalDivider(uint32_t f_pll, uint32_t f, uint32_t &P1,
-		uint32_t &P2, uint32_t &P3) {
+		uint32_t &P2, uint32_t &P3, bool exact) {
 	// see https://www.silabs.com/documents/public/application-notes/AN619.pdf (page 3/6)
 	uint32_t a = f_pll / f;
 	int32_t f_rem = f_pll - f * a;
-	// always using the highest modulus divider results in less than 1Hz deviation for all frequencies, that is good enough
+	uint32_t best_c;
-	uint32_t best_c = (1UL << 20) - 1;
+	uint32_t best_b;
-	uint32_t best_b = (uint64_t) f_rem * best_c / f;
+	if(exact) {
+		// spend the effort to find the best divider possible
+		Algorithm::RationalApproximation ratio;
+		ratio.num = f_rem;
+		ratio.denom = f;
+		auto approx = Algorithm::BestRationalApproximation(ratio, (1UL << 20) - 1);
+		best_c = approx.denom;
+		best_b = approx.num;
+	} else {
+		// just use the highest denominator possible, this is good enough if no exact frequency is required
+		best_c = (1UL << 20) - 1;
+		best_b = (uint64_t) f_rem * best_c / f;
+	}
+
 	// convert to Si5351C parameters
 	uint32_t floor = 128 * best_b / best_c;
 	P1 = 128 * a + floor - 512;

Software/VNA_embedded/Application/Drivers/Si5351C.hpp

@@ -31,8 +31,8 @@ public:
 	};
 	bool Init(uint32_t clkin_freq = 0);
 	bool ConfigureCLKIn(uint32_t clkin_freq);
-	bool SetPLL(PLL pll, uint32_t frequency, PLLSource src);
+	bool SetPLL(PLL pll, uint32_t frequency, PLLSource src, bool exactFrequency=true);
-	bool SetCLK(uint8_t clknum, uint32_t frequency, PLL source, DriveStrength strength = DriveStrength::mA2, uint32_t PLLFreqOverride = 0);
+	bool SetCLK(uint8_t clknum, uint32_t frequency, PLL source, DriveStrength strength = DriveStrength::mA2, uint32_t PLLFreqOverride = 0, bool exactFrequency=true);
 	bool SetBypass(uint8_t clknum, PLLSource source, DriveStrength strength = DriveStrength::mA2);
 	bool SetCLKtoXTAL(uint8_t clknum);
 	bool SetCLKToCLKIN(uint8_t clknum);
@@ -48,7 +48,7 @@ public:
 	bool WriteRawCLKConfig(uint8_t clknum, const uint8_t *config);
 	bool ReadRawCLKConfig(uint8_t clknum, uint8_t *config);
 private:
-	void FindOptimalDivider(uint32_t f_pll, uint32_t f, uint32_t &P1, uint32_t &P2, uint32_t &P3);
+	void FindOptimalDivider(uint32_t f_pll, uint32_t f, uint32_t &P1, uint32_t &P2, uint32_t &P3, bool exact);
 	enum class Reg : uint8_t {
 		DeviceStatus = 0,
 		InterruptStatusSticky = 1,

Software/VNA_embedded/Application/Drivers/algorithm.cpp

@@ -2,6 +2,7 @@
 
 #include "stm.hpp"
 #include <cmath>
+#include <algorithm>
 
 Algorithm::RationalApproximation Algorithm::BestRationalApproximation(float ratio, uint32_t max_denom) {
 	RationalApproximation result;
@@ -42,3 +43,47 @@ Algorithm::RationalApproximation Algorithm::BestRationalApproximation(float rati
 	}
 	return result;
 }
+
+uint32_t gcd(uint32_t u, uint32_t v) {
+	if(u==0) {
+		return v;
+	} else if(v==0) {
+		return u;
+	}
+
+	uint8_t i = __builtin_ctz(u);
+	u >>= i;
+	uint8_t j = __builtin_ctz(v);
+	v >>= j;
+
+	uint8_t k = i < j ? i : j;
+	while(true) {
+		if(u > v) {
+			std::swap(u, v);
+		}
+		v -= u;
+		if(v==0) {
+			return u << k;
+		}
+		v >>= __builtin_ctz(v);
+	}
+}
+
+Algorithm::RationalApproximation Algorithm::BestRationalApproximation(
+		RationalApproximation ratio, uint32_t max_denom) {
+	if(ratio.denom <= max_denom) {
+		// nothing to do, we can just return the ratio as it is
+		return ratio;
+	}
+	// Try to simplify the ratio.
+	// Find greatest common divider
+	uint32_t div = gcd(ratio.num, ratio.denom);
+	ratio.num /= div;
+	ratio.denom /= div;
+	if(ratio.denom <= max_denom) {
+		// small enough now, can use as is
+		return ratio;
+	}
+	// no good, we need to approximate
+	return Algorithm::BestRationalApproximation((float) ratio.num / ratio.denom, max_denom);
+}

Software/VNA_embedded/Application/Drivers/algorithm.hpp

@@ -9,6 +9,7 @@ using RationalApproximation = struct _rationalapproximation {
 	uint32_t denom;
 };
 
+RationalApproximation BestRationalApproximation(RationalApproximation ratio, uint32_t max_denom);
 RationalApproximation BestRationalApproximation(float ratio, uint32_t max_denom);
 
 }

Software/VNA_embedded/Application/Drivers/max2871.cpp

@@ -189,7 +189,10 @@ bool MAX2871::SetFrequency(uint64_t f) {
 	LOG_DEBUG("Looking for best fractional match");
 	float fraction = (float) rem_f / f_PFD;
 
-	auto approx = Algorithm::BestRationalApproximation(fraction, 4095);
+	Algorithm::RationalApproximation ratio;
+	ratio.num = rem_f;
+	ratio.denom = f_PFD;
+	auto approx = Algorithm::BestRationalApproximation(ratio, 4095);
 
 	if (approx.denom == approx.num) {
 		// got an impossible result due to floating point limitations(?)

Software/VNA_embedded/Application/SpectrumAnalyzer.cpp

@@ -94,7 +94,7 @@ static void StartNextSample() {
 					}
 					attenuator = amplitude.attenuator;
 					if(trackingFreq < HW::BandSwitchFrequency) {
-						Si5351.SetCLK(SiChannel::LowbandSource, trackingFreq, Si5351C::PLL::A, amplitude.lowBandPower);
+						Si5351.SetCLK(SiChannel::LowbandSource, trackingFreq, Si5351C::PLL::A, amplitude.lowBandPower, false);
 						FPGA::Disable(FPGA::Periphery::SourceChip);
 						FPGA::Disable(FPGA::Periphery::SourceRF);
 						trackingLowband = true;
@@ -198,8 +198,8 @@ static void StartNextSample() {
 	// only adjust LO2 PLL if necessary (if the deviation is significantly less than the RBW it does not matter)
 	if((uint32_t) abs(LO2freq - lastLO2) > actualRBW / 100) {
 //		Si5351.SetPLL(Si5351C::PLL::B, LO2freq*HW::LO2Multiplier, HW::Ref::getSource());
-		Si5351.SetCLK(SiChannel::Port1LO2, LO2freq, Si5351C::PLL::B, Si5351C::DriveStrength::mA2);
+		Si5351.SetCLK(SiChannel::Port1LO2, LO2freq, Si5351C::PLL::B, Si5351C::DriveStrength::mA2, false);
-		Si5351.SetCLK(SiChannel::Port2LO2, LO2freq, Si5351C::PLL::B, Si5351C::DriveStrength::mA2);
+		Si5351.SetCLK(SiChannel::Port2LO2, LO2freq, Si5351C::PLL::B, Si5351C::DriveStrength::mA2, false);
 		lastLO2 = LO2freq;
 	}
 	if (s.UseDFT) {

Software/VNA_embedded/Application/VNA.cpp

@@ -210,7 +210,7 @@ bool VNA::Setup(Protocol::SweepSettings s) {
 	Si5351.Enable(SiChannel::Port1LO2);
 	Si5351.Enable(SiChannel::Port2LO2);
 	Si5351.Enable(SiChannel::RefLO2);
-	Si5351.SetPLL(Si5351C::PLL::B, last_LO2*HW::LO2Multiplier, HW::Ref::getSource());
+	Si5351.SetPLL(Si5351C::PLL::B, last_LO2*HW::LO2Multiplier, HW::Ref::getSource(), false);
 	Si5351.ResetPLL(Si5351C::PLL::B);
 	Si5351.WaitForLock(Si5351C::PLL::B, 10);
 
@@ -458,7 +458,7 @@ void VNA::SweepHalted() {
 			}
 
 			// need the Si5351 as Source
-			bool freqSuccess = Si5351.SetCLK(SiChannel::LowbandSource, frequency, Si5351C::PLL::A, driveStrength);
+			bool freqSuccess = Si5351.SetCLK(SiChannel::LowbandSource, frequency, Si5351C::PLL::A, driveStrength, false);
 			static bool lowbandDisabled = false;
 			if (pointCnt == 0) {
 				// First point in sweep, switch to correct source
@@ -499,7 +499,7 @@ void VNA::SweepHalted() {
 			Si5351.SetCLK(SiChannel::Source, PLLRefFreqs[sourceRefIndex], Si5351C::PLL::A, Si5351C::DriveStrength::mA8);
 			Si5351.SetCLK(SiChannel::LO1, PLLRefFreqs[LO1RefIndex], Si5351C::PLL::A, Si5351C::DriveStrength::mA8);
 			last_LO2 = HW::getIF1() - HW::getIF2();
-			Si5351.SetPLL(Si5351C::PLL::B, last_LO2*HW::LO2Multiplier, HW::Ref::getSource());
+			Si5351.SetPLL(Si5351C::PLL::B, last_LO2*HW::LO2Multiplier, HW::Ref::getSource(), false);
 			Si5351.ResetPLL(Si5351C::PLL::B);
 			Si5351.WaitForLock(Si5351C::PLL::B, 10);
 			HAL_Delay(2);
@@ -514,7 +514,7 @@ void VNA::SweepHalted() {
 				Si5351.SetCLK(SiChannel::LO1, PLLRefFreqs[LO1RefIndex], Si5351C::PLL::A, Si5351C::DriveStrength::mA8);
 			}
 			if(needs2LOshift(frequency, last_LO2, actualBandwidth, &last_LO2)) {
-				Si5351.SetPLL(Si5351C::PLL::B, last_LO2*HW::LO2Multiplier, HW::Ref::getSource());
+				Si5351.SetPLL(Si5351C::PLL::B, last_LO2*HW::LO2Multiplier, HW::Ref::getSource(), false);
 				Si5351.ResetPLL(Si5351C::PLL::B);
 				Si5351.WaitForLock(Si5351C::PLL::B, 10);
 				// PLL reset causes the 2.LO to turn off briefly and then ramp on back, needs delay before next point