view yeti/test_frequency.yeti @ 67:d6dd6e1cc00e

Plot failing test matrices

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.05;

streamBuilder = sinTestStream sampleRate duration;

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

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:
                // passes test if the correct max bin, or the expected max
                // is out of range, or if all bins are below a threshold
                expected = binForFreq f;
                good =
                   (expected < 0 or expected >= vec.length c) or
                   (vec.max c < threshold) or
                   (vec.maxindex c == binForFreq f);
                if (not good) then
                    println " * bad! maxindex \(vec.maxindex c) != expected \(binForFreq f) for freq \(f) in column \(colno) of \(mat.width m): column is \(vec.list c)";
                println "matrix is:";
                mat.print m;
                chart = plot.plot [Grid m];
                sleep 100;
                chart#dispose();
                else 
                    print "✓"; 
                fi;
                colno := colno + 1;
                good;
            done (mat.asColumns m));
        success;
    ) done
    testFreqs;

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