Learn more about Symbolic Math Toolbox

Integral Transforms and Z-Transforms

On this page…
The Fourier and Inverse Fourier Transforms The Laplace and Inverse Laplace Transforms The Z– and Inverse Z–transforms

The Fourier and Inverse Fourier Transforms

The Fourier transform of a function f(x) is defined as

and the inverse Fourier transform (IFT) as

We refer to this formulation as the Fourier transform of f with respect to x as a function of w. Or, more concisely, the Fourier transform of f with respect to x at w. Mathematicians often use the notation F[f] to denote the Fourier transform of f. In this setting, the transform is taken with respect to the independent variable of f (if f = f(t), then t is the independent variable; f = f(x) implies that x is the independent variable, etc.) at the default variable w. We refer to F[f] as the Fourier transform of f at w and F^–1[f] is the IFT of f at x. See fourier and ifourier in the reference pages for tables that show the Symbolic Math Toolbox commands equivalent to various mathematical representations of the Fourier and inverse Fourier transforms.

For example, consider the Fourier transform of the Cauchy density function, (π(1 + x²))^–1:

syms x
cauchy = 1/(pi*(1+x^2));
fcauchy = fourier(cauchy)

fcauchy = 
((pi*heaviside(w))/exp(w) + pi*heaviside(-w)*exp(w))/pi

fcauchy = expand(fcauchy)

fcauchy = 
heaviside(w)/exp(w) + heaviside(-w)*exp(w)

ezplot(fcauchy)

The Fourier transform is symmetric, since the original Cauchy density function is symmetric.

To recover the Cauchy density function from the Fourier transform, call ifourier:

finvfcauchy = ifourier(fcauchy)

finvfcauchy = 
-(1/(i*x - 1) - 1/(i*x + 1))/(2*pi)

simplify(finvfcauchy)

ans = 
1/(pi*(x^2 + 1))

An application of the Fourier transform is the solution of ordinary and partial differential equations over the real line. Consider the deformation of an infinitely long beam resting on an elastic foundation with a shock applied to it at a point. A "real world" analogy to this phenomenon is a set of railroad tracks atop a road bed.

The shock could be induced by a pneumatic hammer blow.

The differential equation idealizing this physical setting is

Here, E represents elasticity of the beam (rail road track), I is the "beam constant," and k is the spring (road bed) stiffness. The shock force on the right hand side of the differential equation is modeled by the Dirac Delta function δ(x). If you are unfamiliar with δ(x), you may be surprised to learn that (despite its name), it is not a function at all. Rather, δ(x) is an example of what mathematicians call a distribution. The Dirac Delta function (named after the physicist Paul Dirac) has the following important property

A definition of the Dirac Delta function is

where

You can evaluate the Dirac Delta function at a point (say) x = 3, using the commands

syms x
del = sym('dirac(x)');
vpa(subs(del,x,3))

which return

ans =
0.0

Returning to the differential equation, let Y(w) = F[y(x)](w) and Δ(w) = F[δ(x)](w). Indeed, try the command fourier(del,x,w). The Fourier transform turns differentiation into exponentiation, and, in particular,

To see a demonstration of this property, try this

syms w x
fourier(diff(sym('y(x)'), x, 4), x, w)

which returns

ans =
w^4*transform::fourier(y(x), x, -w)

Note that you can call the fourier command with one, two, or three inputs (see the reference pages for fourier). With a single input argument, fourier(f) returns a function of the default variable w. If the input argument is a function of w, fourier(f) returns a function of t. All inputs to fourier must be symbolic objects.

We now see that applying the Fourier transform to the differential equation above yields the algebraic equation

or

Y(w)  =  Δ(w)G(w),

where

for some function g(x). That is, g is the inverse Fourier transform of G:

g(x) = F^–1[G(w)](x)

The Symbolic Math Toolbox counterpart to the IFT is ifourier. This behavior of ifourier parallels fourier with one, two, or three input arguments (see the reference pages for ifourier).

Continuing with the solution of our differential equation, we observe that the ratio

is a relatively "large" number since the road bed has a high stiffness constant k and a rail road track has a low elasticity E and beam constant I. We make the simplifying assumption that

This is done to ease the computation of F ^–1[G(w)](x). Proceeding, we type

