From 8d6f3bcb96043a590dd0ef60faf194066195d1cc Mon Sep 17 00:00:00 2001
From: Marton Balint <cus@passwd.hu>
Date: Sun, 10 Nov 2024 01:14:06 +0100
Subject: [PATCH] avfilter/asrc_sine: increase frequency accuracy

Previously the delta phase was fixed point fractional with 2^32 fractions,
which caused inaccuracies in the output frequency, unless the input
frequency*2^32 was divisable by the sample rate.

This patch improves frequency accuracy by tracking subfractions of the delta
phase fractions. For this we are using a denominator which is a multiple of the
sample rate, making sure that integer frequencies are always accurately
represented.

Signed-off-by: Marton Balint <cus@passwd.hu>
---
 libavfilter/asrc_sine.c | 24 +++++++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/libavfilter/asrc_sine.c b/libavfilter/asrc_sine.c
index 854fba6e37..6a591b551b 100644
--- a/libavfilter/asrc_sine.c
+++ b/libavfilter/asrc_sine.c
@@ -33,6 +33,9 @@
 typedef struct SamplingContext {
     uint32_t phi;  ///< current phase of the sine (2pi = 1<<32)
     uint32_t dphi; ///< phase increment between two samples
+    int phi_rem;   ///< current fractional phase in 1/dphi_den subfractions
+    int dphi_rem;
+    int dphi_den;
 } SamplingContext;
 
 typedef struct SineContext {
@@ -148,12 +151,31 @@ enum {
 
 static void sampling_init(SamplingContext *c, double frequency, int sample_rate)
 {
-    c->dphi = ldexp(frequency, 32) / sample_rate + 0.5;
+    AVRational r;
+    int r_den, max_r_den;
+
+    max_r_den   = INT_MAX / sample_rate;
+    frequency   = fmod(frequency, sample_rate);
+    r           = av_d2q(fmod(frequency, 1.0), max_r_den);
+    r_den       = FFMIN(r.den, max_r_den);
+    c->dphi     = ldexp(frequency, 32) / sample_rate;
+    c->dphi_den = r_den * sample_rate;
+    c->dphi_rem = round((ldexp(frequency, 32) / sample_rate - c->dphi) * c->dphi_den);
+    if (c->dphi_rem >= c->dphi_den) {
+        c->dphi++;
+        c->dphi_rem = 0;
+    }
+    c->phi_rem  = (-c->dphi_den - 1) / 2;
 }
 
 static av_always_inline void sampling_advance(SamplingContext *c)
 {
     c->phi += c->dphi;
+    c->phi_rem += c->dphi_rem;
+    if (c->phi_rem >= 0) {
+        c->phi_rem -= c->dphi_den;
+        c->phi++;
+    }
 }
 
 static av_cold int init(AVFilterContext *ctx)