An easy way to create a transfer function in MATLAB involves using the command
to create s as a variable and then use s in a line of code to make a transfer function.
feedback command in MATLAB takes plant and output sensor transfer functions (G and H in the Nise book's paradigm) and produces the overall transfer function assuming negative feedback. Specifically, if G and H are defined as variables,
is functionally the same as:
T = G / (1 + G*H)
The difference is, in the former case, MATLAB automatically checks for pole-zero cancellation. To make sure that MATLAB always checks for pole-zero cancellation, use the
T = minreal(G / (1 + G*H))
Very Basic Analytical
Assume you have a transfer function for a high-pass circuit where you have a function for the ratio of the Laplace transform of the voltage across the resistor to the voltage across the series combination of a capacitor and a resistor:
There are a few different ways to examine the magnitude and phase content of the Fourier version of this transfer function:
over a range of frequencies. This example will show how to use MATLAB's
tf function to set up and analyze the magnitude and phase of the transfer function of circuit. It will allow the
bode command to generate the plot - including the choice of frequencies over which to plot.
%% Circuit constants R = 10000; C = 22e-9; %% Set up transfer function % Create "s" as a transfer function for use later s = tf('s'); % Use s to generate transfer function for circuit H = (s * R * C) / (s * R * C + 1); % Use bode command to analyze transfer function bode(H); title('Bode plots for System H');
More Advanced Analytical
The next example will show how to use MATLAB's
tf function to set up and analyze the magnitude and phase of the transfer function of circuit. It will also show the code to modify how the information is plotted, including changing the frequency domain over which the information is plotted..
%% Circuit constants R = 10000; C = 22e-9; %% Set up transfer function % Create "s" as a transfer function for use later s = tf('s'); % Use s to generate transfer function for circuit H = (s * R * C) / (s * R * C + 1); % Generate list of frequencies; must use angular frequencies in bode command F = logspace(1, 5, 1000); Omega = 2 * pi * F; % Use bode command to analyze transfer function [HMag, HPhase, HOmega] = bode(H, Omega); HMag = squeeze(HMag); HPhase = squeeze(HPhase); %% Make plot figure(1); clf % Magnitude plot on top subplot(2,1,1) semilogx(HOmega, 20*log10(HMag), 'k-') xlabel('\omega, rad/s') ylabel('|H|, dB') % Phase plot on bottom subplot(2,1,2) semilogx(HOmega, HPhase, 'k-') xlabel('\omega, rad/s') ylabel('\angle H, rad')
Estimates from Time Series Data
If you have a data set and want to find an estimated experimental transfer function between two variables in the set, you can have MATLAB come up with a transfer function estimate using the
tfestimate command. In the example below, assume you are trying to find an estimate for the transfer function:
to generate a magnitude plot and a phase plot of an experimentally
determined transfer function. MATLAB's
tfestimate will produce a numerical
estimate of the magnitude and phase of a transfer function given an
input signal, an output signal, and possibly other information. The specific form of this command is:
[EstH, EstF] = tfestimate(Vin, Vout, , , , Fs) EstMag = abs(EstH); EstPhase = angle(EstH); EstOmega = EstF*2*pi;
- Vin - vector containing the input voltage values
- Vout - vector containing the output voltage values
- Fs - the sampling frequency for the voltage values, in Hz
- EstH - vector containing complex numbers that contain the amplitudes and phases of the estimate of the transfer function
- EstF - vector containing the corresponding frequencies, in Hz, for the magnitudes and phases stored in EstH
- EstMag - magnitude of the transfer function estimate
- EstPhase - phase of the transfer function estimate
- EstOmega - corresponding angular frequencies for the estimate values
The square brackets in the third through fifth arguments are placeholders for parameters whose default values are fine for this experiment. Essentially, this command will determine the frequency content of the input and of the output using subsections of the data; it will then compute the magnitude ratio and phase difference between the input and the output and provide an approximation of the transfer function at particular frequencies contained in EstF.
One issue that comes into play is you must make sure the input signal has energy at as many frequencies as possible to give
tfestimate values to work with. Using a single-frequency cosine as an input, for example, might lead to a disaster - if that frequency happens to exactly hit one of the frequencies
tfestimate is using - since there would only be one input frequency on which to base an estimate for the response to all frequencies. This would be similar to estimating the acoustics of a concert hall by hitting a single tuning fork.
It is important to note, however, that what MATLAB is really doing is assuming that your input is one period's worth of some periodic input - meaning is and . Furthermore, the
tfestimate command breaks the signal up into several windows, further reducing the possibility that a problem of this nature will occur. If, however, the input signal does not have any energy in the higher frequency ranges, the
tfestimate program will not be able to compensate. You should therefore make sure that whatever signal gets used contains a wide range of frequencies. Ideally,
the input signal would be an impulse function since that contains all
frequencies at equal amounts. Barring that - and the Laws of Nature do a pretty nice job of barring that - generating a noise signal or a sweeping frequency signal with energy in the bands of interest works well.
To make a plot, you can use code similar to:
%% Make plot figure(1); clf % Magnitude plot on top subplot(2, 1, 1) semilogx(EstOmega, 20*log10(EstMag), 'k-') xlabel('\omega, rad/s') ylabel('|H|, dB') % Phase plot on bottom subplot(2,1,2) semilogx(EstOmega, EstPhase, 'k-') xlabel('\omega, rad/s') ylabel('\angle H, rad')