G = 1/(w^4 + 1024);
g = ifourier(G, w, x);
g = simplify(g);
pretty(g)

and see

 
   1/2    / pi       \                 1/2                  / pi       \ 
  2    sin| -- + 4 x | heaviside(x)   2    heaviside(-x) sin| -- - 4 x | exp(4 x) 
          \ 4        /                                      \ 4        / 
  --------------------------------- + ------------------------------------------- 
            512 exp(4 x)                                  512

Notice that g contains the Heaviside distribution

Since Y is the product of Fourier transforms, y is the convolution of the transformed functions. That is, F[y] = Y(w) = Δ(w) G(w) = F[δ] F[g] implies

by the special property of the Dirac Delta function. To plot this function, we must substitute the domain of x into y(x), using the subs command.

XX = -3:0.05:3;
YY = double(subs(g, x, XX));
plot(XX, YY)
title('Beam Deflection for a Point Shock')
xlabel('x'); ylabel('y(x)');

The resulting graph

shows that the impact of a blow on a beam is highly localized; the greatest deflection occurs at the point of impact and falls off sharply from there. This is the behavior we expect from experience.

Back to Top

The Laplace and Inverse Laplace Transforms

The Laplace transform of a function f(t) is defined as

while the inverse Laplace transform (ILT) of f(s) is

where c is a real number selected so that all singularities of f(s) are to the left of the line s =  c. The notation L[f] denotes the Laplace transform of f at s. Similarly, L^–1[f] is the ILT of f at t.

The Laplace transform has many applications including the solution of ordinary differential equations/initial value problems. Consider the resistance-inductor-capacitor (RLC) circuit below.

Let Rj and Ij, j = 1, 2, 3 be resistances (measured in ohms) and currents (amperes), respectively; L be inductance (henrys), and C be capacitance (farads); E(t) be the electromotive force, and Q(t) be the charge.

By applying Kirchhoff's voltage and current laws, Ohm's Law, and Faraday's Law, you can arrive at the following system of simultaneous ordinary differential equations.

Let's solve this system of differential equations using laplace. We will first treat the R_j, L, and C as (unknown) real constants and then supply values later on in the computation.

syms R1 R2 R3 L C real
dI1 = sym('diff(I1(t),t)'); dQ = sym('diff(Q(t),t)');
I1 = sym('I1(t)'); Q = sym('Q(t)');
syms t s
E = sin(t);  % Voltage
eq1 = dI1 + R2*dQ/L - (R2 - R1)*I1/L;
eq2 = dQ - (E - Q/C)/(R2 + R3) - R2*I1/(R2 + R3);

At this point, we have constructed the equations in the MATLAB workspace. An approach to solving the differential equations is to apply the Laplace transform, which we will apply to eq1 and eq2. Transforming eq1 and eq2

L1 = laplace(eq1,t,s)
L2 = laplace(eq2,t,s)

returns

L1 = 
s*laplace(I1(t), t, s) - I1(0) 
+ ((R1 - R2)*laplace(I1(t), t, s))/L 
- (R2*(Q(0) - s*laplace(Q(t), t, s)))/L

L2 =
s*laplace(Q(t), t, s) - Q(0) 
- (R2*laplace(I1(t), t, s))/(R2 + R3) - (C/(s^2 + 1) 
- laplace(Q(t), t, s))/(C*(R2 + R3))

Now we need to solve the system of equations L1 = 0, L2 = 0 for laplace(I1(t),t,s) and laplace(Q(t),t,s), the Laplace transforms of I₁ and Q, respectively. To do this, we need to make a series of substitutions. For the purposes of this example, use the quantities R1 = 4 Ω (ohms), R2 = 2 Ω, R3 =  3 Ω, C = 1/4 farads, L = 1.6 H (henrys), I1(0) = 15 A (amperes), and Q(0) = 2 A/sec. Substituting these values in L1

syms LI1 LQ
NI1 = subs(L1,{R1,R2,R3,L,C,'I1(0)','Q(0)'}, ...
      {4,2,3,1.6,1/4,15,2})

returns

NI1 =
s*laplace(I1(t), t, s) + (5*s*laplace(Q(t), t, s))/4
 + (5*laplace(I1(t), t, s))/4 - 35/2

The substitution

