How to plot the continuous convolution result

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)

[Torsten]

How do you want to plot a function that is 0 everywhere and Inf at t=0 ?

[YAMENG SONG]

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:

Answer 1

[Paul]

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

[YAMENG SONG]

thank you so much!

[YAMENG SONG]

can i also ask what does the 2nd return f1 represent?

jammer1(f,f1) = rectangularPulse(f1_L,f1_H,f);

[Paul]

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.