Regarding plotting transfer function of a system equation

1,090 views
Skip to first unread message

Ram Prasadh

unread,
Sep 25, 2012, 10:17:43 AM9/25/12
to fre...@googlegroups.com
hi 

I just started using freemat. I want to know if there is a possibility to find the transfer function of any system and plot it on a bode plot. In matlab ,We have some in-built functions like 'tf' , 'bode' for making the above said things. I think I did not find anything like that here. Is there any way I can do it?? 

Jonathan Weaver

unread,
Sep 26, 2012, 6:16:34 AM9/26/12
to fre...@googlegroups.com
It's not built in, but it's one of the first things I wanted to do.  So, I wrote three functions: tf.m, bode.m, and rlocus.m to create a transfer function response given the numerator and denominator polynomial and a vector of the frequencies I was interested in, and then plot a bode plot and a root locus plot.  These are not as fully functional as MATLAB, but they did what I needed at the time.  To use, try

% plot bode for 1/(s + 1)
N = [0 1]
D = [1 1]
w = logspace(1,6);
resp = tf(N, D, w)
bode(resp, w)

% plot root locus for 1/(s + 1)
k = logspace(1,2)
rlocus(N, D, k)

Functions are inline below:

function [TF] = tf(N, D, w)
    jw = j.*w;
    if size(N) == size(D)
        a = size(D);
        if(a(1) == 1)
            Wv = jw.^0;
            for i = 2:a(2)
                Wv = [jw.^(i-1); Wv];
            end
            TF = (N*Wv)./(D*Wv);
        else
            error('N and D should be column vectors');
        end
    else
        error('N and D should be the same size');
    end
end


function [A] = bode(TF,w)
    subplot(2,1,1);
    semilogx(w./2./pi, 20*log10(abs(TF)));
    title('Magnitude', 'fontsize', 14);
    ylabel('Gain (dB)');
    a = gca
    set(a, 'xminorgrid','on');
    grid on;
    subplot(2,1,2);
    semilogx(w./2./pi, angle(TF));
    title('Phase', 'fontsize', 14);
    ylabel('Angle (Radians)');
    xlabel('Frequency (Hz)');
    a = gca
    set(a, 'xminorgrid','on');
    grid on;
end

function [A] = rlocus(N, D, k)
    plot(real(roots(D)), imag(roots(D)), 'bx');
    hold on;
    plot(real(roots(N)), imag(roots(N)), 'ro');
    C = ones(size(k'))*D + k'*N;
    B = size(C);
    A = zeros(B(1), B(2)-1);
    for i = 1:B(1)
        A(i,:) = roots(C(i,:));
    end
    for i = 1:B(2)-1
        plot(real(A(:,i)), imag(A(:,i)),'.:');
    end
    hold off;
end
Reply all
Reply to author
Forward
0 new messages