view src/may/stream/test/test_resample.yeti @ 568:e9716e9a44e4

Fixes to tests for frequency preservation in resampling

module may.stream.test.test_resample;

vec = load may.vector;
mat = load may.matrix;
syn = load may.stream.syntheticstream;
manip = load may.stream.manipulate;
channels = load may.stream.channels;
resample = load may.stream.resample;
win = load may.signal.window;
waves = load may.stream.waves;
fft = load may.transform.fft;

{ compare, compareUsing, compareMatrices, assert, time } = load may.test;

eps = 1e-14;

measureFrequency cycles stream =
   (case stream.available of
    Known n:
       (v = channels.mixedDown (stream.read n);
        var firstPeak = -1;
        var lastPeak = -1;
        var nPeaks = 0;
        // count +ve peaks
        len = vec.length v;
        for [int(len/4) .. len - 2] do i:
            // allow some fuzz
            x0 = int (10000 * vec.at v (i-1));
            x1 = int (10000 * vec.at v (i));
            x2 = int (10000 * vec.at v (i+1));
            if nPeaks < (cycles + 1) and
                x1 > 0 and x1 > x0 and x1 >= x2 then
                if firstPeak < 0 then firstPeak := i; fi;
                lastPeak := i;
                nPeaks := nPeaks + 1;
            fi;
        done;
        if nPeaks < 2 then 0
        else
            cycle = (lastPeak - firstPeak) / (nPeaks - 1);
            stream.sampleRate / cycle;
        fi);
   _: failWith "Expected stream duration to be known";
   esac);

testFrequencyIntegrityWith sort freq sourceRate forward backward =
   (// Test that downsample and then upsample again (or vice versa,
    // depending on forward and backward manipulators) produces a
    // signal at the same frequency as the original (i.e. the overall
    // speed is retained exactly, regardless of the SNR on the way).
    nCycles = 200;
    duration = int (nCycles * sourceRate / freq);
    originals = manip.duplicated 3
       (manip.withDuration duration (syn.sinusoid sourceRate freq));
    there = manip.duplicated 3 (forward originals[0]);
    back = manip.duplicated 3 (backward there[0]);
    origFreq = measureFrequency (nCycles/2) originals[1];
    backFreq = measureFrequency (nCycles/2) back[1];
    if not (compare backFreq origFreq) then
        eprintln "** note: rate conversion \(sourceRate) -> \(there[1].sampleRate) -> \(sourceRate)";
        examples = array (map do s:
            channels.mixedDown (s[2].read 500)
        done [ originals, there, back ]);
        false
    else
        true
    fi);

testFrequencyIntegrity freq sourceRate targetRate =
    if sourceRate / targetRate == int (sourceRate / targetRate) then
        testFrequencyIntegrityWith "decimated" freq sourceRate
           (resample.decimated (sourceRate / targetRate))
           (resample.interpolated (sourceRate / targetRate))
    elif targetRate / sourceRate == int (targetRate / sourceRate) then
        testFrequencyIntegrityWith "interpolated" freq sourceRate
           (resample.interpolated (targetRate / sourceRate))
           (resample.decimated (targetRate / sourceRate))
    else
        true
    fi and
        testFrequencyIntegrityWith
           (if sourceRate > targetRate then "downsampled" else "upsampled" fi)
            freq sourceRate
           (resample.resampledTo targetRate)
           (resample.resampledTo sourceRate);

