diff --git a/ex-6/fft.c b/ex-6/fft.c index 8a5de4b..59a0877 100644 --- a/ex-6/fft.c +++ b/ex-6/fft.c @@ -115,11 +115,32 @@ gsl_histogram* fft_deconvolve( * 3. copy the kernel to the view */ gsl_vector *vpadded = gsl_vector_calloc(data->n); - gsl_vector center = gsl_vector_subvector( - vpadded, // vector: padded kernel (n+m-1) - orig_size/2, // offset: half of original data size (m/2) - kernel->n).vector; // length: kernel size (n) - gsl_vector_memcpy(¢er, &vkernel); + + /* Copy the first half (origin + positive values) + * into the leftmost side of the padded kernel vector. + */ + gsl_vector pad_left = gsl_vector_subvector( + vpadded, // source + 0, // offset + kernel->n/2 + 1).vector; // size + gsl_vector ker_left = gsl_vector_subvector( + &vkernel, + kernel->n/2, + kernel->n/2 + 1).vector; + gsl_vector_memcpy(&pad_left, &ker_left); + + /* Copy the second half (negative values) + * into the rightmost side of the padded kernel vector. + */ + gsl_vector pad_right = gsl_vector_subvector( + vpadded, // source + data->n - kernel->n/2, // offset + (kernel->n/2)).vector; // size + gsl_vector ker_right = gsl_vector_subvector( + &vkernel, + 0, + kernel->n/2).vector; + gsl_vector_memcpy(&pad_right, &ker_right); /* Compute the DFT of the data and * divide it by the DFT of the kernel. @@ -147,20 +168,6 @@ gsl_histogram* fft_deconvolve( htable, // wavetable (complex) wspace); // workspace - /* Restore the natural order of frequency - * in the result of the DFT (negative→positive). - * To do that roll over the array by a half: - * 1. create a temp array of half size - * 2. copy first half over to temp - * 3. copy second half to the first - * 4. copy back the first to the second half - */ - size_t half_size = data->n/2 * sizeof(double); - double *temp = malloc(half_size); - memcpy(temp, res, half_size); - memcpy(res, res + data->n/2, half_size); - memcpy(res + data->n/2, temp, half_size); - /* Create a histogram with the same edges * as `data`, but with the original size, * to return the cleaned result @@ -190,7 +197,6 @@ gsl_histogram* fft_deconvolve( gsl_fft_halfcomplex_wavetable_free(htable); gsl_fft_real_workspace_free(wspace); free(res); - free(temp); return hist; } diff --git a/ex-6/main.c b/ex-6/main.c index a204cc8..d3b9831 100644 --- a/ex-6/main.c +++ b/ex-6/main.c @@ -68,7 +68,7 @@ gsl_histogram* gaussian_for(gsl_histogram *hist, double sigma) { * gaussian in the middle. */ size_t n = ((double) hist->n * 6) / 100; - n = n % 2 ? n+1 : n; + n = n % 2 ? n : n+1; /* Calculate the (single) bin width assuming * the ranges are uniformely spaced @@ -82,11 +82,11 @@ gsl_histogram* gaussian_for(gsl_histogram *hist, double sigma) { /* The histogram will be such that * the maximum falls in the central bin. */ - long int offset = res->n/2; + long int offset = (res->n)/2; for (long int i = 0; i < (long int)res->n; i++) res->bin[i] = gsl_ran_gaussian_pdf( - ((double)(i - offset) + 0.5) * dx, sigma * dx); + ((double)(i - offset)) * dx, sigma * dx); return res; }