c@10: 
c@19: program cqt;
c@10: 
c@10: cqtkernel = load cqtkernel;
c@10: resample = load may.stream.resample;
c@10: manipulate = load may.stream.manipulate;
c@10: syn = load may.stream.syntheticstream;
c@10: cm = load may.matrix.complex;
c@10: mat = load may.matrix;
c@10: framer = load may.stream.framer;
c@10: cplx = load may.complex;
c@10: fft = load may.transform.fft;
c@10: vec = load may.vector;
c@17: af = load may.stream.audiofile;
c@10: 
c@10: { pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc;
c@10: 
c@10: cqt str =
c@10:    (sampleRate = str.sampleRate;
c@10:     maxFreq = sampleRate/2;
c@10:     minFreq = 40;
c@10:     binsPerOctave = 24;
c@10: 
c@16: eprintln "Here";
c@10: 
c@10:     octaves = ceil (log2 (maxFreq / minFreq));
c@10: 
c@16: eprintln "Here: about to calculate stuff with \(octaves)";
c@10: 
c@10:     actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave));
c@10: 
c@16:     eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave)";
c@10: 
c@10:     kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave };
c@10: 
c@16:     eprintln "atomsPerFrame = \(kdata.atomsPerFrame)";
c@11: 
c@10:     streams = manipulate.duplicated octaves str;
c@10: 
c@10:     //!!! can't be right!
c@10:     kernel = cm.transposed (cm.conjugateTransposed kdata.kernel);
c@10: 
c@16:     eprintln "have kernel";
c@10: 
c@10:     fftFunc = fft.forward kdata.fftSize;
c@10: 
c@10:     cqblocks =
c@10:         map do octave:
c@10:             frames = framer.monoFrames //!!! mono for now
c@10:                 { framesize = kdata.fftSize, hop = kdata.fftHop }
c@10:                 (resample.decimated (pow 2 octave) streams[octave]);
c@10:             map do frame:
c@10:                 freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize));
c@10:                 cm.product kernel (cm.newComplexColumnVector freq);
c@10:             done frames;
c@10:         done [0..octaves-1];
c@10: 
c@13:     // The cqblocks list is a list<list<matrix>>. Each top-level list
c@11:     // corresponds to an octave, from highest to lowest, each having
c@11:     // twice as many elements in its list as the next octave. The
c@11:     // sub-lists are sampled in time with an effective spacing of
c@11:     // fftSize * 2^(octave-1) audio frames, and the matrices are row
c@11:     // vectors with atomsPerFrame * binsPerOctave complex elements.
c@13:     //
c@13:     // ***
c@13:     // 
c@13:     // In a typical constant-Q structure, each (2^(octaves-1) *
c@13:     // fftHop) input frames gives us an output structure conceptually
c@13:     // like this:
c@10:     //
c@10:     // [][][][][][][][]   <- fftHop frames per highest-octave output value
c@10:     // [][][][][][][][]      layered as many times as binsPerOctave (here 2)
c@10:     // [--][--][--][--]   <- fftHop*2 frames for the next lower octave
c@10:     // [--][--][--][--]      etc
c@10:     // [------][------]
c@10:     // [------][------]
c@10:     // [--------------]
c@10:     // [--------------]
c@10:     //
c@13:     // ***
c@13:     //
c@13:     // But the kernel we're using here has more than one temporally
c@13:     // spaced atom; each individual cell is a row vector with
c@13:     // atomsPerFrame * binsPerOctave elements, but that actually
c@13:     // represents a rectangular matrix of result cells with width
c@13:     // atomsPerFrame and height binsPerOctave. The columns of this
c@13:     // matrix (the atoms) then need to be spaced by 2^(octave-1)
c@13:     // relative to those from the highest octave.
c@10: 
c@15:     // Reshape each row vector into the appropriate rectangular matrix
c@19: /*
c@17:     cqblocks = array (map do octlist:
c@14:         map do rv:
c@15:             cm.generate do row col:
c@15:                 cm.at rv ((row * kdata.atomsPerFrame) + col) 0
c@15:             done {
c@15:                 rows = kdata.binsPerOctave,
c@15:                 columns = kdata.atomsPerFrame
c@15:             }
c@14:         done octlist
c@17:     done cqblocks);
c@19: */
c@19: 
c@19: //!!! how then do we arrange to drop a certain number of atoms (rather
c@19: //than of atoms+bins chunks)?
c@14: 
c@17:     assembleBlock bits =
c@19:        (eprintln "assembleBlock: structure of bits is:";
c@19:         eprintln (map length bits);
c@19: 
c@19:         rows = octaves * kdata.binsPerOctave;
c@19:         columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
c@19: 
c@18:         cm.generate do row col:
c@19: 
c@19:             // bits structure: [1,2,4,8,...]
c@19: 
c@19:             // each elt of bits is a list of the chunks that should
c@19:             // make up this block in that octave (lowest octave first)
c@19: 
c@19:             // each chunk has atomsPerFrame * binsPerOctave elts in it
c@19: 
c@19:             // row is disposed with 0 at the top, highest octave (in
c@19:             // both pitch and index into bits structure)
c@19: 
c@19: //            oct = octaves - int (row / binsPerOctave) - 1;
c@18:             oct = int (row / binsPerOctave);
c@19:             binNo = row % kdata.binsPerOctave;
c@19:             chunks = pow 2 oct;
c@19:             colsPerChunk = int (columns / chunks);
c@19:             colsPerAtom = int (colsPerChunk / kdata.atomsPerFrame);
c@19:             chunkNo = int (col / colsPerChunk);
c@19:             atomNo = int ((col % colsPerChunk) / colsPerAtom);
c@18: 
c@19: //            eprintln "row \(row) of \(rows), col \(col) of \(columns): oct \(oct), bin \(binNo), chunk \(chunkNo) of \(chunks), atom \(atomNo) of \(kdata.atomsPerFrame)";
c@19: 
c@19:             cm.at bits[oct][chunkNo] (binNo * kdata.atomsPerFrame + atomNo) 0;
c@19: 
c@19:         done { rows, columns };
c@19:         );
c@15: 
c@17:     processOctaveLists octs =
c@17:         case octs[0] of
c@17:         block::rest:
c@19:            (toAssemble = array 
c@19:                (map do oct:
c@17:                     n = pow 2 oct;
c@17:                     if not empty? octs[oct] then
c@19:                         forBlock = array (take n octs[oct]);
c@17:                         octs[oct] := drop n octs[oct];
c@17:                         forBlock
c@17:                     else
c@19:                         array []
c@17:                     fi
c@19:                 done (keys octs));
c@17:             assembleBlock toAssemble :. \(processOctaveLists octs));
c@17:          _: []
c@15:         esac;
c@15: 
c@19: eprintln "cqblocks has \(length cqblocks) entries";
c@15: 
c@17:     octaveLists = [:];
c@19: 
c@19:     cqblocks = array cqblocks;
c@17:     for [1..octaves] do oct:
c@17:         octaveLists[octaves - oct] := cqblocks[oct-1];
c@17:     done;
c@17: /*
c@17:     \() (map2 do octlist octave:
c@17: println "oct \(octaves) - \(octave) = \(octaves - octave)";
c@17:              octaveLists[octaves - octave] := octlist 
c@17:          done cqblocks [1..octaves]);
c@17: */
c@19: eprintln "octaveLists keys are: \(keys octaveLists)";
c@17:     
c@17:     processOctaveLists octaveLists;
c@15: 
c@10:     );
c@10: 
c@17: //testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500);
c@17: //testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4);
c@17: testStream = af.open "sweep.wav";
c@10: 
c@16: eprintln "have test stream";
c@10: 
c@16: c = cqt testStream;
c@10: 
c@19: //m = take 1 (drop 2 c);
c@19: 
c@19: //thing = take 50 (drop 200 c);
c@19: 
c@19: //m = cm.newComplexMatrix (ColumnMajor ()) thing;
c@19: mm = cm.magnitudes (head c);
c@10: 
c@16: for (mat.asColumns mm) (println . strJoin "," . vec.list);
c@10: 
c@16: ()
c@16: 
c@16: