changeset 94:fa1709ed4a3c

Inverse runs; now merely needs to run correctly
author Chris Cannam <c.cannam@qmul.ac.uk>
date Fri, 09 May 2014 15:16:48 +0100
parents 908be1d06bd2
children 8d04787ca7dc
files Makefile.inc cpp-qm-dsp/CQInverse.cpp cpp-qm-dsp/CQInverse.h cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/processfile.cpp cpp-qm-dsp/test.cpp vamp/CQVamp.cpp
diffstat 7 files changed, 262 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile.inc	Fri May 09 11:41:44 2014 +0100
+++ b/Makefile.inc	Fri May 09 15:16:48 2014 +0100
@@ -16,9 +16,11 @@
 LDFLAGS := $(LDFLAGS) 
 PLUGIN_LDFLAGS := $(LDFLAGS) $(PLUGIN_LDFLAGS)
 TEST_LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework
+PF_LDFLAGS := $(LDFLAGS) -lsndfile
 
 PLUGIN	:= cqvamp$(PLUGIN_EXT)
 TEST	:= $(LIB_DIR)/test
+PF	:= $(LIB_DIR)/processfile
 
 LIB_HEADERS := $(LIB_DIR)/CQBase.h $(LIB_DIR)/CQKernel.h $(LIB_DIR)/ConstantQ.h $(LIB_DIR)/CQSpectrogram.h $(LIB_DIR)/CQInverse.h
 LIB_SOURCES := $(LIB_DIR)/CQKernel.cpp $(LIB_DIR)/ConstantQ.cpp $(LIB_DIR)/CQSpectrogram.cpp $(LIB_DIR)/CQInverse.cpp
@@ -30,12 +32,15 @@
 SOURCES	     := $(LIB_SOURCES) $(VAMP_SOURCES)
 OBJECTS	     := $(SOURCES:.cpp=.o)
 
-TEST_SOURCES := $(LIB_DIR)/test.cpp $(SOURCES)
-TEST_OBJECTS := $(TEST_SOURCES:.cpp=.o)
+TEST_SOURCES := $(LIB_DIR)/test.cpp 
+TEST_OBJECTS := $(TEST_SOURCES:.cpp=.o) $(OBJECTS)
+
+PF_SOURCES := $(LIB_DIR)/processfile.cpp
+PF_OBJECTS := $(PF_SOURCES:.cpp=.o) $(OBJECTS)
 
 LIBS	:= $(QMDSP_DIR)/libqm-dsp.a $(VAMPSDK_DIR)/libvamp-sdk.a -lpthread
 
-all: $(PLUGIN) $(TEST)
+all: $(PLUGIN) $(TEST) $(PF)
 
 $(PLUGIN):	$(OBJECTS)
 	$(CXX) -o $@ $^ $(LIBS) $(PLUGIN_LDFLAGS)
@@ -43,14 +48,17 @@
 $(TEST):	$(TEST_OBJECTS)
 	$(CXX) -o $@ $^ $(LIBS) $(TEST_LDFLAGS)
 
+$(PF):	$(PF_OBJECTS)
+	$(CXX) -o $@ $^ $(LIBS) $(PF_LDFLAGS)
+
 clean:		
-	rm -f $(OBJECTS)
+	rm -f $(OBJECTS) $(TEST_OBJECTS) $(PF_OBJECTS)
 
 distclean:	clean
 	rm -f $(PLUGIN)
 
 depend:
-	makedepend -Y -fMakefile.inc $(SOURCES) $(HEADERS)
+	makedepend -Y -fMakefile.inc $(SOURCES) $(TEST_SOURCES) $(PF_SOURCES) $(HEADERS)
 
 # DO NOT DELETE
 
@@ -59,16 +67,21 @@
 cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/CQKernel.h
 cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h
 cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
-cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQInverse.h cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQInverse.h cpp-qm-dsp/CQBase.h
+cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQKernel.h
 vamp/CQVamp.o: vamp/CQVamp.h cpp-qm-dsp/CQSpectrogram.h
 vamp/CQVamp.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
 vamp/CQVamp.o: cpp-qm-dsp/CQKernel.h
 vamp/libmain.o: vamp/CQVamp.h cpp-qm-dsp/CQSpectrogram.h
 vamp/libmain.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
 vamp/libmain.o: cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/test.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h
