view yeti/test_frequency.yeti @ 68:b75c0eaaa6dd

Allow wrong max in alternate columns where expected max would be in lower octave (which is not present)

module test_frequency;

mat = load may.matrix;
vec = load may.vector;
win = load may.signal.window;
mm = load may.mathmisc;
cm = load may.matrix.complex;
syn = load may.stream.syntheticstream;
plot = load may.plot;

{ cqt } = load cqt;

// Test with a single windowed sinusoid, repeating at various frequencies

sinTestStream sampleRate duration signalFreq = // duration is in samples
   (sin = syn.sinusoid sampleRate signalFreq;
    chunk = mat.getRow 0 (sin.read duration);
    syn.precalculatedMono sampleRate (win.windowed win.hann chunk));

// We want to make a CQ transform spanning more than one octave, but
// not going all the way to fs/2 so we can test it also with
// frequencies above and below its extents

sampleRate = 100;

// fs/2 = 50 so 10->40 gives us 2 octaves
cqmin = 10;
cqmax = 40;
bpo = 4; // fairly arbitrary

testFreqs = map (* 5) [ 0..10 ];
duration = sampleRate * 2;

threshold = 0.08;

streamBuilder = sinTestStream sampleRate duration;

binForFreq f =
    mm.round (bpo * mm.log2 (f / cqmin)) - 1;

report message matrix =
   (eprintln message;
    eprintln "matrix is:";
    mat.eprint matrix;
    chart = plot.plot [Grid matrix];
    sleep 100;
    chart#dispose());

tests = mapIntoHash
    do f: "freq_\(f)" done
    do f: \(
        str = streamBuilder f;
        cq = cqt { maxFreq = cqmax, minFreq = cqmin, binsPerOctave = bpo } str;
        m = mat.concatHorizontal (map cm.magnitudes cq.output);
//    println "binFrequencies = \(cq.kernel.binFrequencies)";
//    println "binForFreq \(f) = \(binForFreq f)";
        var colno = 0;
        success = all id
           (map do c:
                // The test passes for this column if:
                //
                //  * the max bin is the expected one, or
                //
                //  * the expected max is out of range entirely (but
                //  we need to test _something_ in this case --
                //  what?), or
                //
                //  * all bins are below a threshold, or
                //
                //  * this is an odd column and the expected max is in
                //  the lower octave
                //
                // We should also check that all values in the lower
                // octave are zero for odd columns.
                //
                expected = binForFreq f;
                good =
                   (expected < 0 or expected >= vec.length c) or
                   ((colno % 2 == 1) and expected < (vec.length c / 2)) or
                   (vec.max c < threshold) or
                   (vec.maxindex c == binForFreq f);
                if (not good) then
                    report " * bad! maxindex \(vec.maxindex c) != expected \(binForFreq f) for freq \(f) in column \(colno) of \(mat.width m): column is \(vec.list c)" m;
                fi;
                colno := colno + 1;
                good;
            done (mat.asColumns m));
        success;
    ) done
    testFreqs;

tests is hash<string, () -> boolean>