[CM] Snd/CLM: experimental version of src and src-channel

[Tito Latini]

Hi,

the attached short patch changes mus_src() and mus_src_to_buffer()
in snd/clm.c. Probably it is compatible with clm-5/clm.c (no tested).

It seems good with constant rates or envelopes. The aliasing is
remarkably reduced and the amplitude of the resampled sound is very good.

The conversion is performed through the convolution between the input,
upsampled by rate, and the polyphase filters obtained from the windowed
sinc table. The delay between the polyphase filters is not one sample,
but it depends on the rate. For example, if rate is 1.9, the phase shift
is pi/5:

    2 pi (1 - frac(1.9)) = 2 pi 0.1 = pi/5

therefore the index of the first coefficient for the next polyphase
filter is

    SINC-TABLE-DENSITY * 0.1 = 2000 * 0.1 = 200

In practice, sinc_loc (the table index) and sinc_incr are reduced to

    sinc_loc = SRC_SINC_DENSITY * (1 - srpx) + 4;
    sinc_incr = SRC_SINC_DENSITY;

The increment sinc_incr is the consequence of the polyphase filtering:
a coefficient after (SRC_SINC_DENSITY - 1) zeroes.

I'm getting also better results with a Kaiser window (not a surprise).
It is possible to apply that window if a window parameter, for example
window-beta, is set for make-src (not in patch), otherwise the default
remains the Hanning window. Kaiser window is not particularly
interesting with the current implementation (but we increase the table
length), because a compressed interpoled sinc for non-integer rates >1
doesn't start where the window begins, so the sinc is often truncated).
On the contrary, the index of a polyphase filter is always near the
left side (less than SRC_SINC_DENSITY in the upsampled FIR).

What do you think?
-------------- next part --------------
--- snd/clm.c~	2020-09-10 16:06:48.774551471 +0200
+++ snd/clm.c	2020-09-10 16:05:39.163555210 +0200
@@ -13298,8 +13298,8 @@
 mus_float_t mus_src(mus_any *srptr, mus_float_t sr_change, mus_float_t (*input)(void *arg, int direction))
 {
   sr *srp = (sr *)srptr;
-  mus_float_t sum, zf, srx, factor;
-  int lim, loc, xi;
+  mus_float_t sum, srx;
+  int lim, loc;
   bool int_ok;
   mus_float_t *data, *sinc_table;
 
@@ -13351,9 +13351,10 @@
   if (srx < 0.0) srx = -srx;
   if (srx > 1.0) 
     {
-      factor = 1.0 / srx;
+      mus_float_t zf;
+      int xi;
       /* this is not exact since we're sampling the sinc and so on, but it's close over a wide range */
-      zf = factor * (mus_float_t)SRC_SINC_DENSITY; 
+      zf = SRC_SINC_DENSITY / srx;
       xi = (int)(zf + 0.5);
 
       /* (let ((e (make-env '(0 1 1 1.1) :length 11))) (src-channel e))
@@ -13363,20 +13364,16 @@
     }
   else 
     {
-      factor = 1.0;
-      zf = (mus_float_t)SRC_SINC_DENSITY;
-      xi = SRC_SINC_DENSITY;
       int_ok = true;
     }
 
   sum = 0.0;
   if (int_ok)
     {
-      int sinc_loc, sinc_incr, last, last10, xs;
+      int sinc_loc, sinc_incr, last, last10;
       
-      xs = (int)(zf * (srp->width_1 - srp->x));
-      sinc_loc = xs + srp->sinc4;
-      sinc_incr = xi;
+      sinc_loc = SRC_SINC_DENSITY - (int)(SRC_SINC_DENSITY * srp->x) + 4;
+      sinc_incr = SRC_SINC_DENSITY;
       last = loc + lim;
       last10 = last - 10;
 
@@ -13398,12 +13395,11 @@
     }
   else
     {
-      mus_float_t sinc_loc, sinc_incr, x;
+      mus_float_t sinc_loc, sinc_incr;
       int last, last10;
 
-      x = zf * (srp->width_1 - srp->x);
-      sinc_loc = x + srp->sinc4;
-      sinc_incr = zf;
+      sinc_loc = SRC_SINC_DENSITY * (1 - srp->x) + 4;
+      sinc_incr = SRC_SINC_DENSITY;
       last = loc + lim;
 
       last10 = last - 10;
@@ -13426,7 +13422,7 @@
     }
 
   srp->x += srx;
-  return(sum * factor);
+  return(sum);
 }
 
 
@@ -13435,8 +13431,8 @@
   /* sr_change = 0.0
    */
   sr *srp = (sr *)srptr;
-  mus_float_t x, zf, srx, factor, sincx, srpx;
-  int lim, i, xi, xs, dir = 1;
+  mus_float_t srx, sincx, srpx;
+  int lim, i, dir = 1;
   bool int_ok;
   mus_long_t k;
   mus_float_t *data, *sinc_table;
@@ -13454,17 +13450,15 @@
     }
   if (srx > 1.0) 
     {
-      factor = 1.0 / srx;
+      mus_float_t zf;
+      int xi;
       /* this is not exact since we're sampling the sinc and so on, but it's close over a wide range */
-      zf = factor * sincx;
-      xi = (int)zf;
+      zf = sincx / srx;
+      xi = (int)(zf + 0.5);
       if (fabs((xi - zf) * lim) > 2.0) int_ok = false; else int_ok = true;
     }
   else 
     {
-      factor = 1.0;
-      zf = sincx;
-      xi = SRC_SINC_DENSITY;
       int_ok = true;
     }
 
@@ -13502,10 +13496,9 @@
       if (int_ok)
 	{
 	  int sinc_loc, sinc_incr, last, last10;
-	  
-	  xs = (int)(zf * (srp->width_1 - srpx));
-	  sinc_loc = xs + srp->sinc4;
-	  sinc_incr = xi;
+
+	  sinc_loc = SRC_SINC_DENSITY - (int)(SRC_SINC_DENSITY * srpx) + 4;
+	  sinc_incr = SRC_SINC_DENSITY;
 	  last = loc + lim;
 	  last10 = last - 10;
 	  
@@ -13529,10 +13522,9 @@
 	{
 	  mus_float_t sinc_loc, sinc_incr;
 	  int last, last10;
-	  
-	  x = zf * (srp->width_1 - srpx);
-	  sinc_loc = x + srp->sinc4;
-	  sinc_incr = zf;
+
+	  sinc_loc = SRC_SINC_DENSITY * (1 - srpx) + 4;
+	  sinc_incr = SRC_SINC_DENSITY;
 	  last = loc + lim;
 	  
 	  last10 = last - 10;
@@ -13554,7 +13546,7 @@
 	    sum += data[loc] * sinc_table[(int)sinc_loc];
 	}
       srpx += srx;
-      out_data[k] = sum * factor;
+      out_data[k] = sum;
     }
   srp->x = srpx;
 }