+cpp-qm-dsp/test.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/processfile.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
+cpp-qm-dsp/processfile.o: cpp-qm-dsp/CQKernel.h cpp-qm-dsp/CQInverse.h
 cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
 cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h
 cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQKernel.h
-cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQKernel.h
+cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
 vamp/CQVamp.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h
 vamp/CQVamp.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h
--- a/cpp-qm-dsp/CQInverse.cpp	Fri May 09 11:41:44 2014 +0100
+++ b/cpp-qm-dsp/CQInverse.cpp	Fri May 09 15:16:48 2014 +0100
@@ -127,13 +127,15 @@
     m_outputLatency = maxLatency; //!!! for now
 
     for (int i = 0; i < m_octaves; ++i) {
-
 	// Calculate the difference between the total latency applied
 	// across all octaves, and the existing latency due to the
 	// upsampler for this octave
+        m_buffers.push_back(RealSequence(m_outputLatency - latencies[i], 0.0));
+    }
 
-        m_buffers.push_back
-            (vector<double>(m_outputLatency - latencies[i], 0.0));
+    for (int i = 0; i < m_octaves; ++i) {
+        // Fixed-size buffer for IFFT overlap-add
+        m_olaBufs.push_back(RealSequence(m_p.fftSize, 0.0));
     }
 
     m_fft = new FFTReal(m_p.fftSize);
@@ -148,10 +150,16 @@
     // be the case for data that came directly from our ConstantQ
     // implementation.
 
+    int widthProvided = block.size();
+
+    if (widthProvided == 0) {
+        return drawFromBuffers();
+    }
+
+    cerr << "process: given " << widthProvided << " cols" << endl;
+
     int blockWidth = m_p.atomsPerFrame * int(pow(2, m_octaves - 1));
 
-    int widthProvided = block.size();
-
     if (widthProvided % blockWidth != 0) {
         cerr << "ERROR: CQInverse::process: Input block size ("
              << widthProvided
@@ -180,14 +188,14 @@
     // 4. Overlap-add each octave's resynthesised blocks (unwindowed)
     //
     // 5. Resample each octave's overlap-add stream to the original
-    // rate, and sum.
+    // rate
+    //
+    // 6. Sum the resampled streams and return
     
-    // We will, for now, do all but the last step in sequence, one
-    // octave at a time, and push the results to m_buffers for summing
-    // and return.
-
     for (int i = 0; i < m_octaves; ++i) {
         
+        // Step 1
+
         ComplexBlock oct;
 
         for (int j = 0; j < widthProvided; ++j) {
@@ -200,10 +208,51 @@
             oct.push_back(col);
         }
 
+        // Steps 2, 3, 4, 5
         processOctave(i, oct);
     }
-            
-#error need to return something
+    
+    // Step 6
+    return drawFromBuffers();
+}
+
+CQInverse::RealSequence
+CQInverse::drawFromBuffers()
+{
+    // 6. Sum the resampled streams and return
+
+    int available = 0;
+
+    for (int i = 0; i < m_octaves; ++i) {
+        if (i == 0 || int(m_buffers[i].size()) < available) {
+            available = m_buffers[i].size();
+        }
+    }
+
+    RealSequence result(available, 0);
+
+    if (available == 0) {
+        return result;
+    }
+
+    for (int i = 0; i < m_octaves; ++i) {
+        for (int j = 0; j < available; ++j) {
+            result[j] += m_buffers[i][j];
+        }
+        m_buffers[i] = RealSequence(m_buffers[i].begin() + available,
+                                    m_buffers[i].end());
+    }
+
+    cerr << "process: returning " << available << endl;
+
+    return result;
+}
+
+CQInverse::RealSequence
+CQInverse::getRemainingOutput()
+{
+    //!!! for now -- may have to prompt the resamplers?
+    return drawFromBuffers();
 }
 
 void
@@ -245,6 +294,16 @@
     // kernel (which is the conjugate of the forward kernel) and
     // perform an inverse FFT
 
+    if ((int)column.size() != m_p.atomsPerFrame * m_binsPerOctave) {
+        cerr << "ERROR: CQInverse::processOctaveColumn: Height of column ("
+             << column.size() << ") in octave " << octave
+             << " must be atoms-per-frame * bins-per-octave ("
+             << m_p.atomsPerFrame << " * " << m_binsPerOctave << " = "
+             << m_p.atomsPerFrame * m_binsPerOctave << ")" << endl;
+        throw std::invalid_argument
+            ("Column height must match atoms-per-frame * bins-per-octave");
+    }
+
     ComplexSequence transformed = m_kernel->processInverse(column);
 
     int halfLen = m_p.fftSize/2 + 1;