NQ = subs(L2,{R1,R2,R3,L,C,'I1(0)','Q(0)'},{4,2,3,1.6,1/4,15,2})

returns

NQ =
s*laplace(Q(t), t, s) - 1/(5*(s^2 + 1)) 
+ (4*laplace(Q(t), t, s))/5 - (2*laplace(I1(t), t, s))/5 - 2

To solve for laplace(I1(t),t,s) and laplace(Q(t),t,s), we make a final pair of substitutions. First, replace the strings 'laplace(I1(t),t,s)' and 'laplace(Q(t),t,s)' by the syms LI1 and LQ, using

NI1 =...
 subs(NI1,{'laplace(I1(t),t,s)','laplace(Q(t),t,s)'},{LI1,LQ})

to obtain

NI1 =
(5*LI1)/4 + LI1*s + (5*LQ*s)/4 - 35/2

Collecting terms

NI1 = collect(NI1,LI1)

gives

NI1 =
(s + 5/4)*LI1 + (5*LQ*s)/4 - 35/2

A similar string substitution

NQ = ...
subs(NQ,{'laplace(I1(t),t,s)','laplace(Q(t),t,s)'},{LI1,LQ})

yields

NQ =
(4*LQ)/5 - (2*LI1)/5 + LQ*s - 1/(5*(s^2 + 1)) - 2

which, after collecting terms,

NQ = collect(NQ,LQ)

gives

NQ =
(s + 4/5)*LQ - (2*LI1)/5 - 1/(5*(s^2 + 1)) - 2

Now, solving for LI1 and LQ

[LI1, LQ] = solve(NI1, NQ, LI1, LQ)

we obtain

LI1 = 
(300*s^3 + 280*s^2 + 295*s + 280)/(20*s^4 + 51*s^3 + 40*s^2 + 51*s + 20)
 
LQ =
(40*s^3 + 190*s^2 + 44*s + 195)/(20*s^4 + 51*s^3 + 40*s^2 + 51*s + 20)

To recover I1 and Q we need to compute the inverse Laplace transform of LI1 and LQ. Inverting LI1

I1 = ilaplace(LI1, s, t)

produces

I1 =
(15*(cosh((1001^(1/2)*t)/40) 
- (293*1001^(1/2)*sinh((1001^(1/2)*t)/40))/21879))/exp((51*t)/40) 
- (5*sin(t))/51

Inverting LQ

Q = ilaplace(LQ, s, t)

yields

Q = 
(4*sin(t))/51 - (5*cos(t))/51 + (107*(cosh((1001^(1/2)*t)/40) 
+ (2039*1001^(1/2)*sinh((1001^(1/2)*t)/40))/15301))/(51*exp((51*t)/40))

Now let's plot the current I1(t) and charge Q(t) in two different time domains, 0 ≤ t ≤ 10 and 5 ≤ t ≤ 25. The statements

subplot(2,2,1); ezplot(I1,[0,10]);                       
title('Current'); ylabel('I1(t)'); grid                  
subplot(2,2,2); ezplot(Q,[0,10]);                        
title('Charge'); ylabel('Q(t)'); grid                    
subplot(2,2,3); ezplot(I1,[5,25]);                       
title('Current'); ylabel('I1(t)'); grid
text(7,0.25,'Transient'); text(16,0.125,'Steady State');
subplot(2,2,4); ezplot(Q,[5,25]);                        
title('Charge'); ylabel('Q(t)'); grid  
text(7,0.25,'Transient'); text(15,0.16,'Steady State'); 

generate the desired plots

Note that the circuit's behavior, which appears to be exponential decay in the short term, turns out to be oscillatory in the long term. The apparent discrepancy arises because the circuit's behavior actually has two components: an exponential part that decays rapidly (the "transient" component) and an oscillatory part that persists (the "steady-state" component).

Back to Top

The Z– and Inverse Z–transforms

The (one-sided) z-transform of a function f(n) is defined as

The notation Z[f] refers to the z-transform of f at z. Let R be a positive number so that the function g(z) is analytic on and outside the circle |z| = R. Then the inverse z-transform (IZT) of g at n is defined as

