Checking inverse of convolution theorem

Question

Areeba Fatima 2022-1-27

2
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1637555-checking-inverse-of-convolution-theorem

编辑： Matt J 2022-1-28

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')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2022-1-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1637555-checking-inverse-of-convolution-theorem#answer_883200

编辑：Matt J 2022-1-27

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Areeba Fatima 2022-1-27

编辑：Areeba Fatima 2022-1-27

Thanks! A related question then: I see this arises from the periodicity for DFT, but is there any way to not use cconv as eventually I would want to move on to 2D signals.

Matt J 2022-1-27

编辑：Matt J 2022-1-28

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Answer 2

William Rose 2022-1-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1637555-checking-inverse-of-convolution-theorem#answer_883355

在 MATLAB Online 中打开

convolutionTheorem.m

@Areeba Fatima,

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