view yeti/cqt.yeti @ 31:01a3e110bf8d

Don't forget top octave in buffer list

program cqt;

cqtkernel = load cqtkernel;
resample = load may.stream.resample;
manipulate = load may.stream.manipulate;
syn = load may.stream.syntheticstream;
cm = load may.matrix.complex;
mat = load may.matrix;
framer = load may.stream.framer;
cplx = load may.complex;
fft = load may.transform.fft;
vec = load may.vector;
af = load may.stream.audiofile;
plot = load may.plot;

{ pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc;

cqt str =
   (sampleRate = str.sampleRate;
    maxFreq = sampleRate/2;
    minFreq = 50;
    binsPerOctave = 24;

eprintln "Here";

    octaves = ceil (log2 (maxFreq / minFreq));

eprintln "Here: about to calculate stuff with \(octaves)";

    actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave));

    eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave)";

    kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave };

    eprintln "atomsPerFrame = \(kdata.atomsPerFrame)";

    streams = manipulate.duplicated octaves str;

    //!!! can't be right!
    kernel = cm.transposed (cm.conjugateTransposed kdata.kernel);

    eprintln "have kernel";

    fftFunc = fft.forward kdata.fftSize;

    cqblocks =
        map do octave:
            frames = framer.monoFrames //!!! mono for now
                { framesize = kdata.fftSize, hop = kdata.fftHop }
                (resample.decimated (pow 2 octave) streams[octave]);
            map do frame:
                freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize));
                cm.product kernel (cm.newComplexColumnVector freq);
            done frames;
        done [0..octaves-1];

    // The cqblocks list is a list<list<matrix>>. Each top-level list
    // corresponds to an octave, from highest to lowest, each having
    // twice as many elements in its list as the next octave. The
    // sub-lists are sampled in time with an effective spacing of
    // fftSize * 2^(octave-1) audio frames, and the matrices are row
    // vectors with atomsPerFrame * binsPerOctave complex elements.
    //
    // ***
    // 
    // In a typical constant-Q structure, each (2^(octaves-1) *
    // fftHop) input frames gives us an output structure conceptually
    // like this:
    //
    // [][][][][][][][]   <- fftHop frames per highest-octave output value
    // [][][][][][][][]      layered as many times as binsPerOctave (here 2)
    // [--][--][--][--]   <- fftHop*2 frames for the next lower octave
    // [--][--][--][--]      etc
    // [------][------]
    // [------][------]
    // [--------------]
    // [--------------]
    //
    // ***
    //
    // But the kernel we're using here has more than one temporally
    // spaced atom; each individual cell is a row vector with
    // atomsPerFrame * binsPerOctave elements, but that actually
    // represents a rectangular matrix of result cells with width
    // atomsPerFrame and height binsPerOctave. The columns of this
    // matrix (the atoms) then need to be spaced by 2^(octave-1)
    // relative to those from the highest octave.

    // Reshape each row vector into the appropriate rectangular matrix
    // and split into single-atom columns

    emptyHops = kdata.firstCentre / kdata.atomSpacing;
    maxDrop = emptyHops * (pow 2 (octaves-1)) - emptyHops;
    eprintln "maxDrop = \(maxDrop)";

    cqblocks = map do octlist:
        concat
           (map do rv:
                cm.asColumns
                   (cm.generate do row col:
                        cm.at rv ((row * kdata.atomsPerFrame) + col) 0
                    done {
                        rows = kdata.binsPerOctave,
                        columns = kdata.atomsPerFrame
                    })
            done octlist)
    done cqblocks;

    cqblocks = array (map2 do octlist octave:
        d = emptyHops * (pow 2 (octaves-octave)) - emptyHops;

        d = 0; //!!!

        eprintln "dropping \(d)";
        drop d octlist;
    done cqblocks [1..octaves]);

    assembleBlock bits =
       (eprintln "assembleBlock: structure of bits is:";
        eprintln (map length bits);

        rows = octaves * kdata.binsPerOctave;
        columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;

        cm.generate do row col:

            // bits structure: [1,2,4,8,...]

            // each elt of bits is a list of the chunks that should
            // make up this block in that octave (lowest octave first)

            // each chunk has atomsPerFrame * binsPerOctave elts in it

            // row is disposed with 0 at the top, highest octave (in
            // both pitch and index into bits structure)

            oct = int (row / binsPerOctave);
            binNo = row % kdata.binsPerOctave;

            chunks = pow 2 oct;
            colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
            atomNo = int (col / colsPerAtom);
            atomOffset = col % colsPerAtom;

            if /*!!! atomOffset == 0 and */ atomNo < length bits[oct] then
                bits[oct][atomNo][binNo];
            else
                cplx.zero
            fi;

        done { rows, columns };
        );

    processOctaveLists octs =
        case octs[0] of
        block::rest:
           (toAssemble = array 
               (map do oct:
                    n = kdata.atomsPerFrame * pow 2 oct;
                    if not empty? octs[oct] then
                        forBlock = array (take n octs[oct]);
                        octs[oct] := drop n octs[oct];
                        forBlock
                    else
                        array []
                    fi
                done (keys octs));
            assembleBlock toAssemble :. \(processOctaveLists octs));
         _: []
        esac;

eprintln "cqblocks has \(length cqblocks) entries";

    octaveLists = [:];

    cqblocks = array cqblocks;
    for [1..octaves] do oct:
        octaveLists[octaves - oct] := cqblocks[oct-1];
    done;
/*
    \() (map2 do octlist octave:
println "oct \(octaves) - \(octave) = \(octaves - octave)";
             octaveLists[octaves - octave] := octlist 
         done cqblocks [1..octaves]);
*/
eprintln "octaveLists keys are: \(keys octaveLists)";
    
    processOctaveLists octaveLists;

    );

//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500);
//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4);
testStream = af.open "sweep-48000.wav";
//testStream = af.open "sweep.wav";

eprintln "have test stream";

cq = cqt testStream;

for cq do c:
    mm = cm.magnitudes c;
    for (mat.asColumns mm) (println . strJoin "," . vec.list);
done;

bigM = mat.concatHorizontal (map cm.magnitudes cq);

//\() (plot.plot [Contour bigM]);

()