Giu/analistica
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

ex-6: fix fft deconvolution

Browse Source 
The Fourier transform of the kernel wasn't implemented
correctly and resulted in 0 division when the bins are
odd-numbered, moreover the whole histogram was shifted
by one bin to the right.
This commit is contained in:
Michele Guerini Rocco 2020-05-05 22:26:20 +00:00
parent c6939be6d0
commit 67c9d46188
2 changed files with 29 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
ex-6/fft.c
View File

@ -115,11 +115,32 @@ gsl_histogram* fft_deconvolve(
* 3. copy the kernel to the view * 3. copy the kernel to the view
*/ */
gsl_vector *vpadded = gsl_vector_calloc(data->n); gsl_vector *vpadded = gsl_vector_calloc(data->n);
gsl_vector center = gsl_vector_subvector(
vpadded, // vector: padded kernel (n+m-1) /* Copy the first half (origin + positive values)
orig_size/2, // offset: half of original data size (m/2) * into the leftmost side of the padded kernel vector.
kernel->n).vector; // length: kernel size (n) */
gsl_vector_memcpy(&center, &vkernel); 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 /* Compute the DFT of the data and
* divide it by the DFT of the kernel. * divide it by the DFT of the kernel.
@ -147,20 +168,6 @@ gsl_histogram* fft_deconvolve(
htable, // wavetable (complex) htable, // wavetable (complex)
wspace); // workspace 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 /* Create a histogram with the same edges
* as `data`, but with the original size, * as `data`, but with the original size,
* to return the cleaned result * to return the cleaned result
@ -190,7 +197,6 @@ gsl_histogram* fft_deconvolve(
gsl_fft_halfcomplex_wavetable_free(htable); gsl_fft_halfcomplex_wavetable_free(htable);
gsl_fft_real_workspace_free(wspace); gsl_fft_real_workspace_free(wspace);
free(res); free(res);
free(temp);
return hist; return hist;
} }

6
ex-6/main.c
View File

@ -68,7 +68,7 @@ gsl_histogram* gaussian_for(gsl_histogram *hist, double sigma) {
* gaussian in the middle. * gaussian in the middle.
*/ */
size_t n = ((double) hist->n * 6) / 100; 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 /* Calculate the (single) bin width assuming
* the ranges are uniformely spaced * 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 histogram will be such that
* the maximum falls in the central bin. * 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++) for (long int i = 0; i < (long int)res->n; i++)
res->bin[i] = gsl_ran_gaussian_pdf( res->bin[i] = gsl_ran_gaussian_pdf(
((double)(i - offset) + 0.5) * dx, sigma * dx); ((double)(i - offset)) * dx, sigma * dx);
return res; return res;
} }