# How to plot the continuous convolution result

4 views (last 30 days)
YAMENG SONG on 21 Feb 2023
Edited: Paul on 22 Feb 2023
Hi,
I was trying to plot my fourier transform of a continuous function using fplot, but it doesn't give me the expected result (either a straight line or blank),
can anyone check where I went wrong?
thanks a lot!
there's my code:
clear; clc; clear all;
syms f f1 t
%jammer 1
f1_L = f1 - 2.5*10^6;
f1_H = f1 + 2.5*10^6;
jammer1 = rectangularPulse(f1_L,f1_H,f);
jammer1_time_domain = simplify(ifourier(jammer1, f, t))
jammer1_time_domain =
% jammer 2
f2 = f1 + 12.5*10^6;
f2_L = f2 - 2.5*10^6;
f2_H = f2 + 2.5*10^6;
jammer2 = rectangularPulse(f2_L,f2_H,f);
jammer2_time_domain = simplify(ifourier(jammer2, f, t))
jammer2_time_domain =
% signal
f_signal = f1 - 12.5*10^6;
% convolution J1 and J1 on fre domain = multiplicaiton on time domian
fv_1_1_time_domain = jammer1_time_domain * jammer1_time_domain
fv_1_1_time_domain =
fv_1_1_freq_domain = simplify(ifourier(jammer1_time_domain * jammer1_time_domain, f, t));
% convolution J1 and J1 and J2
fv_1_1_2 = simplify(fourier(jammer2_time_domain * fv_1_1_time_domain, f, t))
fv_1_1_2 =
fplot(fv_1_1_freq_domain)
##### 2 CommentsShow 1 older commentHide 1 older comment
YAMENG SONG on 21 Feb 2023
I think I might do the convolution wrong?
I'd like to convolute jammer_1 * jammer_1*jammer_2 (which are rectangular functions) on freq domain, they why I define:
clear; clc; clear all;
syms f f1 t
f1_L = f1 - 2.5*10^6;
f1_H = f1 + 2.5*10^6;
f2 = f1 + 12.5*10^6;
f2_L = f2 - 2.5*10^6;
f2_H = f2 + 2.5*10^6;
jammer1 = rectangularPulse(f1_L,f1_H,f);
jammer2 = rectangularPulse(f2_L,f2_H,f);
since convolution on freq domian equlals to multiplicaiton on time domian, so I used ifourier to convert them back.
jammer1_time_domain = simplify(ifourier(jammer1, f, t))
jammer2_time_domain = simplify(ifourier(jammer2, f, t))
then I multipliy them together, and converted back to freq domain.
% convolution J1 and J1 on fre domain = multiplicaiton on time domian
fv_1_1_time_domain = jammer1_time_domain * jammer1_time_domain
fv_1_1_2 = simplify(fourier(jammer2_time_domain * fv_1_1_time_domain, f, t))
It shouldn't be a delta function as I got, that's why I'm so confused about.
It should be like what circled out in the picture:

Paul on 22 Feb 2023
Hi YAMENG,
Have you considered computing the convolution integrals directly?
Note scaling to get them all to fit on the same plot.
syms f f1 t
f1_L = f1 - 2.5*10^6;
f1_H = f1 + 2.5*10^6;
f2 = f1 + 12.5*10^6;
f2_L = f2 - 2.5*10^6;
f2_H = f2 + 2.5*10^6;
jammer1(f,f1) = rectangularPulse(f1_L,f1_H,f);
jammer2(f,f1) = rectangularPulse(f2_L,f2_H,f);
syms u
j11(f,f1) = int(jammer1(u,f1)*jammer1(f-u,f1),u,-inf,inf);
j211(f,f1) = int(jammer2(u,f1)*j11(f-u,f1),u,-inf,inf);
figure
hold on
fplot(jammer1(f,10e6),[0 6e7])
fplot(jammer2(f,10e6),[0 6e7])
fplot(j11(f,10e6)/1e6,[0 6e7])
fplot(j211(f,10e6)/1e12,[0 6e7])
Paul on 22 Feb 2023
Edited: Paul on 22 Feb 2023
The value of f1 was not specified, so I carried it along as an additional argument. Doing so made it easier to make the plots for a specific value of f1, which I arbitrarily chose as 10e6. There are alternative approaches.
Also, the fourier approach can work as well, if you give Matlab a little bit of help via rewrite. Here, I won't use f1 as an additional argument.
syms f f1 t real
f1_L = f1 - 2.5*10^6;
f1_H = f1 + 2.5*10^6;
f2 = f1 + 12.5*10^6;
f2_L = f2 - 2.5*10^6;
f2_H = f2 + 2.5*10^6;
jammer1(f) = rectangularPulse(f1_L,f1_H,f);
jammer2(f) = rectangularPulse(f2_L,f2_H,f);
J1(t) = rewrite(simplify(ifourier(jammer1(f),f,t)),'exp')
J1(t) =
J2(t) = rewrite(simplify(ifourier(jammer2(f),f,t)),'exp')
J2(t) =
j211(f) = fourier(J2(t)*J1(t)*J1(t),t,f)
j211(f) =
Divide by 1e12 to compare to previous result.
fplot(subs(j211(f),f1,10e6)/1e12,[0 6e7],'Color',[4.940000000000000e-01 1.840000000000000e-01 5.560000000000000e-01])
The reason this result is smaller is that the default paramters for fourier as defined in sympref assume that the frequency is in rad/sec so the ifourier include a 1/(2*pi) out front. Because we used ifourier twice we need to multiply this result by (2*pi)^2, which I think will make the peak of the plot match the peak of the result computed via the convolution integral.
Also, the original code had an error in the last line. Repeating from the Question:
clear all;
syms f f1 t
f1_L = f1 - 2.5*10^6;
f1_H = f1 + 2.5*10^6;
f2 = f1 + 12.5*10^6;
f2_L = f2 - 2.5*10^6;
f2_H = f2 + 2.5*10^6;
jammer1 = rectangularPulse(f1_L,f1_H,f);
jammer2 = rectangularPulse(f2_L,f2_H,f);
jammer1_time_domain = simplify(ifourier(jammer1, f, t))
jammer1_time_domain =
jammer2_time_domain = simplify(ifourier(jammer2, f, t))
jammer2_time_domain =
% convolution J1 and J1 on fre domain = multiplicaiton on time domian
fv_1_1_time_domain = jammer1_time_domain * jammer1_time_domain;
This line is incorrect. The f and t arguments are backwards
fv_1_1_2 = simplify(fourier(jammer2_time_domain * fv_1_1_time_domain, f, t))
fv_1_1_2 =
It should be
fv_1_1_2 = simplify(fourier(jammer2_time_domain * fv_1_1_time_domain, t, f))
fv_1_1_2 =
But, as can be seen, the result is not useful when jammer*_time_domain are expressed in terms of sin and cos instead of in terms of complex exponentials.
It looks like all of the terms of fourier(sin(sigma)/t^3) could be resolved. See this Answer.