Page 456 - Numerical Methods for Chemical Engineering

1-D Fourier transforms in MATLAB 445

Thus, we can extract all of the independent information of the Fourier transform of a sampled
real signal f (t) from only one half of the discrete Fourier transform values, say the first
half.
Even for a real signal f (t), the Fourier transform F(ω) is complex; therefore, it is
2
common practice to compute the real-valued, nonnegative power spectrum |F(ω)| vs.
2
ω. The quantity |F(ω)| ≥ 0 is a measure of the contribution to f (t) from dynamics with
2
a frequency ω. The symmetry F ∗ = F m leads to the property |F(2ω max − ω n )| =
N−m
2
|F(ω n )| , and thus all independent information in the power spectrum is contained in the
lower-frequency half 0 ≤ ω ≤ ω max .

1-D Fourier transforms in MATLAB

The discrete Fourier transform and its inverse (9.54) are computed in MATLAB using fft and
ifft respectively. The following code demonstrates their use for the signal f (t) = sin(t) +
2 cos(2t), which has two peaks in the power spectrum at ω = 1 and ω = 2. The signal is
sampled at uniform times over the period [0,l2π], such that P = lπ and t = 2P/N,
k
where N = 2 for some integer e.
% generate sampled f(t) data
N = 2ˆ8; P = 10*pi; dt = (2*P)/N;
t val = linspace(0,2*P-dt,N);
f val = sin(t val) + 2*cos(2*t val);
% compute discrete FT
F DFT = fft(f val);
% generate sampled frequencies
omega max = pi/dt; d omega = pi/P;
omega val = linspace(0,2*omega max-d omega,N);
% compute continuous FT and power spectrum
F FT = dt/sqrt(2*pi)*F DFT; F PS = abs(F FT);
% plot signal and power spectrum
figure; subplot(2,1,1); plot(t val,f val);
xlabel(‘t’); ylabel(‘f(t)’);
title(‘f(t) and |F(\omega)|’);
subplot(2,1,2); plot(omega val,F PS);
xlabel(‘\omega’); ylabel(‘|F(\omega)|’);
Plots of the time signal f (t) and |F(ω)| are shown in Figure 9.2. There appear two peaks
at ω = 1 and ω = 2 with a 2:1 relative intensity. These peaks are repeated again at high
frequencies due to the ω-symmetry for a real signal, F ∗ = F m . Thus, we obtain all the
N−m
useful information from the low-frequency half ω ∈ [0,ω max ]. From the discrete Fourier
transform, the time signal can be reconstituted from the inverse discrete FT function ifft,
f 2 = real(ifft(F DFT));

The real() operation is done to remove any near-zero imaginary contributions that are
introduced due to numerical error.