c@69: /*
c@69:     Constant-Q library
c@69:     Copyright (c) 2013-2014 Queen Mary, University of London
c@69: 
c@69:     Permission is hereby granted, free of charge, to any person
c@69:     obtaining a copy of this software and associated documentation
c@69:     files (the "Software"), to deal in the Software without
c@69:     restriction, including without limitation the rights to use, copy,
c@69:     modify, merge, publish, distribute, sublicense, and/or sell copies
c@69:     of the Software, and to permit persons to whom the Software is
c@69:     furnished to do so, subject to the following conditions:
c@69: 
c@69:     The above copyright notice and this permission notice shall be
c@69:     included in all copies or substantial portions of the Software.
c@69: 
c@69:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
c@69:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
c@69:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
c@69:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
c@69:     CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
c@69:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
c@69:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
c@69: 
c@69:     Except as contained in this notice, the names of the Centre for
c@69:     Digital Music; Queen Mary, University of London; and Chris Cannam
c@69:     shall not be used in advertising or otherwise to promote the sale,
c@69:     use or other dealings in this Software without prior written
c@69:     authorization.
c@69: */
c@65: 
c@65: module test_frequency;
c@65: 
c@65: mat = load may.matrix;
c@65: vec = load may.vector;
c@65: win = load may.signal.window;
c@65: mm = load may.mathmisc;
c@65: cm = load may.matrix.complex;
c@65: syn = load may.stream.syntheticstream;
c@67: plot = load may.plot;
c@65: 
c@81: { compare } = load may.test;
c@81: 
c@65: { cqt } = load cqt;
c@65: 
c@65: // Test with a single windowed sinusoid, repeating at various frequencies
c@65: 
c@65: sinTestStream sampleRate duration signalFreq = // duration is in samples
c@65:    (sin = syn.sinusoid sampleRate signalFreq;
c@65:     chunk = mat.getRow 0 (sin.read duration);
c@65:     syn.precalculatedMono sampleRate (win.windowed win.hann chunk));
c@65: 
c@65: // We want to make a CQ transform spanning more than one octave, but
c@65: // not going all the way to fs/2 so we can test it also with
c@65: // frequencies above and below its extents
c@65: 
c@65: sampleRate = 100;
c@65: 
c@65: // fs/2 = 50 so 10->40 gives us 2 octaves
c@65: cqmin = 10;
c@65: cqmax = 40;
c@65: bpo = 4; // fairly arbitrary
c@65: 
c@65: testFreqs = map (* 5) [ 0..10 ];
c@65: duration = sampleRate * 2;
c@65: 
c@68: threshold = 0.08;
c@67: 
c@65: streamBuilder = sinTestStream sampleRate duration;
c@65: 
c@65: binForFreq f =
c@65:     mm.round (bpo * mm.log2 (f / cqmin)) - 1;
c@65: 
c@68: report message matrix =
c@68:    (eprintln message;
c@68:     eprintln "matrix is:";
c@81:     mat.eprint matrix);
c@81: //    chart = plot.plot [Grid matrix];
c@81: //    sleep 100;
c@81: //    chart#dispose());
c@68: 
c@65: tests = mapIntoHash
c@65:     do f: "freq_\(f)" done
c@65:     do f: \(
c@65:         str = streamBuilder f;
c@65:         cq = cqt { maxFreq = cqmax, minFreq = cqmin, binsPerOctave = bpo } str;
c@81:         spec = cq.cqSpectrogram;
c@81:         rightSize = all id
c@81:            (map do s:
c@81:             compare (mat.size s) { 
c@81:                 rows = cq.kernel.binsPerOctave * cq.octaves,
c@81:                 columns = cq.kernel.atomsPerFrame * mm.pow 2 (cq.octaves - 1)
c@81:             }
c@81:             done spec);
c@81:         m = mat.concatHorizontal spec;
c@65: //    println "binFrequencies = \(cq.kernel.binFrequencies)";
c@65: //    println "binForFreq \(f) = \(binForFreq f)";
c@67:         var colno = 0;
c@65:         success = all id
c@81:            (rightSize :: map do c:
c@68:                 // The test passes for this column if:
c@68:                 //
c@68:                 //  * the max bin is the expected one, or
c@68:                 //
c@68:                 //  * the expected max is out of range entirely (but
c@68:                 //  we need to test _something_ in this case --
c@68:                 //  what?), or
c@68:                 //
c@68:                 //  * all bins are below a threshold, or
c@68:                 //
c@68:                 //  * this is an odd column and the expected max is in
c@68:                 //  the lower octave
c@68:                 //
c@68:                 // We should also check that all values in the lower
c@68:                 // octave are zero for odd columns.
c@68:                 //
c@65:                 expected = binForFreq f;
c@65:                 good =
c@65:                    (expected < 0 or expected >= vec.length c) or
c@68:                    ((colno % 2 == 1) and expected < (vec.length c / 2)) or
c@67:                    (vec.max c < threshold) or
c@65:                    (vec.maxindex c == binForFreq f);
c@65:                 if (not good) then
c@68:                     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;
c@65:                 fi;
c@67:                 colno := colno + 1;
c@65:                 good;
c@65:             done (mat.asColumns m));
c@65:         success;
c@65:     ) done
c@65:     testFreqs;
c@65: 
c@65: tests is hash<string, () -> boolean>