Ploting the magnitude of the laplace transformed signal against the frequency range

Hello.

I'm trying to plot an amplitude spectre of a signal using a Laplace transform. In the example code given below I choose variable "s" to have a real value of 0 and an imaginary part that falls into the range of -15i to 15i (from -15 rad/s to 15 rad/s).

The code yields an error message "Division by zero" even though I thought I handled all the cases where a division by zero could occur.

Any help would be appreaciated.

syms t s
x(t) = 3*cos(5*t);
X_Laplace(s) = laplace(x(t));

% (3*s)/(s^2 + 25)
% Numerator root 0
% Denumerator roots +5i and -5i

N = 300; %Number of points in the plot
Frequencies = zeros(1,300);
Laplace_Values = zeros(1,300);


Current_Frequency = -15;
Increment = 0.1;
for i = 1:1:N
    Frequencies(1,i) = Current_Frequency;
    if Frequencies(1,i) == -5
        Laplace_Values(1,i) = 1;
    elseif Frequencies(1,i) == 5
        Laplace_Values(1,i) = 1;
    else
        Laplace_Values(1,i) = abs(X_Laplace(Frequencies(1,i)*1i));
    end

    Current_Frequency = Current_Frequency+Increment;
end

figure
plot(Frequencies,Laplace_Values)
grid on

The error messages are:

Error using symengine
Division by zero.

Error in sym/subs>mupadsubs (line 168)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);

Error in sym/subs (line 153)
    G = mupadsubs(F,X,Y);

Error in symfun/feval>evalScalarFun (line 42)
y = subs(formula(Ffun), Fvars, inds);

Error in symfun/feval (line 28)
    varargout{1} = evalScalarFun(F, Fvars, varargin);

Error in indexing (line 210)
            B = feval(A, inds{:});

Error in untitled (line 23)
        Laplace_Values(1,i) = abs(X_Laplace(Frequencies(1,i)*1i));

Answer 1

[Star Strider]

Since you are doing actual numeric calculations, one option is to use matlabFunction to convert ‘X_Laplace’ from a symbolic function to an anonymous function. This has the advantage of being more efficient, and also allowing NaN or Inf values to exist in the ‘Laplace_Values’ vector.

Try this —

syms t s
x(t) = 3*cos(5*t);
X_Laplace(s) = laplace(x(t))

X_Laplace(s) = 

X_Laplace = matlabFunction(X_Laplace)

X_Laplace = function_handle with value:
    @(s)(s.*3.0)./(s.^2+2.5e+1)

% (3*s)/(s^2 + 25)
% Numerator root 0
% Denumerator roots +5i and -5i

N = 300; %Number of points in the plot
Frequencies = zeros(1,300);
Laplace_Values = zeros(1,300);


Current_Frequency = -15;
Increment = 0.1;
for i = 1:1:N
    Frequencies(1,i) = Current_Frequency;
    if Frequencies(1,i) == -5
        Laplace_Values(1,i) = 1;
    elseif Frequencies(1,i) == 5
        Laplace_Values(1,i) = 1;
    else
        Laplace_Values(1,i) = abs(X_Laplace(Frequencies(1,i)*1i));
    end

    Current_Frequency = Current_Frequency+Increment;
end

figure
plot(Frequencies,Laplace_Values)
grid on

EDIT — Corrected typographical errors.

.

[Arvydas Kromelis]

Hey, that really works. Cheers!

[Star Strider]

Thank you!