Show Equivalence Between Multiplication in Time Domain to Convolution in Frequency Domain

Question

My goal is to compute the Fourier of the product between two discrete time signals, y1 and y2. This can be done by computing the convolution between the fourier transform of y1, f1 and the fourier transform of y2, f2 (This is what I understood from the Wikipedia Page).

I tried the following code in MATLAB, it shows that the convolution between y1 and y2 is equal to ifft(f1.*f2) as shown on the Wikipedia page, the code worked perfectly:

y1=sin(x);
y2=sin(3*x);
n=length(y1);
convol=conv(y1,y2);
f1=fft(y1,n*2-1);
f2=fft(y2,n*2-1);
plot(ifft(f1.*f2))
hold on
plot(convol,'x')

Now I want to be able to compute the fourier transform of the product of y1 by y2 using the convolution theorem, this can be done by computing the convolution between f1 and f2, so I wrote the following code:

x=0:0.01:10;
y1=sin(x);
y2=sin(3*x);
n=length(y1);
f1=fft(y1,n*2-1);
f2=fft(y2,n*2-1);
convol=conv(f1,f2);
ffty1y2=fft(y1.*y2,length(convol));
plot(abs(ffty1y2))
hold on
plot(abs(convol)/length(convol),'x')

As you can see the results do not match (neither in location of the peaks nor in their values), what could be wrong here? How do I fix this?

Answer 1

Let's assume you have 2 signals: vX and vY.
So:

clear();

numSamplesX = length(vX);
numSamplesY = length(vY);

numSamplesConv = numSamplesX + numSamplesY - 1;

vTimeDomainConv = conv(vX, vY);

vFrequencyDomainConv = ifft(fft(vX, numSamplesConv) .* fft(vY, numSamplesConv), 'symmetric');

max(abs(vTimeDomainConv - vFrequencyDomainConv)) %<! Should be < 1e-12

This is the Convolution Theorem for Discrete Signals to show convolution in time domain is equivalent to element wise multiplication in frequency domain.

If you want to show element wise multiplication in time domain can be done using the convolution in frequency domain you need to either interpolate the time domain signal to length of linear convolution or just use circular convolution in frequency domain:

clear();

numSamples = 6;

vX = randn(numSamples, 1);
vY = randn(numSamples, 1);

vXY = vX .* vY; %<! Time Domain Mulaiplication
vDxy = fft(vXY); %<! The DFT of the Time Domain Multiplication

vDx = fft(vX);
vDy = fft(vY);

% Applying circular convolution
vXYC = cconv(vDx, vDy, numSamples) / numSamples; %<! Normalization

% Comparing the result in Frequency Domain
max(abs(vXYC - vDxy)) %<! Should be < 1e-12

% Time Domain
vXYFromFrequency = ifft(vXYC, 'symmetric');

% Comparing the result in Time Domain
max(abs(vXY - vXYFromFrequency)) %<! Should be < 1e-12

Answer 2

The problem with your code is that you need to get the indexing and time gating right. Matlab represents the frequencies in the FFT vector somewhat out of order , i.e. from $[0,N-1]$ instead of $[-N/2+1,N/2]. Once you rotate them so that DC is in the center and grab the correct frequencies after convolution, this does indeed work.

%% Create two signals
n = 128;
x1 = cos(2*pi*(0:n-1)'/n*6);
x2 = cos(2*pi*(0:n-1)'/n*17);
% reference in the time domain
yReference = x1.*x2;

%% frequency domain convolution
fx1 = fft(x1);
fx2 = fft(x2);
% rotate so DC is in the center, convolve and scale
fy = 1/n.*conv(circshift(fx1,n/2),circshift(fx2,n/2));
% cut out the middle part and rotate back to DC in the front
fyCut = circshift(fy(n/4+(1:n)),n/4);
% back to the time domain
y = real(ifft(fyCut));

%% calculate and report error
err = sum((y-yReference).^2)./sum(yReference.^2);
fprintf('Error = %6.2f dB\n',10*log10(err));