@@ -261,7 +320,50 @@
 
     m_fft->inverse(ri.data(), ii.data(), timeDomain.data());
 
-
-    //...
+    overlapAddAndResample(octave, timeDomain);
 }
 
+void
+CQInverse::overlapAddAndResample(int octave, const RealSequence &seq)
+{
+    // 4. Overlap-add each octave's resynthesised blocks (unwindowed)
+    //
+    // and
+    //
+    // 5. Resample each octave's overlap-add stream to the original
+    // rate
+
+    if (seq.size() != m_olaBufs[octave].size()) {
+        cerr << "ERROR: CQInverse::overlapAdd: input sequence length ("
+             << seq.size() << ") is expected to match OLA buffer size ("
+             << m_olaBufs[octave].size() << ")" << endl;
+        throw std::invalid_argument
+            ("Input sequence length should match OLA buffer size");
+    }
+
+    RealSequence toResample(m_olaBufs[octave].begin(),
+                            m_olaBufs[octave].begin() + m_p.fftHop);
+
+    RealSequence resampled = 
+        octave > 0 ?
+        m_upsamplers[octave]->process(toResample.data(), toResample.size()) :
+        toResample;
+
+    m_buffers[octave].insert(m_buffers[octave].end(),
+                             resampled.begin(),
+                             resampled.end());
+    
+    m_olaBufs[octave] = RealSequence(m_olaBufs[octave].begin() + m_p.fftHop,
+                                     m_olaBufs[octave].end());
+    
+    RealSequence pad(m_p.fftHop, 0);
+
+    m_olaBufs[octave].insert(m_olaBufs[octave].end(),
+                             pad.begin(),
+                             pad.end());
+
+    for (int i = 0; i < m_p.fftSize; ++i) {
+        m_olaBufs[octave][i] += seq[i];
+    }
+}
+
--- a/cpp-qm-dsp/CQInverse.h	Fri May 09 11:41:44 2014 +0100
+++ b/cpp-qm-dsp/CQInverse.h	Fri May 09 15:16:48 2014 +0100
@@ -44,7 +44,7 @@
     CQInverse(double sampleRate,
 	      double minFreq, double maxFreq,
 	      int binsPerOctave);
-    ~CQInverse();
+    virtual ~CQInverse();
 
     virtual double getSampleRate() const { return m_sampleRate; }
     virtual int getBinsPerOctave() const { return m_binsPerOctave; }
@@ -76,6 +76,7 @@
 
     std::vector<Resampler *> m_upsamplers;
     std::vector<RealSequence> m_buffers;
+    std::vector<RealSequence> m_olaBufs; // fixed-length, for overlap-add
     
     int m_outputLatency;
 
@@ -84,6 +85,8 @@
     void initialise();
     void processOctave(int octave, const ComplexBlock &block);
     void processOctaveColumn(int octave, const ComplexColumn &column);
+    void overlapAddAndResample(int octave, const RealSequence &);
+    RealSequence drawFromBuffers();
 };
 
 #endif
--- a/cpp-qm-dsp/CQSpectrogram.h	Fri May 09 11:41:44 2014 +0100
+++ b/cpp-qm-dsp/CQSpectrogram.h	Fri May 09 15:16:48 2014 +0100
@@ -47,7 +47,7 @@
                   double minFreq, double maxFreq,
                   int binsPerOctave,
                   Interpolation interpolation);
