working source extraction

fdNMF.cpp

@@ -3,20 +3,43 @@
 
 #include "SC_PlugIn.h"
 #include <vector>
+#include <algorithm>
 #include "stft.h"
+#include "nmf.h"
+#include <Eigen/Dense>
+
+
+//Using statements for eigenmf. These will change
+using stft::STFT;
+using stft::ISTFT;
+using stft::Spectrogram;
+using stft:: audio_buffer_t;
+using stft:: magnitude_t;
+using nmf::NMF;
+using nmf::NMFModel;
+using Eigen::MatrixXcd;
+using Eigen::MatrixXd;
+using std::complex;
+using util::stlVecVec2Eigen;
+using util::Eigen2StlVecVec;
+using std::numeric_limits;
 
 static InterfaceTable *ft;
 
+//namespace to gaffatape
+namespace gaffatape {
+
 void BufNMF(World *world, struct SndBuf *dstBuf, struct sc_msg_iter *msg)
 {
 	int dstFrameCount = dstBuf->frames;
 	int dstChanCount = dstBuf->channels;
 
 	uint32 srcBufNum = msg->geti();
-	uint32 rank = msg->geti();
-	uint32 fftSize = msg->geti();
-	uint32 windowSize = msg->geti();
-	uint32 hopSize = msg->geti();
+	long rank = msg->geti();
+	long iterations = msg->geti();
+	long fftSize = msg->geti();
+	long windowSize = msg->geti();
+	long hopSize = msg->geti();
 
 	if (srcBufNum >= world->mNumSndBufs){
 		Print("fdNMF is not happy because the source buffer does not exist.\n");
@@ -33,32 +56,56 @@ void BufNMF(World *world, struct SndBuf *dstBuf, struct sc_msg_iter *msg)
 	int srcFrameCount = srcBuf->frames;
 	int srcChanCount = srcBuf->channels;
 
-	// if (dstChanCount < rank) {
-	// 	Print("fdNMF is not happy because the destination buffer has a lower channel count than the number of ranks.\n");
-	// 	return;
-	// }
+	if (dstChanCount < rank) {
+		Print("fdNMF is not happy because the destination buffer has a lower channel count than the number of ranks.\n");
+		return;
+	}
+
+	// check size of dstBuff
+	if (dstFrameCount < srcFrameCount) {
+		Print("fdNMF is not happy because the destination buffer shorter than the source buffer.\n");
+		return;
+	}
 
 	// make a vector of doubles for the samples
-	std::vector<double> tmp_vec(srcFrameCount);
+	std::vector<double> audio_in(srcFrameCount);
 
-	//construct STFT processors
-	stft::STFT stft(windowSize, fftSize, hopSize);
-	stft::ISTFT istft(windowSize, fftSize, hopSize);
+	//copied as is from max source (setting up the different variables and processes)
+	STFT stft(windowSize, fftSize, hopSize);
+	NMF nmfProcessor(rank, iterations);
+	ISTFT istft(windowSize, fftSize, hopSize);
 
 	//for each channels
-	for (int j=0;j<srcChanCount;j++){
+	// for (int j=0;j<srcChanCount;j++){
+	// just processing the first input channel instead of iterating through each channel, yet keeping the mechanism in there.
+	for (int j=0;j<1;j++){
 		//copies and casts to double the source samples
 		for (int i=0;i<srcFrameCount;i++){
-				tmp_vec[i] = srcBuf->data[(i*srcChanCount)+j];
+				audio_in[i] = srcBuf->data[(i*srcChanCount)+j];
 		}
 
-		// fft and ifft for fun
-		stft::Spectrogram spec = stft.process(tmp_vec);
-		std::vector<double>  out_vec = istft.process(spec);
+		Spectrogram spec = stft.process(audio_in);
+		magnitude_t mag = spec.magnitude();
+		NMFModel decomposition = nmfProcessor.process(mag);
+		MatrixXd W = stlVecVec2Eigen<double>(decomposition.W);
+		MatrixXd H = stlVecVec2Eigen<double>(decomposition.H);
 
-		//writes the output
-		for (int i=0;i<srcFrameCount;i++){
-				dstBuf->data[(i*srcChanCount)+j] = (float)out_vec[i];
+		MatrixXd V = W * H;
+
+		for (int i = 0; i < rank; i++)
+		{
+				MatrixXd source = W.col(i) * H.row(i);
+				MatrixXd filter = source.cwiseQuotient(V);
+				MatrixXcd specMatrix = stlVecVec2Eigen(spec.mData);
+				specMatrix = specMatrix.cwiseProduct(filter);
+				Spectrogram resultS(Eigen2StlVecVec<complex<double>>(specMatrix));
+
+				audio_buffer_t result = istft.process(resultS);
+
+				//writes the output
+				for (int k=0;k<srcFrameCount;k++){
+						dstBuf->data[(k*rank)+i] = (float)result[k];
+				}
 		}
 	}
 }
@@ -67,3 +114,4 @@ PluginLoad(OfflineFluidDecompositionUGens) {
 	ft = inTable;
 	DefineBufGen("BufNMF", BufNMF);
 }
+}

fdfNMF.sc

@@ -1,7 +1,7 @@
 // adds an instance method to the Buffer class
 + Buffer {
-	fdNMF { arg dstBuf, rank = 1, fftSize = 2048, windowSize = 2048, hopSize = 512;
+	fdNMF { arg dstBuf, rank = 1, iterations = 100, fftSize = 2048, windowSize = 2048, hopSize = 512;
 		if(bufnum.isNil) { Error("Cannot call % on a % that has been freed".format(thisMethod.name, this.class.name)).throw };
-		server.listSendMsg([\b_gen, dstBuf.bufnum, "BufNMF", bufnum, rank, fftSize, windowSize, hopSize])
+		server.listSendMsg([\b_gen, dstBuf.bufnum, "BufNMF", bufnum, rank, iterations, fftSize, windowSize, hopSize])
 	}
 }

tests.scd

@@ -1,13 +1,26 @@
 s.quit
 s.boot
 // allocates a 16-sample buffer and fills it with ascending values
-b = Buffer.alloc(s,50000);
-b.sine1([0,0,0,0,0,0,0,0,0,1],true,false);
+b = Buffer.read(s,"/Applications/Max.app/Contents/Resources/C74/media/msp/jongly.aif");
+b.play;
 b.plot;
-c = Buffer.alloc(s,50000);
+c = Buffer.alloc(s,b.numFrames,2);
 
 // run the code
-b.fdNMF(c, 5, 2048, 1024, 256);
+b.fdNMF(c,2)
 
 // read to check
-c.plot
+c.plot;
+{PlayBuf.ar(2,c.bufnum,doneAction:2)}.play
+
+
+// play a sinewave in the background to see if the server freaks out
+d = {SinOsc.ar(100,0,0.1)}.play
+e = Buffer.alloc(s,b.numFrames,5);
+
+// run the code
+b.fdNMF(e,5, 150, 2048,1024,256)
+
+// read to check
+e.plot;
+{Splay.ar(PlayBuf.ar(5,e.bufnum,doneAction:2))}.play