changeset 396:0430226a0ea2 resample

Start work on direct resampler

--- a/.hgsubstate	Tue Sep 24 17:31:14 2013 +0100
+++ b/.hgsubstate	Wed Sep 25 18:32:55 2013 +0100
@@ -1,1 +1,1 @@
-43584910bba1d16db21b2ce6d8feafe87c9361fa ext
+4b91cd40a51b64283ba6b5550709e89f6ec492db ext

--- a/src/may/signal/window.yeti	Tue Sep 24 17:31:14 2013 +0100
+++ b/src/may/signal/window.yeti	Wed Sep 25 18:32:55 2013 +0100
@@ -116,7 +116,7 @@
             done [0..n-1]));
 
 /** 
-  Kaiser window with sidelobe attenuation of 𝛼 dB and window length n
+  Kaiser window with sidelobe attenuation of alpha dB and window length n
 */
 kaiserForAttenuation alpha n =
    (beta =

--- a/src/may/stream/filter.yeti	Tue Sep 24 17:31:14 2013 +0100
+++ b/src/may/stream/filter.yeti	Wed Sep 25 18:32:55 2013 +0100
@@ -11,6 +11,8 @@
 syn = load may.stream.syntheticstream;
 cplx = load may.complex;
 
+plot = load may.plot;
+
 load may.stream.type;
 load may.vector.type;
 load may.matrix.type;
@@ -426,11 +428,24 @@
 /**
   Produce a sinc window of width equal to zc zero crossings per half,
   with n samples from peak to first zero crossing, multiplied by a
-  Kaiser window with attenuation 𝛼 dB
+  Kaiser window with attenuation alpha dB
  */
-kaiserSincWindow zc n 𝛼 =
+kaiserSincWindow zc n alpha =
    (sw = win.sinc (n*2) (n*zc*2 + 1);
-    kw = win.kaiserForAttenuation 𝛼 (vec.length sw);
+    kw = win.kaiserForAttenuation alpha (vec.length sw);
+    bf.multiply sw kw);
+
+/**
+ * Produce a sinc window with n samples from peak to first zero
+ * crossing, multiplied by a Kaiser window with attenuation alpha dB
+ * and transition bandwidth Hz at the given sampleRate. The filter
+ * will contain an odd number of samples.
+ */
+kaiserSincFilterFor { n, alpha, bandwidth, sampleRate } =
+   (kw = win.kaiserForBandwidth alpha bandwidth sampleRate;
+    kw = if (vec.length kw) % 2 == 0 then vec.concat [vec.zeros 1, kw] else kw fi;
+    sw = win.sinc (n*2) (vec.length kw);
+    println "kaiserSincFilterFor: filter length is \(vec.length sw)";
     bf.multiply sw kw);
 
 bandpassed f0 f1 attenuation bandwidth s =
@@ -499,10 +514,9 @@
 interpolated factor s = //!!! factor must be an integer [how to enforce this??]
     if factor == 1 then s
     else
-        nzc = 13; //!!! where does this magic come from?
+        nzc = 13;
         attenuation = 80;
         filter = kaiserSincWindow nzc factor attenuation;
-println "int: filter length for nzc \(nzc), attn \(attenuation), factor \(factor) is \(vec.length filter)";
         out = delayedBy (- (nzc * factor))
            (convolvedWith [Framesize 1024]
                (mat.newMatrix (RowMajor ()) (map \filter [1..s.channels]))
@@ -533,7 +547,6 @@
         nzc = 13;
         attenuation = 80;
         filter = kaiserSincWindow nzc factor attenuation;
-println "dec: filter length for nzc \(nzc), attn \(attenuation), factor \(factor) is \(vec.length filter)";
         filtered =
            (convolvedWith [Framesize 1024]
                (mat.newMatrix (RowMajor ()) (map \filter [1..s.channels])) s);
@@ -546,11 +559,122 @@
 gcd a b = (c = a % b; if c == 0 then b else gcd b c fi);
 
 //!!! slow -- can't in this form be used for e.g. 44100 -> 48000 conversion
-resampledTo rate s =
+resampledSlowlyTo rate s =
    (g = gcd rate s.sampleRate;
    println "orig = \(s.sampleRate), target = \(rate), g = \(g), up = \(rate/g), down = \(s.sampleRate/g)";
     decimated (s.sampleRate / g) (interpolated (rate/g) s));
 
+resampledSplendidlyTo targetRate s =
+   (
+    // We need a low-pass filter with a cutoff of the lower of the two
+    // (nyquist frequencies of the) sample rates. This filter needs to
+    // be sampled at a sufficient resolution that we can get values
+    // from it at every sample point in the source rate and in the
+    // target.
+
+    sourceRate = s.sampleRate;
+    higher = max sourceRate targetRate;
+    lower = min sourceRate targetRate;
+    g = gcd higher lower;
+
+    // Example: ratio of 3:4 in sample periods, corresponding to
+    // e.g. resampling from 48kHz to 36kHz. The lower rate has the
+    // longer sample period. We need a filter with n=4 values from
+    // peak to first pole -- this is the longer sample period relative
+    // to the shorter and is obtained by dividing the *higher* rate by
+    // the gcd (because the ratio of sample periods is the reciprocal
+    // of the ratio of sample rates).
+
+    n = higher / g; // would be 4 in the above example
+
+    println "for target rate \(targetRate) and source rate \(sourceRate), gcd = \(g) and n = \(n), would be \(lower / g) if we used min";
+
+    // Our filter is a sinc function with peak-to-pole length n,
+    // multiplied by a Kaiser window of transition bandwidth narrow
+    // enough at the effective sample rate n to exclude... what?
+
+    //!!! ponder
+//    filt = kaiserSincFilterFor { n, alpha = 80, bandwidth = 10, sampleRate = lower };
+    filt = kaiserSincFilterFor { n, alpha = 80, bandwidth = 1, sampleRate = n };
+
+//\() (    plot.plot [Vector filt] );
+
+    // Now we have a filter of (odd) length flen in which the lower
+    // sample rate corresponds to every n'th point and the higher rate
+    // to every m'th where n and m are higher and lower rates divided
+    // by their gcd respectively. So if x coordinates are on the same
+    // scale as our filter resolution, then source sample i is at i *
+    // (targetRate / gcd) and target sample j is at j * (sourceRate /
+    // gcd).
+
+    // To reconstruct a single target sample, we want a buffer (real
+    // or virtual) of flen values formed of source samples spaced at
+    // intervals of (targetRate / gcd), in our example case 3.  This
+    // is initially formed with the first sample at the filter peak.
+    //
+    // 0  0  0  0  a  0  0  b  0
+    //
+    // and of course we have our filter
+    //
+    // f1 f2 f3 f4 f5 f6 f7 f8 f9
+    //
+    // We take the sum of products of non-zero values from this buffer
+    // with corresponding values in the filter
+    //
+    // a * f5 + b * f8
+    //
+    // Then we drop (sourceRate / gcd) values, in our example case 4,
+    // from the start of the buffer and fill until it has flen values
+    // again
+    //
+    // a  0  0  b  0  0  c  0  0
+    //
+    // repeat to reconstruct the next target sample
+    //
+    // a * f1 + b * f4 + c * f7
+    //
+    // and so on.
+
+    flen = vec.length filt;
+    halflen = int (flen/2); // actual filter length is halflen + halflen + 1
+    spacedInput = spaced (targetRate / g) s;
+    var buffer = mat.toRowMajor
+       (mat.concat (Horizontal ())
+           [mat.zeroMatrix { rows = s.channels, columns = halflen + 1 },
+            spacedInput.read halflen]);
+    var pos = 0;
+
+    reconstructOne ch =
+       (var out = 0;
+        var series = mat.getRow ch buffer;
+        for [0..(vec.length series)-1] do i:
+            x = vec.at series i;
+            if x != 0 then
+                out := out + x * (vec.at filt i);
+            fi;
+        done;
+        out);
+
+    readOne () =
+       (result = mat.newColumnVector
+           (vec.fromList (map reconstructOne [0 .. s.channels-1]));
+        m = sourceRate / g;
+        buffer := mat.concat (Horizontal ())
+           [mat.columnSlice buffer m (flen - m),
+            spacedInput.read m];
+        result);
+    
+    s with
+    {
+        get sampleRate () = targetRate,
+        get position () = pos,
+        get available () = Unknown (), //!!! todo
+        read n = mat.toRowMajor
+           (mat.concat (Horizontal ()) (map do _: readOne () done [1..n])),
+    });
+
+resampledTo = resampledSplendidlyTo;
+
 {
     withDuration is number -> stream -> stream,
     delayedBy is number -> stream -> stream,
@@ -564,7 +688,7 @@
     repeated is stream -> stream,
     duplicated is number -> stream -> array<stream>,
     convolvedWith is list<Fast boolean | Framesize number> -> matrix -> stream -> stream,
-    kaiserSincWindow,
+    kaiserSincWindow, kaiserSincFilterFor,
     lowpassed, bandpassed, highpassed,
     spaced, interpolated, decimated, picked,
     resampledTo,