conv:conv
This is an old revision of the document!
Table of Contents
Convolution implementation pages
Links
Code
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); }
conv/conv.1199959110.txt.gz · Last modified: 2008/01/10 10:58 by deva