[7.9 MATLAB: FOURIER](#page-13-0) TRANSFORM TOPICS

← Back to LINEAR SYSTEMS AND SIGNALS Overview

7.9 MATLAB: FOURIER TRANSFORM TOPICS

MATLAB is useful for investigating a variety of Fourier transform topics. In this section, a rectangular pulse is used to investigate the scaling property, Parseval’s theorem, essential bandwidth, and spectral sampling. Kaiser window functions are also investigated.

^† In addition to truncation, we need to delay the truncated function by T/2 to render it causal. However, the time delay only adds a linear phase to the spectrum without changing the amplitude spectrum. Thus, to simplify our discussion, we shall ignore the delay.

Figure 7.48 Window-based filter design.

7.9-1 The Sinc Function and the Scaling Property

As shown in Ex. 7.2, the Fourier transform of x(t) = rect(t/τ ) is X(ω) = τ sinc (ωτ/2). To represent X(ω) in MATLAB, a sinc function is first required. As an alternative to the signal processing toolbox function sinc, which computes sinc(x) as sin(πx)/πx, we create our own function that follows the conventions of this book and defines sinc(x) = sin(x)/x.

function [y] = CH7MP1(x) % CH7MP1.m : Chapter 7, MATLAB Program 1 % Function M-file computes the sinc function, y = sin(x)/x.

y(x==0) = 1; y(x~=0) = sin(x(x~=0))./x(x~=0);

The computational simplicity of sinc (x) = sin(x)/x is somewhat deceptive: sin(0)/0 results in a divide-by-zero error. Thus, program CH7MP1 assigns sinc (0) = 1 and computes the remaining values according to the definition. Notice that CH7MP1 cannot be directly replaced by an anonymous function. Anonymous functions cannot have multiple lines or contain certain commands such as =, if, or for. M-files, however, can be used to define an anonymous function. For example, we can represent X(ω) as an anonymous function that is defined in terms of CH7MP1.

>> X = @(omega,tau) tau*CH7MP1(omega*tau/2);

Once we have defined X(ω), it is simple to investigate the effects of scaling the pulse width τ . Consider the three cases τ = 1.0, τ = 0.5, and τ = 2.0.

>> omega = linspace(-4*pi,4*pi,200);
>> plot(omega,X(omega,1),'k-',omega,X(omega,0.5),'k-.',omega,X(omega,2),'k--');
>> grid; axis tight; xlabel('\omega'); ylabel('X(\omega)');
>> legend('Baseline (\tau = 1)','Compressed (\tau = 0.5)',...
>> 'Expanded (\tau = 2.0)');

Figure 7.49 confirms the reciprocal relationship between signal duration and spectral bandwidth: time compression causes spectral expansion, and time expansion causes spectral compression. Additionally, spectral amplitudes are directly related to signal energy. As a signal is compressed, signal energy and thus spectral magnitude decrease. The opposite effect occurs when the signal is expanded.

Figure 7.49 Spectra X(ω) = τ sinc (ωτ/2) for τ = 1.0, τ = 0.5, and τ = 2.0.

7.9-2 Parseval’s Theorem and Essential Bandwidth

Parseval’s theorem concisely relates energy between the time domain and the frequency domain:

$$\int_{-\infty}^{\infty} |x(t)|^2 dt = \frac{1}{2\pi} \int_{-\infty}^{\infty} |X(\omega)|^2 d\omega$$

This too is easily verified with MATLAB. For example, a unit amplitude pulse x(t) with duration τ has energy Ex = τ . Thus,

$$\int_{-\infty}^{\infty} |X(\omega)|^2 \, d\omega = 2\pi \, \tau$$

Letting τ = 1, the energy of X(ω) is computed by using the quad function.

>> X_squared = @(omega, tau) (tau*CH7MP1(omega*tau/2)).^2;
>> quad(X_squared,-1e6,1e6,[],[],1)
    ans = 6.2817

Although not perfect, the result of the numerical integration is consistent with the expected value of 2π ≈ 6.2832. For quad, the first argument is the function to be integrated, the next two arguments are the limits of integration, the empty square brackets indicate default values for special options, and the last argument is the secondary input τ for the anonymous function X_squared. Full format details for quad are available from MATLAB’s help facilities.

A more interesting problem involves computing a signal’s essential bandwidth. Consider, for example, finding the essential bandwidth W, in radians per second, that contains fraction β of the energy of the square pulse x(t). That is, we want to find W such that

$$\frac{1}{2\pi} \int_{-W}^{W} |X(\omega)|^2 d\omega = \beta \tau$$

Program CH7MP2 uses a guess-and-check method to find W.

function [W,E_W] = CH7MP2(tau,beta,tol)
% CH7MP2.m : Chapter 7, MATLAB Program 2
% Function M-file computes essential bandwidth W for square pulse.
% INPUTS: tau = pulse width
% beta = fraction of signal energy desired in W
% tol = tolerance of relative energy error
% OUTPUTS: W = essential bandwidth [rad/s]
% E_W = Energy contained in bandwidth W
W = 0; step = 2*pi/tau; % Initial guess and step values
X_squared = @(omega,tau) (tau*CH7MP1(omega*tau/2)).^2;
E = beta*tau; % Desired energy in W
relerr = (E-0)/E; % Initial relative error is 100 percent
while(abs(relerr) > tol),
   if (relerr>0), % W too small, so...
       W=W+step; % ... increase W by step
   elseif (relerr<0), % W too large, so...
       step = step/2; % ... decrease step and then W

W = W-step;
    end
    E_W = 1/(2*pi)*quad(X_squared,-W,W,[],[],tau);
    relerr = (E - E_W)/E;
end

Although this guess-and-check method is not the most efficient, it is relatively simple to understand: CH7MP2 sensibly adjusts W until the relative error is within tolerance. The number of iterations needed to converge to a solution depends on a variety of factors and is not known beforehand. The while command is ideal for such situations:

while expression,

statements;

end

While the expression is true, the statements are continually repeated.

To demonstrate CH7MP2, consider the 90% essential bandwidth W for a pulse of 1 second duration. Typing [W,E_W]=CH7MP2(1,0.9,0.001) returns an essential bandwidth W = 5.3014 that contains 89.97% of the energy. Reducing the error tolerance improves the estimate. CH7MP2(1,0.9,0.00005) returns an essential bandwidth W = 5.3321 that contains 90.00% of the energy. These essential bandwidth calculations are consistent with estimates presented after Ex. 7.2.

7.9-3 Spectral Sampling

Consider a signal with finite duration τ . A periodic signal xT⁰ (t) is constructed by repeating x(t) every T⁰ seconds, where T⁰ ≥ τ . From Eq. (7.5), we can write the Fourier series coefficients of xT⁰ (t) as Dn = (1/T0)X(n2π/T0). Put another way, the Fourier series coefficients are obtained by sampling the spectrum X(ω).

By using spectral sampling, it is simple to determine the Fourier series coefficients for an arbitrary duty-cycle, square-pulse periodic signal. The square pulse x(t) = rect(t/τ ) has spectrum X(ω) = τ sinc(ωτ /2). Thus, the nth Fourier coefficient of the periodic extension xT⁰ (t) is Dn = (τ/T0)sinc (nπτ/T0). As in Ex. 6.4, τ = π and T⁰ = 2π provide a square-pulse periodic signal. The Fourier coefficients are determined by

>> tau = pi; T_0 = 2*pi; n = [0:10];

D_n = tau/T_0*MS7P1(n*pi*tau/T_0);

>> stem(n,D_n); xlabel('n'); ylabel('D_n');

axis([-0.5 10.5 -0.2 0.55]);

The results, shown in Fig. 7.50, agree with Fig. 6.6b. Doubling the period to T⁰ = 4π effectively doubles the density of spectral samples and halves the spectral amplitude, as shown in Fig. 7.51.

As T⁰ increases, the spectral sampling becomes progressively finer while the amplitude becomes infinitesimal. An evolution of the Fourier series toward the Fourier integral is seen by allowing the period T⁰ to become large. Figure 7.52 shows the result for T⁰ = 40π.

If T⁰ = τ , the signal xT⁰ is a constant and the spectrum should concentrate energy at dc. In this case, the sinc function is sampled at the zero crossings and Dn = 0 for all n not equal to 0. Only the sample corresponding to n = 0 is nonzero, indicating a dc signal, as expected. It is a simple matter to modify the previous code to verify this case.

Figure 7.50 Fourier spectra for τ = π and T⁰ = 2π.

Figure 7.51 Fourier spectra for τ = π and T⁰ = 4π.

Figure 7.52 Fourier spectra for τ = π and T⁰ = 40π.

7.9-4 Kaiser Window Functions

A window function is useful only if it can be easily computed and applied to a signal. The Kaiser window, for example, is flexible but appears rather intimidating:

$$w_K(t) = \begin{cases} \frac{I_0(\alpha\sqrt{1 - 4(t/T)^2})}{I_0(\alpha)} & |t| < T/2\\ 0 & \text{otherwise} \end{cases}$$

Fortunately, the bark of a Kaiser window is worse than its bite! The function I0(x), a zero-order modified Bessel function of the first kind, can be computed according to

$$I_0(x) = \sum_{k=0}^{\infty} \left(\frac{x^k}{2^k k!}\right)^2$$

or, more simply, by using the MATLAB function besseli(0,x). In fact, MATLAB supports a wide range of Bessel functions, including Bessel functions of the first and second kinds (besselj and bessely), modified Bessel functions of the first and second kinds (besseli and besselk), Hankel functions (besselh), and Airy functions (airy).

Program CH7MP3 computes Kaiser windows at times t by using parameters T and α.

function [w_K] = CH7MP3(t,T,alpha)
% CH7MP3.m : Chapter 7, MATLAB Program 3
% Function M-file computes a width-T Kaiser window using parameter alpha.
% Alpha can also be a string identifier: 'rectangular', 'Hamming', or
% 'Blackman'.
% INPUTS: t = independent variable of the window function
% T = window width
% alpha = Kaiser parameter or string identifier
% OUTPUTS: w_K = Kaiser window function
if strncmpi(alpha,'rectangular',1),
   alpha = 0;
elseif strncmpi(alpha,'Hamming',3),
   alpha = 5.4414;
elseif strncmpi(alpha,'Blackman',1),
   alpha = 8.885;
elseif isa(alpha,'char')
   disp('Unrecognized string identifier.'); return
end
w_K = zeros(size(t)); i = find(abs(t)<T/2);
w_K(i) = besseli(0,alpha*sqrt(1-4*t(i).^2/(T^2)))/besseli(0,alpha);

Recall that α = 0, α = 5.4414, and α = 8.885 correspond to rectangular, Hamming, and Blackman windows, respectively. CH7MP3 is written to allow these special-case Kaiser windows to be identified by name rather than by α value. While unnecessary, this convenient feature is achieved with the help of the strncmpi command.

The strncmpi(S1,S2,N) command compares the first N characters of strings S1 and S2, ignoring case. More completely, MATLAB has four variants of string comparison: strcmp, strcmpi, strncmp, and strncmpi. Comparisons are restricted to the first N characters when n is present; case is ignored when i is present. Thus, CH7MP3 identifies any string alpha that starts with the letter r or R as a rectangular window. To prevent confusion with a Hanning window, the first three characters must match to identify a Hamming window. The isa(alpha,‘char’) command determines whether alpha is a character string. MATLAB help documents the many other classes that isa can identify. In CH7MP3, isa is used to terminate execution if a string identifier alpha has not been recognized as one of the three special cases.

Figure 7.53 Special-case, unit-duration Kaiser windows.

Figure 7.53 shows the three special-case, unit-duration Kaiser windows generated by

>> t = [-0.6:.001:0.6]; T = 1;

>> plot(t,CH7MP3(t,T,'r'),'k-',t,CH7MP3(t,T,'ham'),'k-.',t,CH7MP3(t,T,'b'),'k--');

>> axis([-0.6 0.6 -.1 1.1]); xlabel('t'); ylabel('w_K(t)');

>> legend('Rectangular','Hamming','Blackman','Location','EastOutside');