2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Could you please tell me what is missing from my answer? The code does exactly what you wanted. $\endgroup$ – Royi Mar 28 at 11:46
3
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ My goal is to compute the Fourier transform between the product of y1 and y2, I want to do this using the convolution between the FFT of y1 and the FFT of y2, this is the second code in my question and as you can see it's not working. $\endgroup$ – M.H. Mar 27 at 8:37
  • $\begingroup$ What do you mean? You want to show you can calculate the element wise product in time domain using convolution in the spectrum? $\endgroup$ – Royi Mar 27 at 10:38
  • $\begingroup$ I want to compute the FFT of the element wise product using the convolution in the spectrum, I want to apply this for another problem but I first to need to check that it works but it's not working in the code I wrote. $\endgroup$ – M.H. Mar 27 at 12:54
  • $\begingroup$ See my 2nd code. It works. It shows how to do Circular Convolution in Frequency Domain to match Element Wise Multiplication in the Time Domain. Exactly what you wanted. $\endgroup$ – Royi Mar 27 at 13:40
1
$\begingroup$

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));
| improve this answer | |
$\endgroup$
  • $\begingroup$ How different it is from the code I posted before? No need to shift anything. Have a look on my code. $\endgroup$ – Royi Mar 27 at 13:40
  • $\begingroup$ Why do you think there is no need to shift ? $\endgroup$ – Hilmar Mar 27 at 17:39
  • $\begingroup$ Have a look at my code. No need to shift. Just use Circular Convolution. $\endgroup$ – Royi Mar 27 at 18:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.