Compute Fourier and Inverse Fourier Transforms

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

F[f](w)=f(x)eiwxdx,

and the inverse Fourier transform (IFT) as

F1[f](x)=12πf(w)eiwxdw.

This documentation refers 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 indicate 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. This documentation refers 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 + x2))–1:

syms x
cauchy = 1/(pi*(1+x^2));
fcauchy = fourier(cauchy)
fcauchy = 
(pi*exp(-w)*heaviside(w) + pi*heaviside(-w)*exp(w))/pi
fcauchy = expand(fcauchy)
fcauchy = 
exp(-w)*heaviside(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/(x*i - 1) - 1/(x*i + 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

d4ydx4+kEIy=1EIδ(x), <x<.

Here, E represents elasticity of the beam (railroad track), I is the "beam constant," and k is the spring (road bed) stiffness. The shock force on the right side of the differential equation is modeled by the Dirac Delta function δ(x). The Dirac Delta function has the following important property:

f(xy)δ(y)dy=f(x).

A definition of the Dirac Delta function is

δ(x)=limnnχ(1/2n,1/2n)(x),

where

χ(1/2n,1/2n)(x)={1for 12n<x<12n0otherwise.

Let Y(w) = F[y(x)](w) and Δ(w) = F[δ(x)](w). Indeed, try the command fourier(dirac(x), x, w). The Fourier transform turns differentiation into exponentiation, and, in particular,

F[d4ydx4](w)=w4Y(w).

To see a demonstration of this property, try this

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

which returns

ans =
w^4*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.

Applying the Fourier transform to the differential equation above yields the algebraic equation

(w4+kEI)Y(w)=Δ(w),

or

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

where

G(w)=1w4+kEI=F[g(x)](w)

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 the differential equation, observe that the ratio

KEI

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

KEI=1024.

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

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

and see

/                   /    / pi       \                   / pi
| sqrt(2) exp(-4 x) | sin| -- + 4 x | heaviside(x) - cos| --
\                   \    \  4       /                   \  4

         \                             \ \
   + 4 x | exp(8 x) (heaviside(x) - 1) | |/512
         /                             / /

Notice that g contains the Heaviside distribution

H(x)={1forx>00forx<01/2forx=0.

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

y(x)=(δg)(x)=g(xy)δ(y)dy=g(x).

by the special property of the Dirac Delta function. To plot this function, substitute the domain of x into y(x), using the subs command. 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.

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

Was this topic helpful?