Convolution and Deconvolution in C

Question

I understand the theoretical foundations of convolution, but now that I'm trying to program it I'm having some issues conceptually.

Say I have two blocks of 64 audio samples each. I also have access to fft and inverse fft functions. I want to take these two blocks and convolve them. Is this just taking the Fourier transform of each block and multiplying the resulting elements pairwise?

Say now that I want to deconvolve these two (e.g., its a signal and its convolution). Would I just divide the elements pairwise?

For example, the output of my Fourier transform is this pair of amplitudes.

Real:
11.09  -0.2055  1.214  1.101  1.66  4.594  -3.809  -0.3608  
4.198  3.572  0.3175  -0.5848  -5.129  2.26  -5.641  -1.039  
0.2445  1.905  -6.607  -1.083  1.726  -6.772  2.04  5.563  
-6.907  -4.122  0.276  4.894  -0.8119  1.064  -3.989  -0.1623  
-1.712  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  

Imaginary:
0  -3.147  0.7265  -2.434  -1.125  3.208  -3.066  -1.637  
-4.558  2.981  2.04  4.246  1.041  -0.1337  -5.165  1.238  
-1.223  -1.29  4.51  -4.716  -2.676  2.557  1.385  -1.173  
-3.938  -5.14  2.763  -3.018  -1.727  7.329  4.278  2.936  
0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  

for a block size of 64. If I know that other signal has a Fourier transform given by imp_real[64] and imp_imag[64], how do I practically convolve and deconvolve the signals?

Answer 1

Convolution, using FFT, is much faster for very long sequences. Circular convolution is based on FFT and Matlab's fftfilt() uses FFT on convolution step. It is explained very well when it is faster on its documentation.

Example code for convolution:

L = length(x)+length(y)-1;
c = ifft(fft(x, L) .* fft(y, L));

---

Deconvolution can be calculated using FFT as well. However, element-wise division may result NaN outputs when both elements are equal to zero. You either can use isnan() function, or add a very small number on both elements to fix this issue.

Example code for deconvolution:

epsilon = 2e-17;
y = ifft((epsilon+fft(c))./(epsilon+fft(x, length(c))));

Here c is the convolution of x and y. Leftover elements can be cut as; y(1:length(c)-length(x)+1).

---

Results from Matlab:

a = 1:3; b = 4:7; c1 = conv(a,b)   %Matlab's convolution

c1 =

     4    13    28    34    32    21

d1 = deconv(c1,b)   %Matlab's deconvolution

d1 =

     1     2     3

d2 = ifft((eps + fft(c1)./(eps + fft(b, length(c1)))))

d2 =

    1.0000    2.0000    3.0000   -0.0000         0   -0.0000

First 3 elements will be a value where you can cut after length(c1)-length(b)+1. You can implent it in C language easily.

Answer 2

Your signals are in the time domain, to convole:

• Inverse FFT both sets of time domain data

• You now have two sets of frequency domain data

• Multiply the two sets of frequency domain data together (pair-wise)

• You now have one combined array of frequency data

• FFT the combined frequency data

• You now have the convolved time domain data

Answer 3

You can't multiply just the FFT amplitudes, or just the FFT real parts separately. Instead you need to use complex multiplication of the complex FFT result elements.

Also, unless you zero-pad sufficiently, you will end up with a circular convolution, not an equivalent to a linear convolution.