Aasimon.org Wiki

• http://archive.cnmat.berkeley.edu/~adrian/Csigproc.html
• http://sunsite.univie.ac.at/Linux-soundapp/dsp.html
• http://www.dspguide.com/ch9/3.htm
• http://people.scs.fsu.edu/~burkardt/c_src/fftw3/fftw3.html
• http://www.ludd.luth.se/~torger/brutefir.html#bruteconv_3

Compile with:

g++ conv.cc -o conv -lfftw

#include <stdlib.h>
#include <fftw3.h>
 
//typedef short int sample_t;
typedef double sample_t;
 
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)
{
	size_t size = 3;
	sample_t *data = (sample_t*)malloc(sizeof(sample_t) * size);
	data[0] = 1;
	data[1] = 1;
	data[2] = 1;
	samples_t *samples = new_samples(data, size);
	return samples;
}
 
 
void write_data(const char *filename, samples_t *samples)
{
	for(int i = 0; i < samples->size; i++) {
		printf("s[%d] = %f\n", i, samples->data[i]);
	}
}
 
#define RE 0
#define IM 1
 
samples_t *conv_fftw(samples_t *x, samples_t *h)
{
	int size = x->size + h->size - 1;
	sample_t *data = (sample_t*)calloc(sizeof(sample_t), size);
 
	samples_t *y = new_samples(data, size);
 
	// FFTW
	fftw_plan p;
 
	int N = x->size;
 
	fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
	fftw_complex *xfft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
	fftw_complex *hfft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
 
	p = fftw_plan_dft_r2c_1d(x->size, x->data, xfft, 0);
	fftw_execute(p);
 
	/*
		for(int i = 0; i < N; i++) {
		printf("fftw[%d] = %f + i%f\n", i, xfft[i][0], xfft[i][1]);
		}
	*/
 
	p = fftw_plan_dft_r2c_1d(h->size, h->data, hfft, 0);
	fftw_execute(p);
 
	for(int i = 0; i < N; i++ ) {
		out[i][RE] = hfft[i][RE] * xfft[i][RE] - hfft[i][IM] * xfft[i][IM];
		out[i][IM] = xfft[i][IM] * hfft[i][RE] + xfft[i][RE] * hfft[i][IM];
	}
 
	// invers FFTW
 
	p = fftw_plan_dft_c2r_1d(N, out, data, 0);
	fftw_execute(p);
 
	for(int i = 0; i < N; i++ ) {
		data[i] /= N;
	}
 
	fftw_destroy_plan(p);
	fftw_free(xfft);
	fftw_free(hfft);
	// end of FFTW
 
	return y;
}
 
 
samples_t *conv_simple(samples_t *x, samples_t *h)
{
	int size = x->size + h->size - 1;
	sample_t *data = (sample_t*)calloc(sizeof(sample_t), size);
 
	samples_t *y = new_samples(data, size);
 
 
	for(int xn = 0; xn < x->size; xn++) {
		for(int hn = 0; hn < h->size; hn++) {
			y->data[xn + hn] += h->data[hn] * x->data[xn];
		}
	}
 
	return y;
}
 
 
int main(int argc, char *argv[])
{
	//	y(n)=\sum_{k=-\infty}^\infty h(k)x(n-k)
 
	int xsz, hsz;
 
	samples_t *h = read_data("h.dat");
	samples_t *x = read_data("x.dat");
 
	printf("\nInput Side Algorithm:\n");
	samples_t *y = conv_simple(x, h);
	write_data("y_simple.dat", y);
	free(y);
 
 
	printf("\nFFTW Algorithm:\n");
	y = conv_fftw(x, h);
	write_data("y_fftw.dat", y);
	free(y);
 
	free(h);
	free(x);
}