Differences

This shows you the differences between two versions of the page.

--- conv:conv [2008/01/24 12:08] – deva
+++ conv:conv [2008/09/04 12:09] (current) – deva
@@ Line 1: / Line 1: @@
-=====Convolution implementation pages (Spring 2008 studies group under Ole Caprani)=====
+=====Experiments with realtime convolution reverb - (Eksperimenter med realtids foldningsrumklang)=====
+Spring 2008 studies group under Ole Caprani\\
+\\
 Our goal is to implement a realtime plugin to do convolution.\\
 The basic "Input Side" is way too slow (convolution implemented with two for loops), so we will be experimenting with a dft-multiply-idft algorithm (sometimes called "The Overlap-and-Save algorithm") instead.\\
 To do the actual fft we will use the fftw3 library (http://www.fftw.org), which implements the "Fastest Fourier Transform in the West" algorithm.\\
-The next step is to implement at test the "Partitioned Overlap-and-Save" algorithm, where the filter is split into partitions prior to convolution.
+Currently we have implemented and tested convolution in the frequency domain where the input sound is partitioned which makes real time usage possible.
 ====Group====
@@ Line 28: / Line 30: @@
 As it is seen on the graphs the speed improves with the buffer size increasing, up to 1024 samples. At this point the iteration time starts fluctuating. This may be caused by internal cache size on the given CPU.
-====Links====
+====Literature====
-  * http://www.fftw.org - The fftw library.
+  * T. G. Stockham Jr. High-speed convolution and correlation. AFIPS Proc. 1966 Spring Joint Computer Conf., Vol 28, Spartan Books, 1966, pp. 229 - 233.
-  * http://archive.cnmat.berkeley.edu/~adrian/Csigproc.html - Some C optimizations that can be used when programming DSP.
+  * The Scientist and Engineer's Guide to Digital Signal Processing, copyright ©1997-1998 by Steven W. Smith. For more information visit the book's website at: [[http://www.DSPguide.com]]
-  * http://sunsite.univie.ac.at/Linux-soundapp/dsp.html - A list of links to various sites about audio and DSP programming.
+  * Real-time partitioned convolution for Ambiophonics surround sound. Torger, A.; Farina, A. Applications of Signal Processing to Audio and Acoustics, 2001 IEEE Workshop on the - Volume , Issue , 2001 Page(s):195 - 198
-  * http://www.dspguide.com - A very good online book about audio convolution.
+Further online references can be found in the links section.
-  * http://people.scs.fsu.edu/~burkardt/c_src/fftw3/fftw3.html - Some info about the fftw3 library.
-  * http://www.ludd.luth.se/~torger/brutefir.html#bruteconv_3 - Some info about partitioned convolution.
-  * http://www.aurora-plugins.it/Public/Papers/164-Mohonk2001.PDF - More on partitioned convolution.
-  * http://www.linuxdevcenter.com/lpt/a/586 - A LADSPA tutorial
 ====Code====
-Compile with:
+You can get the source code for the software used to create these tests here: {{:conv:conv-01242008.tar.gz}}\\
+The source code depends on [[http://www.mega-nerd.com/libsndfile/|libsndfile]] and [[http://www.fftw.org|fftw3]]. Compile it with:
 <code>
-g++ -O2 -mmmx -msse -msse2 conv_opt.cc -o conv_opt -lfftw3 -lsndfile
+make
 </code>
-or
+and run it with
 <code>
-g++ -O2 -mmmx -msse -msse2 conv_opt.cc -o conv_opt -DFFTW_FLOAT -lfftw3f -lsndfile
+./conv_opt [filterfile] [soundfile] [outputfile] [buffer size]
 </code>
-to use floating point fftw instead of double.\\
-The actual code is shown below:
+The input and filter files must be mono wav files.
-<code c>
-#include <stdlib.h>
-#include <fftw3.h>
-#include <sndfile.hh>
-#include <string.h>
-#include <assert.h>
-#include "sys/timeb.h"
-//#define BUFFER_SIZE_MS 113 //ms 113,  // Find a size with small prime factors
+====Examples====
-#define BUFFER_SIZE_MS 1000 //ms 113,  // Find a size with small prime factors
+Here are a dry input file called ''drums.wav'' and a wet output file called ''output.wav'': {{:conv:sounds.tar.gz}}\\
-#define SAMPLE_RATE_HZ 44100 //hz
+Get impulse responses here [[http://www.voxengo.com/impulses]].\\
-#define BUFFER_SIZE_SAMPLES 4096
+**NOTE**: Remember to convert them into mono before usage. Running the application on stereo files has undefined results.
-//(SAMPLE_RATE_HZ * BUFFER_SIZE_MS / 1000)
-//#define MEASURE
+====Realtime LADSPA plugin====
+We have experimented with implementing our code as a realtime ladspa plugin.\\
+This however is a very fragile implementation that only works on a hardcoded filter file and jack buffersize.\\
+If still interested, it can be fetched here: {{:conv:conv_opt_ladspa-0.1.tar.gz}}
-#ifdef FFTW_FLOAT
+====Links====
-#define FFTW_DATATYPE fftwf_complex
+  * http://www.fftw.org - The fftw library.
-#define FFTW_PLAN_DFT fftwf_plan_dft_1d
+  * http://archive.cnmat.berkeley.edu/~adrian/Csigproc.html - Some C optimizations that can be used when programming DSP.
-#define FFTW_EXECUTE fftwf_execute
+  * http://sunsite.univie.ac.at/Linux-soundapp/dsp.html - A list of links to various sites about audio and DSP programming.
-#define FFTW_PLAN fftwf_plan
+  * http://www.dspguide.com - A very good online book about audio convolution.
-#define FFTW_DESTROY_PLAN fftwf_destroy_plan
+  * http://people.scs.fsu.edu/~burkardt/c_src/fftw3/fftw3.html - Some info about the fftw3 library.
-#define FFTW_FREE fftwf_free
+  * http://www.ludd.luth.se/~torger/brutefir.html#bruteconv_3 - Some info about partitioned convolution.
-#define FFTW_MALLOC fftwf_malloc
+  * http://www.aurora-plugins.it/Public/Papers/164-Mohonk2001.PDF - More on partitioned convolution.
-typedef float sample_t;
+  * http://www.linuxdevcenter.com/lpt/a/586 - A LADSPA tutorial
-#else
-#define FFTW_DATATYPE fftw_complex
-#define FFTW_PLAN_DFT fftw_plan_dft_1d
-#define FFTW_EXECUTE fftw_execute
-#define FFTW_PLAN fftw_plan
-#define FFTW_DESTROY_PLAN fftw_destroy_plan
-#define FFTW_FREE fftw_free
-#define FFTW_MALLOC fftw_malloc
-typedef double sample_t;
-#endif
-typedef struct {
-	sample_t *data;
-	size_t size;
-} samples_t;
-samples_t *new_samples(sample_t *data, size_t size)
-{
-	samples_t *samples = (samples_t*)malloc(sizeof(samples_t));
-	samples->data = data;
-	samples->size = size;
-	return samples;
-}
-samples_t *read_data(const char *filename, double scale)
-{
-	printf("Loading %s ... ", filename); fflush(stdout);
-	SndfileHandle sfh(filename);
-	size_t size = sfh.seek(0, SEEK_END);
-	sfh.seek(0, SEEK_SET);
-	//	size_t size = 3;
-	sample_t *data = (sample_t*)malloc(sizeof(sample_t) * size);
-	sfh.read(data, size);
-	printf("[%d samples] ... ", size); fflush(stdout);
-	for(int i = 0; i < size; i++) {
-		data[i] = scale * data[i];
-	}
-  //	for(int i = size / 2; i < (size / 2) +20; i++) printf("\t[%d] %f\n", i, data[i]);
-	samples_t *samples = new_samples(data, size);
-	printf("done\n");
-	return samples;
-}
-void write_data(const char *filename, samples_t *samples)
-{
-	printf("Saving %s ...", filename); fflush(stdout);
-	//	for(int i = samples->size / 2; i < (samples->size / 2) +20; i++)printf("\t[%d] %f\n", i, samples->data[i]);
-	SndfileHandle sfh(filename, SFM_WRITE, SF_FORMAT_WAV | SF_FORMAT_PCM_16, 1, 44100);
-	sfh.write(samples->data, samples->size);
-	printf("done\n");
-}
-#define RE 0
-#define IM 1
-static FFTW_DATATYPE *hfft = NULL;
-static FFTW_DATATYPE *outm = NULL;
-static FFTW_DATATYPE *out = NULL;
-static FFTW_DATATYPE *inx = NULL;
-static FFTW_DATATYPE *xfft = NULL;
-void init_fftw(samples_t *h)
-{
-  FFTW_PLAN p;
-	int size = BUFFER_SIZE_SAMPLES + h->size - 1;
-	hfft = (FFTW_DATATYPE*) FFTW_MALLOC(sizeof(FFTW_DATATYPE) * size);
-  FFTW_DATATYPE *inh = (FFTW_DATATYPE*) FFTW_MALLOC(sizeof(FFTW_DATATYPE) * size);
-  for(int i = 0; i < size; i++) {
-		inh[i][RE] = inh[i][IM] = 0;
-   }
-  for(int i = 0; i < h->size; i++) {
-    inh[i][RE] = h->data[i];
-  }
-	p = FFTW_PLAN_DFT(size, inh, hfft, FFTW_FORWARD, FFTW_ESTIMATE);
-  FFTW_EXECUTE(p);
-	FFTW_DESTROY_PLAN(p);
-	FFTW_FREE(inh);
-	outm = (FFTW_DATATYPE*) FFTW_MALLOC(sizeof(FFTW_DATATYPE) * size);
-  out = (FFTW_DATATYPE*) FFTW_MALLOC(sizeof(FFTW_DATATYPE) * size);
-  inx = (FFTW_DATATYPE*) FFTW_MALLOC(sizeof(FFTW_DATATYPE) * size);
-  for(int i = 0; i < size; i++) {
-    inx[i][RE] = inx[i][IM] = 0;
-  }
-  xfft = (FFTW_DATATYPE*) FFTW_MALLOC(sizeof(FFTW_DATATYPE) * size);
-}
-void deinit_fftw()
-{
-	FFTW_FREE(hfft);
-  FFTW_FREE(out);
-  FFTW_FREE(outm);
-  FFTW_FREE(inx);
-  FFTW_FREE(xfft);
-}
-#ifdef MEASURE
-long static time_used;
-#endif
-samples_t *conv_fftwc_preprocess(samples_t *x, size_t hsize/*samples_t *h*/)
-{
-	assert(hfft); // init_fftw not called!
-	assert(out); // init_fftw not called!
-	assert(outm); // init_fftw not called!
-	assert(inx); // init_fftw not called!
-	assert(xfft); // init_fftw not called!
-#ifdef MEASURE
-  struct timeb systime;
-  long begin = 0;
-  long end = 0;
-  ftime(&systime);
-  begin = systime.time * 1000 + systime.millitm;
-#endif
-  FFTW_PLAN p;
-  int size = x->size + hsize - 1;
-  sample_t *data = (sample_t*)calloc(sizeof(sample_t), size);
-  samples_t *y = new_samples(data, size);
-	for(int i = 0; i < x->size; i++) {
-    inx[i][RE] =  x->data[i];
-  }
-  p = FFTW_PLAN_DFT(size, inx, xfft, FFTW_FORWARD, FFTW_ESTIMATE);
-  FFTW_EXECUTE(p);
-  for(int i = 0; i < size; i++ ) {
-    outm[i][RE] = hfft[i][RE] * xfft[i][RE] - hfft[i][IM] * xfft[i][IM];
-    outm[i][IM] = xfft[i][IM] * hfft[i][RE] + xfft[i][RE] * hfft[i][IM];
-  }
-	p = FFTW_PLAN_DFT(size, outm, out, FFTW_BACKWARD, FFTW_ESTIMATE);
-  FFTW_EXECUTE(p);
-  for(int i = 0; i < size; i++) {
-    data[i] = out[i][RE] / size;
-  }
-  FFTW_DESTROY_PLAN(p);
-#ifdef MEASURE
-  ftime(&systime);
-  end = systime.time * 1000 + systime.millitm;
-  printf("\t Time consumed: %d\n", end - begin);
-  time_used = end - begin;
-#endif
-	return y;
-}
-sample_t *filter(const sample_t *input, size_t filter_size)
-{
-	unsigned int buffer_size =  BUFFER_SIZE_SAMPLES + filter_size - 1;
-	static sample_t *buffer = (sample_t*)calloc(buffer_size, sizeof(sample_t));
-	static unsigned int p = 0;
-	samples_t *input_s = new_samples((sample_t *)input, BUFFER_SIZE_SAMPLES);
-	samples_t *y;
-	y = conv_fftwc_preprocess(input_s, filter_size);
-	for(int i = 0; i < BUFFER_SIZE_SAMPLES; i++) {
-		buffer[(i + p + filter_size - 1) % buffer_size] = 0; //y->data[i];
-	}
-	for(int i = 0/* BUFFER_SIZE_SAMPLES*/; i < buffer_size; i++) {
-		buffer[(i + p) % buffer_size] += y->data[i];
-	}
-	// alloker returbuffer
-	sample_t *output_buffer = (sample_t*)malloc(BUFFER_SIZE_SAMPLES * sizeof(sample_t));
-	// kopier retsultat til output buffer
-	for(int i = 0; i < BUFFER_SIZE_SAMPLES; i++) {
-		output_buffer[i] = buffer[(i + p) % buffer_size];
-	}
-	p += BUFFER_SIZE_SAMPLES;
-	// frigør y
-	free(input_s);
-	return output_buffer;
-}
-sample_t puls = 0.5;
-int main(int argc, char *argv[])
-{
-	//	y(n)=\sum_{k=-\infty}^\infty h(k)x(n-k)
-	if(argc != 4) {
-		printf("Usage %s [filterfile] [soundfile] [outputfile]\n", argv[0]);
-		exit(1);
-	}
-	int xsz, hsz;
-	const char *filterfile = argv[1];
-	const char *inputfile = argv[2];
-	const char *outputfile = argv[3];
-	samples_t *h = read_data(filterfile, 0.08);
-	//	samples_t *h = new_samples(&puls, 0.5);
-	samples_t *x = read_data(inputfile, 0.5);
-	//	samples_t *x = new_samples(&puls, 1);
-	size_t s = h->size + x->size - 1;
-	sample_t *d = (sample_t*)malloc(s  * sizeof(sample_t));
-	samples_t *y = new_samples(d, s);
-  // 0 padding
-  sample_t *pd = (sample_t*)malloc((y->size) * sizeof(sample_t));
-  for(int i = 0; i < x->size; i++) {
-    pd[i] = x->data[i];
-  }
-  for(int i = x->size; i < y->size + (y->size % BUFFER_SIZE_SAMPLES); i++) {
-    pd[i] = 0;
-  }
-  x = new_samples(pd, x->size);
-	init_fftw(h);
-  int e = 0;
-#ifdef MEASURE
-  // find a good delta value (ie. small prime factors)
-  int o = 20000;
-  for(int i = 0; i < 20; i++) {
-    printf("Measuring %d: ", i);
-    filter(&(x->data[0]), h->size - i);
-    if(time_used < o)  { e = i; o = time_used; }
-  }
-#endif
-  for(int i = 0; i < y->size; i += BUFFER_SIZE_SAMPLES) {
-#ifdef MEASURE
-    printf("[c:%d, s:%d](%d, %d)", i, y->size, BUFFER_SIZE_SAMPLES, e); fflush(stdout);
-#endif
-    sample_t *data = filter(&(x->data[i]), h->size /*- e*/);
-		for(int j = 0; j < BUFFER_SIZE_SAMPLES; j++) {
-			y->data[j + i] = data[j];
-		}
-    free(data);
-	}
-	deinit_fftw();
-	write_data(outputfile, y);
-	free(y);
-	free(h);
-	free(x);
-}
-</code>