Fourier transform of circular convolution is not equal to multiplication of two Fourier transforms?

In Fourier theory, regardless discrete or continuous Fourier transform, we have an important property So, the Fourier transform of circular convolution of f and g should be equal to multiplication of the Fourier transforms of f and g. Although I have followed strictly the tutorial from MATLAB, I cannot justify this property. What did I do wrong? Could you please explain this to me?

N = 32;                                 % number of grid points
L = 1;                                  % problem domain
h = L/N;                                % size of 1 cell = x(2) - x(1)
x = 0 : h : L-h;                        % grid points
scale = 2*pi / L;                       % to scale integer wavenumbers to wavenumbers for domain [0, 1]
xi = fftshift(-N/2 : N/2-1) * scale;    % wave numbers
f = exp(cos(2 * pi * x));               % function 1
g = exp(sin(2 * pi * x));               % function 2
dfdx = ifft(1j * xi .* fft(f));         % derivative of f using fft
dgdx = ifft(1j * xi .* fft(g));         % derivative of g using fft
close all
figure(Position=[100, 100, 800, 300])
% The derivatives are perfectly correct, thus the definition of xi is OK.
subplot(1,3,1)
plot(x, dfdx, 'k-', x, - 2*pi* f.*sin(2*pi*x), '--', LineWidth=2), grid on, title("$df/dx$", Interpreter="latex")
subplot(1,3,2)
plot(x, dgdx, 'k-', x,   2*pi* g.*cos(2*pi*x), '--', LineWidth=2), grid on, title("$dg/dx$", Interpreter="latex")
% Compute circular convolution following MATLAB documentation:
% https://uk.mathworks.com/help/signal/ug/linear-and-circular-convolution.html
fpad = [f, zeros(1, length(g)-1)];
gpad = [g, zeros(1, length(f)-1)];
ccirc = ifft(fft(fpad) .* fft(gpad));
f_cconv_g = ccirc(1:N);
subplot(1,3,3)
plot(x, f_cconv_g, LineWidth=2), grid on, title("$f \ast g$", Interpreter="latex")

f_lconv_g = conv(f, g);
% Comparison between linear convolution and circular convolution (following
% MATLAB document) ---> The result is very good, so no question here!
norm(ccirc - f_lconv_g, 'inf')      % comparison is OK

ans = 1.7764e-14

% My question is here!
f_cconv_g_fft = fft(f_cconv_g);
fhat_time_ghat = fft(f) .* fft(g);
% Mathematically, we know that fft(f-convolute-g) should be equal to fft(f)
% .* fft(g). But this is not verified.
figure(Position=[500, 500, 800, 300])
subplot(1,2,1), 
plot(x, real(f_cconv_g_fft), 'k-', x, real(fhat_time_ghat), '--', LineWidth=2), grid on
title("$\mathrm{Re}[\widehat{f \ast g}]$", Interpreter="latex")
subplot(1,2,2)
plot(x, imag(f_cconv_g_fft), 'k-', x, imag(fhat_time_ghat), '--', LineWidth=2), grid on
title("$\mathrm{Im}[\widehat{f \ast g}]$", Interpreter="latex")

