Difference between revisions of "MATLAB:Transfer Functions"

From PrattWiki
Jump to navigation Jump to search
(Very Basic Analytical)
 
(9 intermediate revisions by the same user not shown)
Line 1: Line 1:
== Example ==
+
== Creation ==
===Analytical===
+
=== Basics ===
Assume you have a transfer function <math>\mathbb{H}(s)</math> for a low-pass <math>RC</math> circuit where you have a function for the ratio of the Laplace transform of the voltage across the capacitor to the voltage across the series combination of a capacitor and a resistor:
+
An easy way to create a transfer function in MATLAB involves using the command
 +
<source lang="matlab">
 +
s=tf('s')
 +
</source>
 +
to create s as a variable and then use s in a line of code to make a transfer function.
 +
=== Advanced ===
 +
The <code>feedback</code> 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,
 +
<source lang="matlab">
 +
T=feedback(G, H)
 +
</source>
 +
is functionally the same as:
 +
<source lang="matlab">
 +
T = G / (1 + G*H)
 +
</source>
 +
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 <code>minreal</code> command:
 +
<source lang="matlab">
 +
T = minreal(G / (1 + G*H))
 +
</source>
 +
 
 +
== Examples ==
 +
===Very Basic Analytical===
 +
Assume you have a transfer function <math>\mathbb{H}(s)</math> for a high-pass <math>RC</math> 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:
 
<center>
 
<center>
 
<math>
 
<math>
 
\begin{align}
 
\begin{align}
\mathbb{H}(s)=\frac{\frac{1}{sC}}{R+\frac{1}{sC}}=\frac{1}{(RC)s+1}
+
\mathbb{H}(s)=\frac{R}{R+\frac{1}{sC}}=\frac{(RC)s}{(RC)s+1}
 +
\end{align}
 +
</math></center>
 +
There are a few different ways to examine the magnitude and phase content of the Fourier version of this transfer function:
 +
<center>
 +
<math>
 +
\begin{align}
 +
\mathbb{H}(j\omega)=\frac{(RC)j\omega}{(RC)j\omega+1}
 
\end{align}
 
\end{align}
 
</math></center>
 
</math></center>
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 document will show how to use MATLAB's <code>tf</code> function to set up and analyze the magnitude and phase of the transfer function of circuit.
+
 
 +
over a range of frequencies.  This example will show how to use MATLAB's <code>tf</code> function to set up and analyze the magnitude and phase of the transfer function of circuit.  It will allow the <code>bode</code> command to generate the plot - including the choice of frequencies over which to plot.
  
 
<source lang="matlab">
 
<source lang="matlab">
Line 17: Line 46:
 
%% Set up transfer function
 
%% Set up transfer function
 
%  Create "s" as a transfer function for use later
 
%  Create "s" as a transfer function for use later
s = tf([1 0], [1]);
+
s = tf('s');
 
%  Use s to generate transfer function for circuit
 
%  Use s to generate transfer function for circuit
H = 1 / (s * R * C + 1);
+
H = (s * R * C) / (s * R * C + 1);
 +
%  Use bode command to analyze transfer function
 +
bode(H);
 +
title('Bode plots for System H');
 +
 
 +
%% Print file
 +
print -dpng bode_1_plot
 +
</source>
 +
 
 +
===More Advanced Analytical===
 +
The next example will show how to use MATLAB's <code>tf</code> 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..
 +
 
 +
<source lang="matlab">
 +
%% 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
 
%  Generate list of frequencies; must use angular frequencies in bode command
 
F = logspace(1, 5, 1000);
 
F = logspace(1, 5, 1000);
Line 42: Line 92:
 
</source>
 
</source>
  
===Theoretical===
+
===Estimates from Time Series Data===
If you have a data set and want to find the theoretical transfer function between two variables in the set, you can have MATLAB come up with a transfer function estimate using the <code>tfestimate</code> command. In the example below, assume you are trying to find an estimate for the transfer function:
+
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 <code>tfestimate</code> command. In the example below, assume you are trying to find an estimate for the transfer function:
 
<center><math>
 
<center><math>
 
\begin{align}
 
\begin{align}
Line 92: Line 142:
 
ylabel('\angle H, rad')
 
ylabel('\angle H, rad')
 
</source>
 
</source>
 
===Comparison Between Theoretical and Experimental ===
 
Take the plotting examples and combine the semilogx arguments.  Use different line styles.  Add a legend.  Comment to self about how easy that was!
 

Latest revision as of 10:46, 16 March 2021

Creation

Basics

An easy way to create a transfer function in MATLAB involves using the command

s=tf('s')

to create s as a variable and then use s in a line of code to make a transfer function.

Advanced

The 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,

T=feedback(G, H)

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 minreal command:

T = minreal(G / (1 + G*H))

Examples

Very Basic Analytical

Assume you have a transfer function \(\mathbb{H}(s)\) for a high-pass \(RC\) 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:

\( \begin{align} \mathbb{H}(s)=\frac{R}{R+\frac{1}{sC}}=\frac{(RC)s}{(RC)s+1} \end{align} \)

There are a few different ways to examine the magnitude and phase content of the Fourier version of this transfer function:

\( \begin{align} \mathbb{H}(j\omega)=\frac{(RC)j\omega}{(RC)j\omega+1} \end{align} \)

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');

%% Print file
print -dpng bode_1_plot

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:

\( \begin{align} \mathbb{H}_{\mbox{est}}(j\omega)&=\frac{\mathbb{V}_{out}(j\omega)}{\mathbb{V}_{in}(j\omega)} \end{align} \)

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;

where:

  • 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 \(T\) is \({t}_{\mbox{end}}-{t}_{\mbox{start}}\) and \(\omega_0=2\pi/T\). 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')