Mercurial > hg > constant-q-cpp

changeset 22:701900c371b0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Calculate (so far unnormalised) CQ kernel
author Chris Cannam <c.cannam@qmul.ac.uk>
date Fri, 01 Nov 2013 16:13:22 +0000
parents 5787785f4560
children 4f2f28d5df58
files cpp-qm-dsp/CQKernel.cpp cpp-qm-dsp/CQKernel.h cpp-qm-dsp/Makefile cpp-qm-dsp/test.cpp yeti/cqt.yeti
diffstat 5 files changed, 239 insertions(+), 4 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQKernel.cpp	Fri Nov 01 16:13:22 2013 +0000
@@ -0,0 +1,138 @@
+
+#include "CQKernel.h"
+
+#include "qm-dsp/maths/MathUtilities.h"
+#include "qm-dsp/dsp/transforms/FFT.h"
+#include "qm-dsp/base/Window.h"
+
+#include <cmath>
+#include <cassert>
+#include <vector>
+
+using std::vector;
+
+CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave)
+{
+    m_p.sampleRate = sampleRate;
+    m_p.maxFrequency = maxFreq;
+    m_p.binsPerOctave = binsPerOctave;
+    generateKernel();
+}
+
+CQKernel::~CQKernel()
+{
+    delete m_fft;
+}
+
+void
+CQKernel::generateKernel()
+{
+    double q = 1;
+    double atomHopFactor = 0.25;
+    double thresh = 0.0005;
+
+    double bpo = m_p.binsPerOctave;
+
+    m_p.minFrequency = (m_p.maxFrequency / 2) * pow(2, 1.0/bpo);
+    m_p.Q = q / (pow(2, 1.0/bpo) - 1.0);
+
+    double maxNK = round(m_p.Q * m_p.sampleRate / m_p.minFrequency);
+    double minNK = round
+	(m_p.Q * m_p.sampleRate /
+	 (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo)));
+
+    m_p.atomSpacing = round(minNK * atomHopFactor);
+    m_p.firstCentre = m_p.atomSpacing * ceil(ceil(maxNK / 2.0) / m_p.atomSpacing);
+    m_p.fftSize = MathUtilities::nextPowerOfTwo
+	(m_p.firstCentre + ceil(maxNK / 2.0));
+
+    m_p.atomsPerFrame = floor
+	(1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing);
+
+    int lastCentre = m_p.firstCentre + (m_p.atomsPerFrame - 1) * m_p.atomSpacing;
+
+    m_p.fftHop = (lastCentre + m_p.atomSpacing) - m_p.firstCentre;
+
+    m_fft = new FFT(m_p.fftSize);
+
+    for (int k = 1; k <= m_p.binsPerOctave; ++k) {
+	
+	int nk = round(m_p.Q * m_p.sampleRate /
+		       (m_p.minFrequency * pow(2, ((k-1.0) / bpo))));
+
+	// The MATLAB version uses a symmetric window, but our windows
+	// are periodic. A symmetric window of size N is a periodic
+	// one of size N-1 with the first element stuck on the end
+	Window<double> w(BlackmanHarrisWindow, nk-1);
+	vector<double> win = w.getWindowData();
+	win.push_back(win[0]);
+
+	for (int i = 0; i < (int)win.size(); ++i) {
+	    win[i] = sqrt(win[i]) / nk;
+	}
+	
+	double fk = m_p.minFrequency * pow(2, ((k-1.0) * bpo));
+
+	vector<double> reals, imags;
+	
+	for (int i = 0; i < nk; ++i) {
+	    double arg = 2.0 * M_PI * fk * i / m_p.sampleRate;
+	    reals[i] = win[i] * cos(arg);
+	    imags[i] = win[i] * sin(arg);
+	}
+
+	int atomOffset = m_p.firstCentre - ceil(nk/2);
+
+	for (int i = 0; i < m_p.atomsPerFrame; ++i) {
+
+	    int shift = atomOffset + (i * m_p.atomSpacing);
+
+	    vector<double> rin(m_p.fftSize, 0.0);
+	    vector<double> iin(m_p.fftSize, 0.0);
+
+	    for (int j = 0; j < nk; ++j) {
+		rin[j + shift] = reals[j];
+		iin[i + shift] = imags[j];
+	    }
+
+	    vector<double> rout(m_p.fftSize, 0.0);
+	    vector<double> iout(m_p.fftSize, 0.0);
+		
+	    m_fft->process(false,
+			   rin.data(), iin.data(),
+			   rout.data(), iout.data());
+
+	    int firstNZ = -1, lastNZ = -1;
+
+	    for (int j = 0; j < m_p.fftSize; ++j) {
+		if (sqrt(rout[i] * rout[i] + iout[i] * iout[i]) >= thresh) {
+		    lastNZ = j;
+		    if (firstNZ < 0) firstNZ = j;
+		}
+	    }
+
+	    vector<double> rnz, inz;
+
+	    if (firstNZ >= 0) {
+		for (int j = firstNZ; j <= lastNZ; ++j) {
+		    rnz.push_back(rout[j] / m_p.fftSize);
+		    inz.push_back(iout[j] / m_p.fftSize);
+		}
+		m_kernel.offsets.push_back(firstNZ);
+	    } else {
+		m_kernel.offsets.push_back(0);
+	    }
+
+	    m_kernel.real.push_back(rnz);
+	    m_kernel.imag.push_back(inz);
+	}
+    }
+
+    assert((int)m_kernel.offsets.size() == m_p.binsPerOctave);
+    assert((int)m_kernel.real.size() == m_p.binsPerOctave);
+    assert((int)m_kernel.imag.size() == m_p.binsPerOctave);
+
+    //!!! and normalise
+}
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/CQKernel.h	Fri Nov 01 16:13:22 2013 +0000
@@ -0,0 +1,44 @@
+
+#ifndef CQ_KERNEL_H
+#define CQ_KERNEL_H
+
+#include <vector>
+
+class FFT;
+
+class CQKernel
+{
+public:
+    CQKernel(double sampleRate, double maxFreq, int binsPerOctave);
+    ~CQKernel();
+    
+    struct Properties {
+	double sampleRate;
+	double maxFrequency;
+	double minFrequency;
+	int binsPerOctave;
+	int fftSize;
+	int fftHop;
+	int atomsPerFrame;
+	int atomSpacing;
+	int firstCentre;
+	double Q;
+    };
+
+    Properties getProperties() const { return m_p; }
+
+private:
+    Properties m_p;
+    FFT *m_fft;
+
+    struct KernelMatrix {
+	std::vector<int> offsets;
+	std::vector<std::vector<double> > real;
+	std::vector<std::vector<double> > imag;
+    };
+    KernelMatrix m_kernel;
+
+    void generateKernel();
+};
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/Makefile	Fri Nov 01 16:13:22 2013 +0000
@@ -0,0 +1,30 @@
+
+DEFINES := -DUSE_PTHREADS
+
+CFLAGS := -I../.. $(CFLAGS) $(DEFINES)
+
+CXXFLAGS := -I../.. -Wall -g $(CXXFLAGS) $(DEFINES)
+#CXXFLAGS := -I../.. -Wall -O3 -ffast-math -ftree-vectorize $(CXXFLAGS) $(DEFINES)
+
+LDFLAGS := $(LDFLAGS)
+
+#VG	:= valgrind
+
+LIBS 		:= ../../qm-dsp/libqm-dsp.a -lpthread
+PROGRAM_LIBS	:= -lsndfile
+
+SOURCES	:= CQKernel.cpp test.cpp
+
+OBJECTS := $(SOURCES:.cpp=.o)
+OBJECTS := $(OBJECTS:.c=.o)
+
+PROGRAM	:= test
+
+all: $(PROGRAM)
+
+test:	$(OBJECTS)
+	$(CXX) -o $@ $^ $(LDFLAGS) $(LIBS) $(PROGRAM_LIBS)
+
+clean: 
+	rm -f *.o 
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/test.cpp	Fri Nov 01 16:13:22 2013 +0000
@@ -0,0 +1,13 @@
+
+#include "CQKernel.h"
+
+#include <iostream>
+
+int main(int argc, char **argv)
+{
+    CQKernel k(96000, 48000, 24);
+
+    std::cerr << "Q = " << k.getProperties().Q << std::endl;
+
+}
+
--- a/yeti/cqt.yeti	Wed Oct 30 18:26:04 2013 +0000
+++ b/yeti/cqt.yeti	Fri Nov 01 16:13:22 2013 +0000
@@ -12,13 +12,14 @@
 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 = 40;
+    minFreq = 50;
     binsPerOctave = 24;
 
 eprintln "Here";
@@ -109,6 +110,9 @@
 
     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]);
@@ -140,7 +144,7 @@
             atomNo = int (col / colsPerAtom);
             atomOffset = col % colsPerAtom;
 
-            if atomOffset == 0 and atomNo < length bits[oct] then
+            if /*!!! atomOffset == 0 and */ atomNo < length bits[oct] then
                 bits[oct][atomNo][binNo];
             else
                 cplx.zero
@@ -190,14 +194,20 @@
 //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";
 
-for (cqt testStream) do c:
+cq = cqt testStream;
+
+for cq do c:
     mm = cm.magnitudes c;
     for (mat.asColumns mm) (println . strJoin "," . vec.list);
 done;
 
+bigM = mat.concat (Horizontal ()) (map cm.magnitudes cq);
+
+//\() (plot.plot [Contour bigM]);
+
 ()
 
-