Checking inverse of convolution theorem

45 views (last 30 days)
Areeba Fatima
Areeba Fatima on 27 Jan 2022
Edited: Matt J on 28 Jan 2022
One of the forms of the convolution theorem says that if a function can be expressed as the product of two functions say f(x) and g(x) then the Fourier transform of it is given by the convolution of the Fourier transform of the two functions i.e.,
FT{f(x).g(x)} = FT{f(x)}*FT{g(x)}
Here, . represents elementwise multiplication, and * represents convolution.
I am trying to check this on MATLAB, but the results for the left handside of the equation and the right hand side of the equation don't match perfectly. I have used the other form of convolution many times (i.e. h(x)*b(x) = IFT[FT{h(x)}.FT{g(x)}] ) and the results match perfectly for well zero padded cases. So, don't understand what am I missing here. I am mostly interested in checking the phase values from both the sides. My code is:
close all
clear all
clc
A = padarray(rand(15,1),10,0,'both');
B = padarray(rand(15,1),10,0,'both');
Prod = A.*B;
Res1 = fftshift(fft2(ifftshift(Prod))); % LHS of the equation
Value1 = fftshift(fft2(ifftshift(A)));
Value2 = fftshift(fft2(ifftshift(B)));
%
Value1_pad = padarray(Value1,25,0,'both');
Value2_pad = padarray(Value2,25,0,'both');
%
Res2 = conv(Value1,Value2,'same'); %RHS of the equation
figure(1)
plot(real(Res1))
hold on
plot((real(Res2))/(15+20));
hold off
title('real plot')
legend('FT of product','conv of FTs')
figure(2)
plot(angle(Res1));
hold on
plot(angle(Res2));
hold off
title('phase plot')
legend('FT of product','conv of FTs')

Accepted Answer

Matt J
Matt J on 27 Jan 2022
Edited: Matt J on 27 Jan 2022
You need to use circulant convolution, rather than linear convolution. Also, your fftshift must be applied to the output of the convolution, not to the inputs.
A = padarray(rand(15,1),10,0,'both');
B = padarray(rand(15,1),10,0,'both');
Prod = A.*B;
Res1 = fftshift( fft2(ifftshift(Prod)) ); % LHS of the equation
Value1 = fft2(ifftshift(A));
Value2 = fft2(ifftshift(B)) ;
N=numel(Value1);
Res2 = fftshift( cconv(Value1,Value2,N)/N ); %RHS of the equation
figure(1)
plot(real(Res1))
hold on
plot((real(Res2)),'x');
hold off
title('real plot')
legend('FT of product','conv of FTs')
figure(2)
plot(angle(Res1));
hold on
plot(angle(Res2),'x');
hold off
title('phase plot')
legend('FT of product','conv of FTs')
  2 Comments
Matt J
Matt J on 27 Jan 2022
Why would you want to use the right hand size fo the equation? The left hand side is usually the more computationally efficient.
In any case, you can implement 2D circulant convolution as follows:
conv2( x( [2:end,1:end], [2:end,1:end] ), y,'valid' );
wen x and y are equal-sized matrices.

Sign in to comment.

More Answers (1)

William Rose
William Rose on 27 Jan 2022
I realize that you have already accepted the excellent answer from @Matt J. Here is a script that demonstrates the equivalence of the two forms of the convoution theorem. This example does not use padding and does not use shifting. Thus it is a simpler demonstartion of the equivalence.
As @Matt J said, the key is to do a circular convolution. The circular convolution is divided by N, to normalize correctly.
I use lower case for time domain and upper case for frequency domain.
The bottom left plot shows time domain: x3=x1.*x2 and x3 computed by applying the convolution theorem in the frequency domian. The two series are identical, showing a match between the two ways of computing x3.
The bottom right plot shows frequency domain: X3=fft(x1.*x2) and the circular convolution of X1, X2. The two series are identical, showing a match between the two ways of computing X3.
The key lines of code are below; see that attached file for the full script.
x1=sin(2*pi*10*t);
X1=fft(x1); %Fourier transform of x1(t)
x2=0.5*(1-cos(2*pi*t));
X2=fft(x2); %Fourier transform of x2(t)
x3=x1.*x2; %x3=point-by-point product
X3=fft(x3); %Fourier transform of x1.*x2
X3cc=cconv(X1,X2,N)/N; %circular convolution of X1,X2
x3ifft=ifft(X3); %inv.FT of cconv(X1,X2)/N

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by