[John D'Errico]

You also posted this exact same question as an answer on another question. I deleted that answer, since it was not an answer to anything.

Answer 1

[Paul]

Hi Khiem,

Consider the following to see the difference between circular and linear convolution.

N = 32;                                 % number of grid points
L = 1;                                  % problem domain
h = L/N;                                % size of 1 cell = x(2) - x(1)
x = 0 : h : L-h;                        % grid points
f = exp(cos(2 * pi * x));               % function 1
g = exp(sin(2 * pi * x));               % function 2

Zero padding the fft input is used to compute the linear convolution, over its duration, of two finite duration signals.

The duration of the linear convolution of f and g is

D = numel(f) + numel(g) - 1;

% this code works
% fpad = [f, zeros(1, length(g)-1)];
% gpad = [g, zeros(1, length(f)-1)];
% but I prefer
fpad = [f, zeros(1, D - numel(f))];
gpad = [g, zeros(1, D - numel(g))];

Now, because of the padding, the ifft computed next is really D elements from n = 0 to n = D-1 of the linear convolution of f and g, as must be the case based on comparison to the conv(f, g) that comes next.

%ccirc = ifft(fft(fpad) .* fft(gpad));
clinear = ifft(fft(fpad) .* fft(gpad));

%f_cconv_g = ccirc(1:N);
f_lconv_g = conv(f, g);
% Comparison between linear convolution computed in the time domain and the linear convolution
% computed via zero-padding in the frequency domain. Compare the entire
% sequence, not just the first N elements.
%norm(ccirc - f_lconv_g, 'inf')      % comparison is OK
norm(clinear - f_lconv_g, 'inf')      % comparison is OK

ans = 1.7764e-14

% My question is here!
%f_cconv_g_fft = fft(f_cconv_g);

Without zero padding, the product of the fft's of f and g is the circular convolution of f and g

fhat_time_ghat = fft(f) .* fft(g);

Or we can compute the circular convolution cconv of f and g, and then take the fft of the result

f_cconv_g_fft = fft(cconv(f,g,numel(f)));

Show equality

norm(f_cconv_g_fft - fhat_time_ghat,'inf')

ans = 4.1388e-14

figure(Position=[500, 500, 800, 300])
subplot(1,2,1)
plot(x, real(f_cconv_g_fft), 'k-', x, real(fhat_time_ghat), '--', LineWidth=2), grid on
title("$\mathrm{Re}[\widehat{f \ast g}]$", Interpreter="latex")
subplot(1,2,2)
plot(x, imag(f_cconv_g_fft), 'k-', x, imag(fhat_time_ghat), '--', LineWidth=2), grid on
title("$\mathrm{Im}[\widehat{f \ast g}]$", Interpreter="latex")

In summary, with f and g both signals of length N:

the product of the zero-padded fft's of f and g is the fft of the linear convolution of f and g (as long as the padding is long enough)

the product of the fft's of f and g is the ftt of the N-circular convolution of f and g

[Khiem Nguyen]

Thank you for your explanation and spot on code demonstration. I can also verify another property, namely by following your suggestion. I am working on a project that needs not only 1D convolution but also 2D and 3D circular convolution. The Signal Processing toolbox provides only 1D version . I would be very grateful if you can give me some hint on how to generalize 1D code to 2D and 3D version. I know we can do to compute the circular convolution . But, I want to do the the convolution of and , namely in the 2D and 3D. The algorithm I am working on is better established in the Fourier space. If I understand correctly, the code does not work in the Fourier space.

[Paul]

I'm not following the notation. Are f and g continous-time signals?

Does the asterisk (*) indicate a convolution integral, linear convolution sum, or circular convolution sum?

Does F[g] indicate the Continuous Time Fourier Transform of g(t) or the Discrete Time Fourier Transform of g[n] or something else?

[Khiem Nguyen]

I am sorry for the confusion. The asterisk denotes the circular convolution sum (as I am working with discrete values). F[f] denotes the fft(f), discrete Fourier transform of f.

Assume that we have two 2D arrays A and B of the same size. These arrays represent periodic functions (just like 1D example I provided in the code). Perhaps, there might be someway to perform “cconv2” just like “cconv” but for 2D array. What I want to ask is how can I generalize the cconv to the 2D arrays A and B so that I can still obtain an array C satisfying fft2(C) == cconv2(A, B, size(A)). Similarly, is there a way to perform “cconv” but for N-dimensional arrays so that fftn(C) = cconvn(A, B, size(A)).

[Paul]

Check out the Answer to this question by @Matt J that provides a function that generalizes cconv to the N-D case (even though the Question was only about the 2D case).

As an aside, why not just define a cconv2 function as

cconv2 = @(f,g) ifft2(fft2(g).*fft2(f));

[Khiem Nguyen]

Hi Paul: Sorry for my late response - I was super busy yesterday. What I wanted to do was to perform circular convolution in the Fourier space. Initially, I thought using

ifft(fft(fhat) .* fft(ghat))

is not possible because fhat and ghat are already in the Fourier space. However, it turns out possible as well -- the math proof is the same. However, as you pointed out @Matt J provided excellent solution using both approaches. I tried @Matt J's suggested code and it worked as a charm. Thank you for your super support.