How to plot the continuous convolution result

4 views (last 30 days)
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 Comments
Torsten
Torsten on 21 Feb 2023
How do you want to plot a function that is 0 everywhere and Inf at t=0 ?
YAMENG SONG
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:

Sign in to comment.

Accepted Answer

Paul
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])
  3 Comments
YAMENG SONG
YAMENG SONG on 22 Feb 2023
can i also ask what does the 2nd return f1 represent?
jammer1(f,f1) = rectangularPulse(f1_L,f1_H,f);
Paul
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!