The notation Z^–1[f] means the IZT of f at n. The Symbolic Math Toolbox commands ztrans and iztrans apply the z-transform and IZT to symbolic expressions, respectively. See ztrans and iztrans for tables showing various mathematical representations of the z-transform and inverse z-transform and their Symbolic Math Toolbox counterparts.

The z-transform is often used to solve difference equations. In particular, consider the famous "Rabbit Problem." That is, suppose that rabbits reproduce only on odd birthdays (1, 3, 5, 7, ...). If p(n) is the rabbit population at year n, then p obeys the difference equation

p(n+2) = p(n+1) + p(n), p(0) = 1, p(1) = 2.

We can use ztrans to find the population each year p(n). First, we apply ztrans to the equations

pn = sym('p(n)');
pn1 = sym('p(n+1)');
pn2 = sym('p(n+2)');
syms n z
eq = pn2 - pn1 - pn;
Zeq = ztrans(eq, n, z)

to obtain

Zeq = 
z*p(0) - z*ztrans(p(n), n, z) - z*p(1) + z^2*ztrans(p(n), n, z) 
   - z^2*p(0) - ztrans(p(n), n, z)

Next, replace 'ztrans(p(n), n, z)' with Pz and insert the initial conditions for p(0) and p(1).

syms Pz
Zeq = subs(Zeq,{'ztrans(p(n), n, z)', 'p(0)', 'p(1)'}, {Pz, 1, 2})

to obtain

Zeq =
Pz*z^2 - z - Pz*z - Pz - z^2

Collecting terms

eq = collect(Zeq, Pz)

yields

eq =
(z^2 - z - 1)*Pz - z^2 - z

Now solve for Pz

P = solve(eq, Pz)

to obtain

P =
-(z^2 + z)/(- z^2 + z + 1)

To recover p(n), we take the inverse z-transform of P.

p = iztrans(P, z, n);
p = simple(p)

The result is a bit complicated, but explicit:

p = 
(3*5^(1/2)*(1/2 - 5^(1/2)/2)^(n - 1))/5 - 
   (3*5^(1/2)*(5^(1/2)/2 + 1/2)^(n - 1))/5 + 
   (4*(-1)^n*cos(n*(pi/2 + i*asinh(1/2))))/i^n

This result can be used as is. But to simplify this result even further, use the MuPAD rewrite function:

p = feval(symengine,'rewrite',p ,'exp');
p = simple(p)

p = 
(3*5^(1/2)*(1/2 - 5^(1/2)/2)^(n - 1))/5 - ...
  (3*5^(1/2)*(5^(1/2)/2 + 1/2)^(n - 1))/5 + ...
  2*exp(pi*i*n)*(2/(5^(1/2) + 1))^n + 2*(5^(1/2)/2 + 1/2)^n

Replace the expression exp(pi*i*n) with (-1)^n:

p = subs(p,exp(pi*i*n),(-1)^n);
p = simple(p)

p = 
((5^(1/2)/2 + 1/2)^n*(3*5^(1/2) + 5) - ...
   (1/2 - 5^(1/2)/2)^n*(3*5^(1/2) + 15) + ...
   20*(-1)^n*(2/(5^(1/2) + 1))^n)/10

pretty(p)

 
  /  1/2     \n                /      1/2 \n 
  | 5      1 |      1/2        | 1   5    |      1/2                n /    2     \n 
  | ---- + - |  (3 5    + 5) - | - - ---- |  (3 5    + 15) + 20 (-1)  | -------- | 
  \  2     2 /                 \ 2    2   /                           |  1/2     | 
                                                                      \ 5    + 1 / 
  --------------------------------------------------------------------------------- 
                                         10

Finally, plot p:

m = 1:10;
y = double(subs(p,n,m));
plot(m,y,'rO')
title('Rabbit Population');
xlabel('years'); ylabel('p');
grid on

to show the growth in rabbit population over time.

References

Andrews, L.C., Shivamoggi, B.K., Integral Transforms for Engineers and Applied Mathematicians, Macmillan Publishing Company, New York, 1986

Crandall, R.E., Projects in Scientific Computation, Springer-Verlag Publishers, New York, 1994

Strang, G., Introduction to Applied Mathematics, Wellesley-Cambridge Press, Wellesley, MA, 1986

Back to Top

Solving Equations

Special Functions of Applied Mathematics

© 1984-2009- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS