[2.7 MATLAB: M-FILES](#page-8-0)

← Back to LINEAR SYSTEMS AND SIGNALS Overview

2.7 MATLAB: M-FILES

M-files are stored sequences of MATLAB commands and help simplify complicated tasks. There are two types of M-file: script and function. Both types are simple text files and require a .m filename extension.

Although M-files can be created by using any text editor, MATLAB’s built-in editor is the preferable choice because of its special features. As with any program, comments improve the readability of an M-file. Comments begin with the % character and continue through the end of the line.

An M-file is executed by simply typing the filename (without the .m extension). To execute, M-files need to be located in the current directory or any other directory in the MATLAB path. New directories are easily added to the MATLAB path by using the addpath command.

^† This follows directly from Eq. (2.49) with ^{λ = σ +}jω⁰ and ^{= σ}.

2.7-1 Script M-Files

Script files, the simplest type of M-file, consist of a series of MATLAB commands. Script files record and automate a series of steps, and they are easy to modify. To demonstrate the utility of a script file, consider the operational amplifier circuit shown in Fig. 2.25.

The system’s characteristic modes define the circuit’s behavior and provide insight regarding system behavior. Using ideal, infinite gain difference amplifier characteristics, we first derive the differential equation that relates output y(t) to input x(t). Kirchhoff’s current law (KCL) at the node shared by R¹ and R³ provides

$$\frac{x(t) - v(t)}{R_3} + \frac{y(t) - v(t)}{R_2} + \frac{0 - v(t)}{R_1} - C_2 \dot{v}(t) = 0$$

KCL at the inverting input of the op amp gives

$$\frac{v(t)}{R_1} + C_1 \dot{y}(t) = 0$$

Combining and simplifying the KCL equations yield

$$\ddot{y}(t) + \frac{1}{C_2} \left\{ \frac{1}{R_1} + \frac{1}{R_2} + \frac{1}{R_3} \right\} \dot{y}(t) + \frac{1}{R_1 R_2 C_1 C_2} y(t) = -\frac{1}{R_1 R_3 C_1 C_2} x(t)$$

which is the desired constant coefficient differential equation. Thus, the characteristic equation is given by

$$\lambda^2 + \frac{1}{C_2} \left\{ \frac{1}{R_1} + \frac{1}{R_2} + \frac{1}{R_3} \right\} \lambda + \frac{1}{R_1 R_2 C_1 C_2} = (a_0 \lambda^2 + a_1 \lambda + a_2) = 0 \tag{2.50}$$

The roots λ¹ and λ² of Eq. (2.50) establish the nature of the characteristic modes e^λ1*^t* and e^λ2*^t* .

As a first case, assign nominal component values of R¹ = R² = R³ = 10 k and C¹ = C² = 1 µF. A series of MATLAB commands allows convenient computation of the roots λ = [λ1;λ2]. Although λ can be determined using the quadratic equation, MATLAB’s roots command is more convenient. The roots command requires an input vector that contains the polynomial coefficients in descending order. Even if a coefficient is zero, it must still be included in the vector.

Figure 2.25 Operation-amplifier circuit.

% CH2MP1.m : Chapter 2, MATLAB Program 1
% Script M-file determines characteristic roots of op-amp circuit.
% Set component values:
R = [1e4, 1e4, 1e4]; C = [1e-6, 1e-6];
% Determine coefficients for characteristic equation:
A = [1, (1/R(1)+1/R(2)+1/R(3))/C(2), 1/(R(1)*R(2)*C(1)*C(2))];
% Determine characteristic roots:
lambda = roots(A);

A script file is created by placing these commands in a text file, which in this case is named CH2MP1.m. While comment lines improve program clarity, their removal does not affect program functionality. The program is executed by typing

>> CH2MP1

After execution, all the resulting variables are available in the workspace. For example, to view the characteristic roots, type

lambda lambda = -261.8034 -38.1966

Thus, the characteristic modes are simple decaying exponentials: e−261.8034*^t* and e−38.1966*^t* .

Script files permit simple or incremental changes, thereby saving significant effort. Consider what happens when capacitor C¹ is changed from 1.0 µF to 1.0 nF. Changing CH2MP1.m so that C = [1e-9, 1e-6] allows computation of the new characteristic roots:

>> CH2MP1
>> lambda
    lambda = 1.0e+003 *
             -0.1500 + 3.1587i
             -0.1500 - 3.1587i

Perhaps surprisingly, the characteristic modes are now complex exponentials capable of supporting oscillations. The imaginary portion of λ dictates an oscillation rate of 3158.7 rad/s or about 503 Hz. The real portion dictates the rate of decay. The time expected to reduce the amplitude to 25% is approximately t = ln 0.25/Re(λ) ≈ 0.01 second.

2.7-2 Function M-Files

It is inconvenient to modify and save a script file each time a change of parameters is desired. Function M-files provide a sensible alternative. Unlike script M-files, function M-files can accept input arguments as well as return outputs. Functions truly extend the MATLAB language in ways that script files cannot.

Syntactically, a function M-file is identical to a script M-file except for the first line. The general form of the first line is

function [output1, …, outputN] = filename(input1, …, inputM)

For example, consider modification of CH2MP1.m to make function CH2MP2.m. Component values are passed to the function as two separate inputs: a length-3 vector of resistor values and a length-2 vector of capacitor values. The characteristic roots are returned as a 2×1 complex vector.

function [lambda] = CH2MP2(R,C)
% CH2MP2.m : Chapter 2, MATLAB Program 2
% Function M-file finds characteristic roots of op-amp circuit.
% INPUTS: R = length-3 vector of resistances
% C = length-2 vector of capacitances
% OUTPUTS: lambda = characteristic roots
% Determine coefficients for characteristic equation:
A = [1, (1/R(1)+1/R(2)+1/R(3))/C(2), 1/(R(1)*R(2)*C(1)*C(2))];
% Determine characteristic roots:
lambda = roots(A);

As with script M-files, function M-files execute by typing the name at the command prompt. However, inputs must also be included. For example, CH2MP2 easily confirms the oscillatory modes of the preceding example.

>> lambda = CH2MP2([1e4, 1e4, 1e4],[1e-9, 1e-6])
    lambda = 1.0e+003 *
             -0.1500 + 3.1587i
             -0.1500 - 3.1587i

Although scripts and functions have similarities, they also have distinct differences that are worth pointing out. Scripts operate on workspace data; either functions must be supplied data through inputs or they must create their own data. Unless passed as an output, variables and data created by functions remain local to the function; variables or data generated by scripts are global and are added to the workspace. To emphasize this point, consider polynomial coefficient vector A, which is created and used in both CH2MP1.m and CH2MP2.m. Following execution of function CH2MP2, the variable A is not added to the workspace. Following execution of script CH2MP1, however, A is available in the workspace. Recall, the workspace is easily viewed by typing either who or whos.

2.7-3 For-Loops

Real resistors and capacitors never exactly equal their nominal values. Suppose that the circuit components are measured as R¹ = 10.322 k, R² = 9.952 k, R³ = 10.115 k, C¹ = 1.120 nF, and C² = 1.320 µF. These values are consistent with the 10 and 25% tolerance resistor and capacitor values commonly and readily available. CH2MP2.m uses these component values to calculate the new values of λ.

216 CHAPTER 2 TIME-DOMAIN ANALYSIS OF CONTINUOUS-TIME SYSTEMS

>> lambda = CH2MP2([10322,9592,10115],[1.12e-9, 1.32e-6])
    lambda = 1.0e+003 *
             -0.1136 + 2.6113i
             -0.1136 - 2.6113i

Now the natural modes oscillate at 2611.3 rad/s or about 416 Hz. Decay to 25% amplitude is expected in t = ln 0.25/(−113.6) ≈ 0.012 second. These values, which differ significantly from the nominal values of 503 Hz and t ≈ 0.01 second, warrant a more formal investigation of the effect of component variations on the locations of the characteristic roots.

It is sensible to look at three values for each component: the nominal value, a low value, and a high value. Low and high values are based on component tolerances. For example, a 10% 1 k resistor could have an expected low value of 1000(1 − 0.1) = 900 and an expected high value of 1000(1+0.1) = 1100 . For the five passive components in the design, 35 = 243 permutations are possible.

Using either CH2MP1.m or CH2MP2.m to solve each of the 243 cases would be very tedious and boring. For-loops help automate repetitive tasks such as this. In MATLAB, the general structure of a for statement is

for variable = expression, statement, ..., statement, end

Five nested for-loops, one for each passive component, are required for the present example.

% CH2MP3.m : Chapter 2, MATLAB Program 3
% Script M-file determines characteristic roots over a range of component
  values.
% Pre-allocate memory for all computed roots:
lambda = zeros(2,243);
% Initialize index to identify each permutation:
p=0;
for R1 = 1e4*[0.9,1.0,1.1],
    for R2 = 1e4*[0.9,1.0,1.1],
        for R3 = 1e4*[0.9,1.0,1.1],
            for C1 = 1e-9*[0.75,1.0,1.25],
                for C2 = 1e-6*[0.75,1.0,1.25],
                    p = p+1;
                    lambda(:,p) = CH2MP2([R1 R2 R3],[C1 C2]);
                end
            end
        end
    end
end
plot(real(lambda(:)),imag(lambda(:)),'kx',...
    real(lambda(:,1)),imag(lambda(:,1)),'kv',...
    real(lambda(:,end)),imag(lambda(:,end)),'k^')
xlabel('Real'),ylabel('Imaginary')
legend('Char. Roots','Min. Val. Roots','Max. Val. Roots','Location','West');

Figure 2.26 Effect of component values on characteristic root locations.

The command lambda = zeros(2,243) preallocates a 2×243 array to store the computed roots. When necessary, MATLAB performs dynamic memory allocation, so this command is not strictly necessary. However, preallocation significantly improves script execution speed. Notice also that it would be nearly useless to call script CH2MP1 from within the nested loop; script file parameters cannot be changed during execution.

The plot instruction is quite long. Long commands can be broken across several lines by terminating intermediate lines with three dots (…). The three dots tell MATLAB to continue the present command to the next line. Black x’s locate roots of each permutation. The command lambda(:) vectorizes the 2 × 243 matrix lambda into a 486 × 1 vector. This is necessary in this case to ensure that a proper legend is generated. Because of loop order, permutation p = 1 corresponds to the case of all components at the smallest values and permutation p = 243 corresponds to the case of all components at the largest values. This information is used to separately highlight the minimum and maximum cases using down-triangles () and up-triangles (), respectively. In addition to terminating each for loop, end is used to indicate the final index along a particular dimension, which eliminates the need to remember the particular size of a variable. An overloaded function, such as end, serves multiple uses and is typically interpreted based on context.

The graphical results provided by CH2MP3 are shown in Fig. 2.26. Between extremes, root oscillations vary from 365 to 745 Hz and decay times to 25% amplitude vary from 6.2 to 12.7 ms. Clearly, this circuit’s behavior is quite sensitive to ordinary component variations.

2.7-4 Graphical Understanding of Convolution

MATLAB graphics effectively illustrate the convolution process. Consider the case of y(t) = x(t)∗ h(t), where x(t) = 1.5 sin(πt)(u(t) − u(t − 1)) and h(t) = 1.5(u(t) − u(t − 1.5)) − u(t − 2) + u(t − 2.5). Program CH2MP4 steps through the convolution over the time interval (−0.25 ≤ t ≤ 3.75).

% CH2MP4.m : Chapter 2, MATLAB Program 4
% Script M-file graphically demonstrates the convolution process.
figure(1) % Create figure window and make visible on screen
u = @(t) 1.0*(t>=0);
x = @(t) 1.5*sin(pi*t).*(u(t)-u(t-1));
h = @(t) 1.5*(u(t)-u(t-1.5))-u(t-2)+u(t-2.5);
dtau = 0.005; tau = -1:dtau:4;
ti = 0; tvec = -.25:.1:3.75;
y = NaN*zeros(1,length(tvec)); % Pre-allocate memory
for t = tvec,
    ti = ti+1; % Time index
    xh = x(t-tau).*h(tau); lxh = length(xh);
    y(ti) = sum(xh.*dtau); % Trapezoidal approximation of convolution integral
    subplot(2,1,1),plot(tau,h(tau),'k-',tau,x(t-tau),'k--',t,0,'ok');
    axis([tau(1) tau(end) -2.0 2.5]);
    patch([tau(1:end-1);tau(1:end-1);tau(2:end);tau(2:end)],...
        [zeros(1,lxh-1);xh(1:end-1);xh(2:end);zeros(1,lxh-1)],...
        [.8 .8 .8],'edgecolor','none');
    xlabel('\tau'); title('h(\tau) [solid], x(t-\tau) [dashed], h(\tau)x(t-\tau) [gray]');
    c = get(gca,'children'); set(gca,'children',[c(2);c(3);c(4);c(1)]);
    subplot(2,1,2),plot(tvec,y,'k',tvec(ti),y(ti),'ok');
    xlabel('t'); ylabel('y(t) = \int h(\tau)x(t-\tau) d\tau');
    axis([tau(1) tau(end) -1.0 2.0]); grid;
    drawnow;
end

At each step, the program plots h(τ ), x(t − τ ), and shades the area h(τ )x(t − τ ) gray. This gray area, which reflects the integral of h(τ )x(t − τ ), is also the desired result, y(t). Figures 2.27, 2.28, and 2.29 display the convolution process at times t of 0.75, 2.25, and 2.85 seconds, respectively. These figures help illustrate how the regions of integration change with time. Figure 2.27 has limits of integration from 0 to (t = 0.75). Figure 2.28 has two regions of integration, with limits (t −1 = 1.25) to 1.5 and 2.0 to (t = 2.25). The last plot, Fig. 2.29, has limits from 2.0 to 2.5.

Several comments regarding CH2MP4 are in order. The command figure(1) opens the first figure window and, more important, makes sure it is visible. Anonymous functions are used to represent the functions u(t), x(t), and h(t). NaN, standing for not-a-number, usually results from operations such as 0/0 or ∞−∞. MATLAB refuses to plot NaN values, so preallocating y(t) with NaNs ensures that MATLAB displays only values of y(t) that have been computed. As its name suggests, length returns the length of the input vector. The subplot(a,b,c) command partitions the current figure window into an a-by-b matrix of axes and selects axes c for use. Subplots facilitate graphical comparison by allowing multiple axes in a single figure window. The patch command is used to create the gray-shaded area for h(τ )x(t − τ ). In CH2MP4, the get and set commands are used to reorder plot objects so that the gray area does not obscure other lines. Details of the patch, get, and set commands, as used in CH2MP4, are somewhat advanced and are not pursued here.† MATLAB also prints most Greek letters if the Greek name is preceded by a backslash () character. For example, \tau in the xlabel command produces the symbol τ in the plot’s axis label. Similarly, an integral sign is produced by \int. Finally, the drawnow

^† Interested students should consult the MATLAB help facilities for further information. Actually, the get and set commands are extremely powerful and can help modify plots in almost any conceivable way.

Figure 2.27 Graphical convolution at step t = 0.75 second.

Figure 2.28 Graphical convolution at step t = 2.25 seconds.

command forces MATLAB to update the graphics window for each loop iteration. Although slow, this creates an animation-like effect. Replacing drawnow with the pause command allows users to manually step through the convolution process. The pause command still forces the graphics window to update, but the program will not continue until a key is pressed.

Figure 2.29 Graphical convolution at step t = 2.85 seconds.