[

// Test for duration of decimated stream (does not test contents, that
// happens in the filters tests below).
// Resampled streams only have accurate duration when their inputs
// have known duration (otherwise they may be rounded up to the
// processing block size).

"dec-duration": \(
    original = manip.withDuration 8 (syn.precalculatedMono 4 (vec.fromList [1,2,3,4,5,6,7,8]));
    str = resample.decimated 2 original;
    compare str.position 0 and
        compare str.channels 1 and
        compare str.sampleRate 2 and
        compare str.available (Known 4) and
        compare str.finished? false and
       (r = str.read 3;
        compare (mat.size r) { rows = 1, columns = 3 }) and
        compare str.position 3 and
        compare str.available (Known 1) and
        compare str.finished? false and
       (r = str.read 3;
        compare (mat.size r) { rows = 1, columns = 1 }) and
        compare str.position 4 and
        compare str.available (Known 0) and
        compare str.finished? true and
        ( str.close (); true )
),

// Test for duration of interpolated stream (does not test contents,
// that happens in the filters tests below)
// Resampled streams only have accurate duration when their inputs
// have known duration (otherwise they may be rounded up to the
// processing block size).

"int-duration": \(
    original = manip.withDuration 8 (syn.precalculatedMono 4 (vec.fromList [1,2,3,4,5,6,7,8]));
    str = resample.interpolated 2 original;
    compare str.position 0 and
        compare str.channels 1 and
        compare str.sampleRate 8 and
        compare str.available (Known 16) and
        compare str.finished? false and
        compare (mat.size (str.read 12)) { rows = 1, columns = 12 } and
        compare str.position 12 and
        compare str.available (Known 4) and
        compare str.finished? false and
        compare (mat.size (str.read 12)) { rows = 1, columns = 4 } and
        compare str.position 16 and
        compare str.available (Known 0) and
        compare str.finished? true and
        ( str.close (); true )
),

// Test for duration of resampled stream (does not test contents,
// that happens in the filters tests below)
// Resampled streams only have accurate duration when their inputs
// have known duration (otherwise they may be rounded up to the
// processing block size).
"resamp-duration": \(
    original = manip.withDuration 8 (syn.precalculatedMono 4 (vec.fromList [1,2,3,4,5,6,7,8]));
    str = resample.resampledTo 6 original;
    compare str.position 0 and
        compare str.channels 1 and
        compare str.sampleRate 6 and
        compare str.available (Known 12) and
        compare str.finished? false and
        compare (mat.size (str.read 9)) { rows = 1, columns = 9 } and
        compare str.position 9 and
        compare str.available (Known 3) and
        compare str.finished? false and
        compare (mat.size (str.read 9)) { rows = 1, columns = 3 } and
        compare str.position 12 and
        compare str.available (Known 0) and
        compare str.finished? true and
        ( str.close (); true )
),

"interpolated-sine": \(
    // Interpolating a sinusoid should give us a sinusoid, once we've
    // dropped the first few samples
    sinusoid = syn.sinusoid 8 2; // 2Hz sine sampled at 8Hz: [ 0, 1, 0, -1 ] etc
    output = resample.interpolated 2 sinusoid;
    \() (output.read 300);
    result = output.read 200;
    reference = syn.sinusoid 16 2;
    \() (reference.read 300);
    expected = reference.read 200;
    compareMatrices 1e-4 result expected;
),

"interpolated-rs-sine": \(
    // Just as above, but using resampledTo instead of interpolated
    sinusoid = syn.sinusoid 8 2; // 2Hz sine sampled at 8Hz: [ 0, 1, 0, -1 ] etc
    output = resample.resampledTo 16 sinusoid;
    \() (output.read 300);
    result = output.read 200;
    reference = syn.sinusoid 16 2;
    \() (reference.read 300);
    expected = reference.read 200;
    compareMatrices 1e-4 result expected;
),

"decimated-sine": \(
    // Decimating a sinusoid should give us a sinusoid, once we've
    // dropped the first few samples
    sinusoid = syn.sinusoid 32 1; // 1Hz sine sampled at 32Hz
    output = resample.decimated 2 sinusoid;
    \() (output.read 300);
    result = output.read 200;
    reference = syn.sinusoid 16 1;
    \() (reference.read 300);
    expected = reference.read 200;
    compareMatrices 1e-4 result expected;
),

"decimated-rs-sine": \(
    // Just as above, but using resampledTo instead of interpolated
    sinusoid = syn.sinusoid 32 1; // 1Hz sine sampled at 32Hz
    output = resample.resampledTo 16 sinusoid;
    \() (output.read 300);
    result = output.read 200;
    reference = syn.sinusoid 16 1;
    \() (reference.read 300);
    expected = reference.read 200;
    compareMatrices 1e-4 result expected;
),

"resample-spectrum": \(
    // Principle: Generate a wave with a lot of harmonics, resample by
    // some ratio, check that the resulting signal has substantially
    // the same magnitude spectrum as the original.
    
    // Read and window a chunk of signal from a stream
    sigOf str n =
       (sig = mat.getRow 0 (str.read n);
        win.windowed win.hann sig);

    // Return magnitude spectrum, running fft at n points and then
    // truncating to m
    specOf sig n m =
        vec.resizedTo m (vec.divideBy n (fft.realForwardMagnitude n sig));

    all id
       (map do { inrate, outrate }:
            freq = 500;
            forms = manip.duplicated 2 (waves.square inrate freq);
            inform = forms[0];
            outform = resample.resampledTo outrate forms[1];
            incount = inrate;
            outcount = outrate;
            speclen = (min incount outcount) / 2;
            // We don't compare bins within 4% of the Nyquist freq
            speclen = speclen - int (speclen / 25);
            inmag = specOf (sigOf inform incount) incount speclen;
            outmag = specOf (sigOf outform outcount) outcount speclen;
            compareMatrices 1e-7
               (mat.newRowVector outmag)
               (mat.newRowVector inmag);
        done [
            { inrate =  8000, outrate = 44100 },
            { inrate = 48000, outrate = 44100 },
            { inrate = 44100, outrate = 22050 },
            { inrate = 48000, outrate =  8000 },
        ]);
),

"interpolated-misc": \(
    // Interpolating any signal by N should give a signal in which
    // every Nth sample is the original signal
    data = vec.fromList [ 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 ];
    data = vec.concat [ data, vec.scaled (5/4) data, vec.scaled (3/4) data, data ];
    data = vec.concat [ data, data ];
    input = manip.withDuration (vec.length data) (syn.precalculatedMono 4 data);
    factor = 3;
    up = resample.interpolated factor input;
    result = mat.getRow 0 (up.read (factor * vec.length data));
    b = vec.fromList
       (map do i: vec.at result (i*factor) done [0..vec.length data - 1]);
    compareMatrices 1e-5 (mat.newRowVector b) (mat.newRowVector data);
),

"interpolated-rs-misc": \(
    // Just as above, but using resampledTo instead of interpolated
    data = vec.fromList [ 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 ];
    data = vec.concat [ data, vec.scaled (5/4) data, vec.scaled (3/4) data, data ];
    data = vec.concat [ data, data ];
    input = manip.withDuration (vec.length data) (syn.precalculatedMono 4 data);
    factor = 3;
    up = resample.resampledTo (factor * input.sampleRate) input;
    result = mat.getRow 0 (up.read (factor * vec.length data));
    b = vec.fromList
       (map do i: vec.at result (i*factor) done [0..vec.length data - 1]);
    compareMatrices 1e-5 (mat.newRowVector b) (mat.newRowVector data);
),

"down-up-2": \(
    testFrequencyIntegrity 441 44100 22050;
),

"down-up-16": \(
    testFrequencyIntegrity 300 48000 3000;
),

"up-down-2": \(
    testFrequencyIntegrity 441 44100 88200;
),

"up-down-16": \(
    testFrequencyIntegrity 300 3000 48000;
),

"same-rate": \(
    // Test that resample to the original rate leaves data unchanged
    data = vec.fromList [ 0, 0.1, -0.3, -0.4, -0.3, 0, 0.5, 0.2, 0.8, -0.1 ];
    data = vec.concat [ data, vec.scaled (5/4) data, vec.scaled (3/4) data, data ];
    data = vec.concat [ data, data ];
    input = syn.precalculatedMono 4 data;
    output = resample.resampledTo (input.sampleRate) input;
    b = output.read (vec.length data);
    compareMatrices 1e-8 b (mat.newRowVector data);
),

] is hash<string, () -> boolean>;