Why the fourier transform of a signal does not produce any graphing?

Hi, I have written a MATLAB code to convert the signal equations from time domain into frequency domain. However, when plotting the graph in frequency domain, there is no output as shown in the screenshot attached. Since my a(t) is a cosine, I believe there should be containing only the real part in frequency domain according to the Fourier Transform Table. Thank you.

%create symbolic functions x, a, b, c, d, e, f_c with independent variable t
syms x a h b c d e f_c f_c1 f_c2 t tau

%create symbolic functions A, B, C, D, and E with independent variable w
syms A B C D E w

x(t) = cos(100*pi*t);
a(t) = x(0.4*t);
h(t) = dirac(t-0.02);
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf);

f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);

c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);

figure
A(w) = fourier(a(t));
subplot(2,2,1), fplot((A(w)))
xlim([-50 50]), ylim([-10 10])
title('Real part of A(\omega)')
xlabel('Frequency, \omega')
ylabel('Real A(\omega)')
grid on

B(w) = fourier(b(t));
subplot(2,2,2), fplot((B(w)))
xlim([-50 50]), ylim([-10 10])
title('Real part of B(\omega)')
xlabel('Frequency, \omega')
ylabel('Real B(\omega)')
grid on

C(w) = fourier(c(t));
subplot(2,2,3), fplot((C(w)))
xlim([-50 50]), ylim([-10 10])
title('Real part of C(\omega)')
xlabel('Frequency, \omega')
ylabel('Real C(\omega)')
grid on

D(w) = fourier(d(t));
subplot(2,2,4), fplot((D(w)))
xlim([-50 50]), ylim([-10 10])
title('Real part of D(\omega)')
xlabel('Frequency, \omega')
ylabel('Real D(\omega)')
grid on

Answer 1

[Santosh Fatale]

Hi Teoh,

The function "dirac" represents the Dirac delta function and not the Kronecker delta function. Thus, it achieves an infinity value at . The function "fplot" does not plot the infinity for the Dirac delta function, and you need to take care of that using the following workaround.

syms t
dirac_t = dirac(t);

t = -5:0.1:5;
dirac_t = dirac(t);
idx = dirac_t == Inf;
dirac_t(idx) = 1; % Assign suitable value instead of Inf.

stem(t, dirac_t); % You can even use plot function.

The modified code to plot real part of variable A is as follows:

%create symbolic functions x, a, b, c, d, e, f_c with independent variable t
syms x(t) a(t) h(t) b(t) c(t) d(t) e f_c(t) f_c1(t) f_c2(t) t tau

%create symbolic functions A, B, C, D, and E with independent variable w
syms A(w) B(w) C(w) D(w) E w

x(t) = cos(100*pi*t)
a(t) = x(0.4*t)
h(t) = dirac(t-0.02)
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf)

f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);

c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);

figure
A(w) = fourier(a(t), w);
w = -45*pi:0.1*pi:45*pi;
subsA = A(w);
% Replace Inf value with suitable value
idx = subsA == Inf;
subsA(idx) = pi; % choose suitable value from the expression of fourier transform.
plot(w,real(subsA));
ylabel("\Re(A(\omega))");
xlabel("\omega");
xticks(-40*pi:20*pi:40*pi);
xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})

Refer documentation of “dirac”, “fourier”, and “real” for more information and illustration examples.

[TEOH CHEE JIN]

Hi, @Santosh Fatale Thank you for your guidance. I believe the same techniques applied to B(w), C(w), D(w) and E(w). Hence, I proceed with the script for B(w) first. However, I encounter the index error in this manner. I have perform manual calculation and I found out that the graph of B(w) should be the same as A(w), but because of the error, the graph is not presentable in the output. How should I correct my script in this cass?

%create symbolic functions x, a, b, c, d, e, f_c with independent variable t
syms x(t) a(t) h(t) b(t) c(t) d(t) e f_c(t) f_c1(t) f_c2(t) t tau

%create symbolic functions A, B, C, D, and E with independent variable w
syms A(w) B(w) C(w) D(w) E(w)

x(t) = cos(100*pi*t);
a(t) = x(0.4*t);
h(t) = dirac(t-0.02);
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf);
f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);
c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);

figure
subplot (2,1,1)
fplot(a(t))
xlim([-0.05 0.05]),ylim([-1.5 1.5])
title ('Time domain of signal a(t)')
xlabel('Time, t')
ylabel('Amplitude, a(t)')
grid on

A(w) = fourier(a(t), w);
w = -45*pi:0.1*pi:45*pi;
subsA = A(w);
% Replace Inf value with suitable value
idx = subsA == Inf;
subsA(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (2,1,2)
plot(w,real(subsA));
ylim([0 3.5])
title ('Time domain of signal A(\omega)')
ylabel("\Re(A(\omega))");
xlabel("\omega");
xticks(-40*pi:20*pi:40*pi);
xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})

figure
subplot(2,1,1)
fplot(b(t))
xlim([-0.05 0.05]),ylim([-1.5 1.5])
title ('Time domain of signal b(t)')
xlabel('Time, t')
ylabel('Amplitude, b(t)')
grid on

B(w) = fourier(b(t), w);
simplify(B(w))
w = -45*pi:0.1*pi:45*pi;
subsB = B(w);
% Replace Inf value with suitable value
idx = subsB == Inf;
subsB(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (2,1,2)
plot(w,real(subsB));
ylim([0 3.5])
title ('Time domain of signal B(\omega)')
ylabel("\Re(B(\omega))");
xlabel("\omega");
xticks(-40*pi:20*pi:40*pi);
xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})

[Santosh Fatale]

Hi Teoh,

You can follow the same procedure for other transforms. Please note that you need to take care of "Infinity" value in resepctive transform.

A few of the pointers to resolve the issue are:

• When we define values for variable "w", it changes its class from "syms" to double. Every time you use a symbolic variable, you must first create it.
• For "A(w)", since there is no complex component in transform, the condition that I considered will work.
• For others, use "abs(subsB) == Inf" condition instead of "subsB == Inf" condition.

The updated sample code is as follows:

syms w;
B(w) = fourier(b(t), w);
simplify(B(w))
w = -45*pi:0.1*pi:45*pi;
subsB = B(w);
% Replace Inf value with suitable value
idx = abs(subsB) == Inf;
subsB(idx) = pi; % choose suitable value from the expression of fourier transform.
subplot (2,1,2)
plot(w,real(subsB));
ylim([0 3.5])
title ('Time domain of signal B(\omega)')
ylabel("\Re(B(\omega))");
xlabel("\omega");
xticks(-40*pi:20*pi:40*pi);
xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})

[Santosh Fatale]

Hi Teoh,

Please check the sample code in my last reply. You need to define the symbolic variable "w" before using it in the "fourier" function.

I also want to bring to your attention that you should assign a suitable value wherever you get the Dirac function in the transform. Weights associated with the respective impulse function could be a suitable value.

[TEOH CHEE JIN]

Hi @Santosh Fatale, for the amplitude of D(w), I have used the command "simplify" and obtained that there are two different amplitudes at different frequency. At w=40pi, the absolute amplitude is 50pi, whereas at w=4960pi and w=5040pi, the amplitudes are 25pi, what should I do to plot the graph showing different amplitudes at different frequency?

%create symbolic functions x, a, b, c, d, e, f_c with independent variable t
syms x(t) a(t) h(t) b(t) c(t) d(t) e f_c(t) f_c1(t) f_c2(t) t tau
%create symbolic functions A, B, C, D, and E with independent variable w
syms A(w) B(w) C(w) D(w) E(w) w
x(t) = cos(100*pi*t);
a(t) = x(0.4*t);
h(t) = dirac(t-0.02);
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf);
f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);
c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);
% f_c1(t)
% f_c2(t)
% figure
% 
% subplot (3,1,1)
% fplot(a(t))
% xlim([-0.05 0.05]),ylim([-1.5 1.5])
% title ('Time domain of signal a(t)')
% xlabel('Time, t')
% ylabel('Amplitude, a(t)')
% grid on
% 
% A(w) = fourier(a(t), w);
% w = -60*pi:0.1*pi:60*pi;
% subsA = A(w);
% % Replace Inf value with suitable value
% idx = subsA == Inf;
% subsA(idx) = pi; % choose suitable value from the expression of fourier transform.
% subplot (3,1,2)
% plot(w,real(subsA));
% ylim([0 4])
% title ('Frequency domain of signal A(\omega)')
% ylabel("\Re(A(\omega))");
% xlabel("\omega");
% xticks(-40*pi:20*pi:40*pi);
% xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})
% 
% subplot (3,1,3)
% plot(w, angle(A(w)))
% title ('Phase spectrum of signal A(\omega)')
% ylabel("\angle(A(\omega))");
% xlabel("\omega");
% 
% 
% figure
% 
% subplot(3,1,1)
% fplot(b(t))
% xlim([-0.05 0.05]),ylim([-1.5 1.5])
% title ('Time domain of signal b(t)')
% xlabel('Time, t')
% ylabel('Amplitude, b(t)')
% grid on
% 
% syms B w
% B(w) = fourier(b(t), w);
% w = -60*pi:0.1*pi:60*pi;
% subsB = B(w);
% % Replace Inf value with suitable value
% idx = abs(subsB) == Inf;
% subsB(idx) = pi; % choose suitable value from the expression of fourier transform.
% subplot (3,1,2)
% plot(w,real(subsB));
% ylim([0 4])
% title ('Frequency domain of signal B(\omega)')
% ylabel("\Re(B(\omega))");
% xlabel("\omega");
% xticks(-60*pi:20*pi:60*pi);
% xticklabels({'-40\pi', '-20\pi', '0', '20\pi', '40\pi'})
% 
% subplot (3,1,3)
% plot(w, angle(B(w)))
% ylim([-4 4])
% title ('Phase spectrum of signal B(\omega)')
% ylabel("\angle(B(\omega))");
% xlabel("\omega");
% figure
% 
% subplot(3,1,1)
% fplot(c(t))
% xlim([-0.05 0.05]),ylim([-12 12])
% title ('Time domain of signal c(t)')
% xlabel('Time, t')
% ylabel('Amplitude, c(t)')
% grid on
% 
% syms C w
% C(w) = fourier(c(t), w);
% simplify(C(w))
% w = -2800*pi:20*pi:2800*pi;
% subsC = C(w);
% % Replace Inf value with suitable value
% idx = abs(subsC) == Inf;
% subsC(idx) = 5*pi; % choose suitable value from the expression of fourier transform.
% subplot (3,1,2)
% plot(w,real(subsC));
% 
% ylim([0 20])
% title ('Frequency domain of signal C(\omega)')
% ylabel("\Re(C(\omega))");
% xlabel("\omega");
% xticks([-2540*pi -2460*pi 0 2460*pi 2540*pi]);
% xticklabels({'-2540\pi', '-2460\pi', '0', '2460\pi', '2540\pi'})
% 
% subplot (3,1,3)
% plot(w, angle(C(w)))
% ylim([-4 4])
% title ('Phase spectrum of signal C(\omega)')
% ylabel("\angle(C(\omega))");
% xlabel("\omega");

figure 
subplot(3,1,1)
fplot(d(t))
xlim([-0.1 0.1]),ylim([-110 110])
title ('Time domain of signal d(t)')
xlabel('Time, t')
ylabel('Amplitude, d(t)')
grid on
syms D w
D(w) = fourier(d(t), w);
w = -5040*pi:10*pi:5040*pi;
simplify(D(w))

ans = 

subsD = D(w);
% Replace Inf value with suitable value
idx = abs(subsD) == Inf;
subsD(idx) = 25*pi; % choose suitable value from the expression of fourier transform.
subplot (3,1,2)
plot(w,real(subsD));
ylim([0 160])
title ('Frequency domain of signal D(\omega)')
ylabel("\Re(D(\omega))");
xlabel("\omega");
subplot (3,1,3)
plot(w, angle(D(w)))
ylim([-4 4])
title ('Phase spectrum of signal D(\omega)')
ylabel("\angle(D(\omega))");
xlabel("\omega");

This is the calculation that I have done to verify the answer.

[TEOH CHEE JIN]

Hi @Santosh Fatale or anyone, any comments to achieve the goal above?

[Paul]

Hi TEOH,

As a general approach to this problem, I'm not sure that evaluating D(w) at a bunch of frequencies and looking for infs is robust. Unless, of course, you know which frequencies to evaluate, in which case you only need to do the evalution at those specific frequencies.

Assuming that path, then I guess the substitution should be whatever the coefficient is on the dirac. But doing an evaluation of D(w) at the dirac frequencies may not be the best approach. For this problem, the symbolic evaluation at w = 40*pi is

%create symbolic functions x, a, b, c, d, e, f_c with independent variable t
syms x(t) a(t) h(t) b(t) c(t) d(t) e f_c(t) f_c1(t) f_c2(t) t tau
%create symbolic functions A, B, C, D, and E with independent variable w
syms A(w) B(w) C(w) D(w) E(w) w
x(t) = cos(100*pi*t);
a(t) = x(0.4*t);
h(t) = dirac(t-0.02);
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf);
f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);
c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);
syms D w
D(w) = fourier(d(t), t, w);
D(40*sym(pi))

ans = 

From there, it may be possible to isolate the -50*pi.

However, the actual coefficient on dirac(w - 40*pi) term is

d40 = int(D(w),w,39*pi,41*pi)

d40 = 

which has both a real and imaginary part

vpa(d40,10)

ans = 

real(d40)

ans = 

imag(d40)

ans = 

Though the magnitude is 50*pi

simplify(abs(d40))

ans = 

I don't think this result is consistent with your by-hand result, which only captures the magnitude of the coefficient.

25*sym(pi)*2*dirac(w-40*sym(pi))*exp(-1j*0.8*w)

ans = 

simplify(ans)

ans = 

[TEOH CHEE JIN]

Hi @Paul @Santosh Fatale, I would like to present this in graphing method, instead of symbolic equation. Symbolic equation is just the way to verify my by-hand calculation. So any solution for this?

[TEOH CHEE JIN]

I am intesrested in the real axis and its magnitude only.

[Paul]

"I am intesrested in the real axis and its magnitude only."

The code posted above includes the lines

plot(w,real(subsD));
plot(w, angle(D(w)))

so I thought you were interested in the real and imaginary parts of D(w).

If you already have the result and just need a graphical representation, maybe use stem ?

Here's an example

syms t w real
f(t) = 5*sin(t);
F(w) = fourier(f(t),t,w)

F(w) = 

We see that the amplitude is 5*pi at w = +- 1

stem([-1 1],5*pi*[1 1])
axis('padded')

[TEOH CHEE JIN]

@Santosh Fatale @Paul there are 50pi and 25pi spectrums respectively, so if I use your stem method, the output graph will still be one amplitude only, either i set it to be 25pi or 50pi. I wish to have both 25pi and 50pi in one graph.

[Paul]

D(w) is the sum of six, scaled dirac delta functions and nothing else. Consequently, for this particular problem, you can use the method I showed above to pick off the coefficients of each delta function and then do whatever you want as far as plotting the results.

syms x(t) a(t) h(t) b(t) c(t) d(t) e f_c(t) f_c1(t) f_c2(t) t tau
%create symbolic functions A, B, C, D, and E with independent variable w
syms A(w) B(w) C(w) D(w) E(w) w
x(t) = cos(100*pi*t);
a(t) = x(0.4*t);
h(t) = dirac(t-0.02);
b(t) = int(a(tau)*h(t-tau), 'tau', -inf, inf);
f_c(t) = 10*cos(2500*pi*t);
f_c1(t) = f_c(t);
f_c2(t) = f_c(t);
c(t) = b(t)*f_c1(t);
d(t) = c(t)*f_c2(t);
syms D w
D(w) = fourier(d(t), t, w)

D(w) = 

% compute the dirac coefficients
wfreq = sym([40 4960 5040]);
wfreq = [wfreq -wfreq]*sym(pi);
for ii = 1:6
    coefficient(ii) = int(D(w),w,wfreq(ii)-1,wfreq(ii)+1);
end

% plot
figure
subplot(211)
stem(wfreq,abs(coefficient))
ylabel('magnitude')
subplot(212);
stem(wfreq,angle(coefficient));
ylabel('phase (rad)')
xlabel('frequency (rad/sec)')

[TEOH CHEE JIN]

@Paul @Santosh Fatale Hi, thanks for your suggestion. I appreciate it very much. Do you have any script that can remain the impulse at w=40pi using LPF via FilterDesign/FDA tool? I wish to remove the impulses on left and right side.

[Paul]

Not familiar with that tool, but I suspect you can find examples in the doc or here on Answers on how to design an LPF.

[TEOH CHEE JIN]

@Paul， do you know any staff who is expert in programming filter tool?

[Paul]

I do not.