-    ~CQSpectrogram();
+    virtual ~CQSpectrogram();
 
     virtual double getSampleRate() const { return m_cq.getSampleRate(); }
     virtual int getBinsPerOctave() const { return m_cq.getBinsPerOctave(); }
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpp-qm-dsp/processfile.cpp	Fri May 09 15:16:48 2014 +0100
@@ -0,0 +1,119 @@
+
+#include "ConstantQ.h"
+#include "CQInverse.h"
+
+#include <sndfile.h>
+
+#include <iostream>
+
+using std::vector;
+using std::cerr;
+using std::endl;
+
+#include <cstring>
+
+int main(int argc, char **argv)
+{
+    if (argc != 3) {
+	cerr << "Usage: " << argv[0] << " infile.wav outfile.wav" << endl;
+	return 2;
+    }
+
+    char *fileName = strdup(argv[1]);
+    char *fileNameOut = strdup(argv[2]);
+
+    SNDFILE *sndfile;
+    SNDFILE *sndfileOut;
+    SF_INFO sfinfo;
+    SF_INFO sfinfoOut;
+    memset(&sfinfo, 0, sizeof(SF_INFO));
+
+    sndfile = sf_open(fileName, SFM_READ, &sfinfo);
+    if (!sndfile) {
+	cerr << "ERROR: Failed to open input file \"" << fileName << "\": "
+	     << sf_strerror(sndfile) << endl;
+	return 1;
+    }
+
+    sfinfoOut.channels = 1;
+    sfinfoOut.format = sfinfo.format;
+    sfinfoOut.frames = sfinfo.frames;
+    sfinfoOut.samplerate = sfinfo.samplerate;
+    sfinfoOut.sections = sfinfo.sections;
+    sfinfoOut.seekable = sfinfo.seekable;
+
+    sndfileOut = sf_open(fileNameOut, SFM_WRITE, &sfinfoOut) ;
+    if (!sndfileOut) {
+	cerr << "ERROR: Failed to open output file \"" << fileNameOut << "\" for writing: "
+	     << sf_strerror(sndfileOut) << endl;
+	return 1;
+    }
+    
+    int ibs = 1024;
+    int channels = sfinfo.channels;
+    float *fbuf = new float[channels * ibs];
+
+    ConstantQ cq(sfinfo.samplerate,
+		 100, sfinfo.samplerate / 3,
+		 60);
+
+    CQInverse cqi(sfinfo.samplerate,
+		 100, sfinfo.samplerate / 3,
+		 60);
+
+    int inframe = 0;
+    int outframe = 0;
+    int latency = cq.getLatency() + cqi.getLatency();
+
+    while (inframe < sfinfo.frames) {
+
+        int count = -1;
+	
+	if ((count = sf_readf_float(sndfile, fbuf, ibs)) < 0) {
+	    break;
+	}
+
+	vector<double> cqin;
+	for (int i = 0; i < count; ++i) {
+	    double v = fbuf[i * channels];
+	    if (channels > 1) {
+		for (int c = 1; c < channels; ++c) {
+		    v += fbuf[i * channels + c];
+		}
+		v /= channels;
+	    }
+	    cqin.push_back(v);
+	}
+	
+	vector<double> cqout = cqi.process(cq.process(cqin));
+
+	if (outframe >= latency) {
+
+	    sf_writef_double(sndfileOut, 
+			     cqout.data(), 
+			     cqout.size());
+
+	} else if (outframe + (int)cqout.size() >= latency) {
+
+	    int offset = latency - outframe;
+	    sf_writef_double(sndfileOut, 
+			     cqout.data() + offset,
+			     cqout.size() - offset);
+	}
+
+	inframe += count;
+	outframe += cqout.size();
+    }
+
+    vector<double> r1 = cqi.process(cq.getRemainingOutput());
+    vector<double> r2 = cqi.getRemainingOutput();
+
+    sf_writef_double(sndfileOut, r1.data(), r1.size());
+    sf_writef_double(sndfileOut, r2.data(), r2.size());
+
+    sf_close(sndfile);
+    sf_close(sndfileOut);
+
+    return 0;
+}
+
--- a/cpp-qm-dsp/test.cpp	Fri May 09 11:41:44 2014 +0100
+++ b/cpp-qm-dsp/test.cpp	Fri May 09 15:16:48 2014 +0100
@@ -50,7 +50,7 @@
 	in.push_back(sin(i * M_PI / 2.0));
     }
 
-    CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::None);
+    CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::InterpolateZeros);
 
     vector<vector<double> > out = k.process(in);
     vector<vector<double> > rest = k.getRemainingOutput();
--- a/vamp/CQVamp.cpp	Fri May 09 11:41:44 2014 +0100
+++ b/vamp/CQVamp.cpp	Fri May 09 15:16:48 2014 +0100
@@ -272,7 +272,7 @@
 
     if (m_cq) {
         char name[20];
-        for (int i = 0; i < d.binCount; ++i) {
+        for (int i = 0; i < (int)d.binCount; ++i) {
             float freq = m_cq->getBinFrequency(i);
             sprintf(name, "%.1f Hz", freq);
             d.binNames.push_back(name);