[6.7 MATLAB: FOURIER](#page-12-0) SERIES APPLICATIONS

← Back to LINEAR SYSTEMS AND SIGNALS Overview

6.7 MATLAB: FOURIER SERIES APPLICATIONS

Computational packages such as MATLAB simplify the Fourier-based analysis, design, and synthesis of periodic signals. MATLAB permits rapid and sophisticated calculations, which promote practical application and intuitive understanding of the Fourier series.

6.7-1 Periodic Functions and the Gibbs Phenomenon

It is sufficient to define any T0-periodic function over the interval (0 ≤ t < T0). For example, consider the 2π-periodic function given by

$$x(t) = \begin{cases} t/A & 0 \le t < A \\ 1 & A \le t < \pi \\ 0 & \pi \le t < 2\pi \\ x(t + 2\pi) & \text{otherwise} \end{cases}$$

Although similar to a square wave, x(t) has a linearly rising edge of width A, where (0 < A < π). As A → 0, x(t) approaches a square wave; as A → π, x(t) approaches a type of sawtooth wave.

In MATLAB, the mod command helps represent periodic functions such as x(t).

x = @(t,A) mod(t,2*pi)/A.*(mod(t,2*pi)<A)+((mod(t,2*pi)>=A)&(mod(t,2*pi)<pi));

Sometimes referred to as the signed remainder after division, mod(t,2*pi) returns the value t modulo 2π. Thought of another way, the mod operator appropriately shifts t into [0,T0), where x(t) is conveniently defined. Notice also that an anonymous function can be defined with multiple input variables.

The exponential Fourier series coefficients for x(t) (see Prob. 6.3-2) are given by

$$D_n = \begin{cases} \frac{2\pi - A}{4\pi} & n = 0\\ \frac{1}{2\pi n} \left( \frac{e^{-jnA} - 1}{nA} + je^{-jn\pi} \right) & \text{otherwise} \end{cases}$$

Since x(t) is real, D−ⁿ = D^{∗ n}. Truncating the Fourier series at |n| = N yields the approximation

$$x(t) \approx x_N(t) = D_0 + \sum_{n=1}^{N} \left( D_n e^{jnt} + D_n^* e^{-jnt} \right)$$

\n(6.57)

For a user-specified N, program CH6MP1 uses Eq. (6.57) to compute xN(t) over (−π/4 ≤ t < 2π +π/4).

function [x_N,t] = CH6MP1(A,N);
% CH6MP1.m : Chapter 6, MATLAB Program 1
% Function M-file approximates x(t) using Fourier series truncated at |n|=N
% INPUTS: A = width of rising edge
% N = largest harmonic of truncated Fourier series
% OUTPUTS: x_N = Nth harmonic truncated Fourier series

% t = time vector for x_N

% Define FS coefficients for signal x(t)
D = @(n) 1/(2*pi*n)*((exp(-1j*n*A)-1)/(n*A) + 1j*exp(-1j*n*pi));
% Construct truncated FS approximation of x(t) using N harmonics
t = linspace(-pi/4,2*pi+pi/4,10000); % Time vector exceeds one period.
x_N = (2*pi-A)/(4*pi)*ones(size(t)); % Compute dc term
for n = 1:N, % Compute N remaining terms
   x_N = x_N+real(D(n)*exp(1j*n*t) + conj(D(n))*exp(-1j*n*t));
end

Although theoretically not required, the real command ensures that small computer round-off errors do not cause a complex-valued result.

Using program CH6MP1 with A = π/2 and N = 20, Fig. 6.29 compares x(t) and x20(t).

A = pi/2; [x_20,t] = CH6MP1(A,20); >> plot(t,x_20,‘k’,t,x(t,A),‘k:’); axis([-pi/4,2*pi+pi/4,-0.1,1.1]); >> xlabel(‘t’); ylabel(‘x_{20}(t)’);

As expected, the falling edge is accompanied by the overshoot that is characteristic of the Gibbs phenomenon.

Increasing N to 100, as shown in Fig. 6.30, improves the approximation but does not reduce the overshoot.

>> [x_100,t] = CH6MP1(A,100);

plot(t,x_100,‘k’,t,x(t,A),‘k:’); axis([-pi/4,2*pi+pi/4,-0.1,1.1]);

>> xlabel('t'); ylabel('x_{100}(t)');

Reducing A to π/64 produces a curious result. For N = 20, both the rising and falling edges are accompanied by roughly 9% of overshoot, as shown in Fig. 6.31. As the number of terms is increased, overshoot persists only in the vicinity of jump discontinuities. For xN(t), increasing N decreases the overshoot near the rising edge but not near the falling edge. Remember that it is a

Figure 6.29 Comparison of x20(t) and x(t) when A = π/2.

Figure 6.30 Comparison of x100(t) and x(t) when A = π/2.

Figure 6.31 Comparison of x20(t) and x(t) when A = π/64.

Figure 6.32 Comparison of x100(t) and x(t) when A = π/64.

true jump discontinuity that causes the Gibbs phenomenon. A continuous signal, no matter how sharply it rises, can always be represented by a Fourier series at every point within any small error by increasing N. This is not the case when a true jump discontinuity is present. Figure 6.32 illustrates this behavior using N = 100.

6.7-2 Optimization and Phase Spectra

Although magnitude spectra typically receive the most attention, phase spectra are critically important in some applications. Consider the problem of characterizing the frequency response of an unknown system. By applying sinusoids one at a time, the frequency response is empirically measured one point at a time. This process is tedious at best. Applying a superposition of many sinusoids, however, allows simultaneous measurement of many points of the frequency response. Such measurements can be taken by a spectrum analyzer equipped with a transfer function mode or by applying Fourier analysis techniques, which are discussed in later chapters.

A multitone test signal m(t) is constructed as a superposition of N real sinusoids

$$m(t) = \sum_{n=1}^{N} M_n \cos{(\omega_n t + \theta_n)}$$

where Mn and θ*ⁿ* establish the relative magnitude and phase of each sinusoidal component. It is sensible to constrain all gains to be equal, Mn = M for all n. This ensures equal treatment at each point of the measured frequency response. Although the value M is normally chosen to set the desired signal power, we set M = 1 for convenience.

While not required, it is also sensible to space the sinusoidal components uniformly in frequency.

$$m(t) = \sum_{n=1}^{N} \cos(n\omega_0 t + \theta_n)$$

\n(6.58)

Another sensible alternative, which spaces components logarithmically in frequency, is treated in Prob. 6.7-4.

Equation (6.58) is now a truncated compact-form Fourier series with a flat magnitude spectrum. Frequency resolution and range are set by ω⁰ and N, respectively. For example, a 2 kHz range with a resolution of 100 Hz requires ω⁰ = 2π100 and N = 20. The only remaining unknowns are the θn.

While it is tempting to set θ*ⁿ* = 0 for all n, the results are quite unsatisfactory. MATLAB helps demonstrate the problem by using ω⁰ = 2π100 and N = 20 sinusoids, each with a peak-to-peak voltage of 1 volt.

>> m = @(theta,t,omega) sum(cos(omega*t+theta*ones(size(t))));
>> N = 20; omega = 2*pi*100*[1:N]'; theta = zeros(size(omega));
>> t = linspace(-0.01,0.01,10000);
>> plot(t,m(theta,t,omega),'k'); xlabel('t [sec]'); ylabel('m(t) [volts]');

As shown in Fig. 6.33, θ*ⁿ* = 0 causes each sinusoid to constructively add. The resulting 20 volt peak can saturate system components, such as operational amplifiers operating with ±12 volt rails. To improve signal performance, the maximum amplitude of m(t) over t needs to be reduced.

One way to reduce maxt(|m(t)|) is to reduce M, the strength of each component. Unfortunately, this approach reduces the system’s signal-to-noise ratio and ultimately degrades measurement quality. Therefore, reducing M is not a smart decision. The phases θn, however, can be adjusted to reduce maxt(|m(t)|) while preserving signal power. In fact, since θ*ⁿ* = 0 maximizes maxt(|m(t)|), just about any other choice of θ*ⁿ* will improve the situation. Even a random choice should improve performance.

Figure 6.33 Test signal m(t) with θ*ⁿ* = 0.

As with any computer, MATLAB cannot generate truly random numbers. Rather, it generates pseudo-random numbers. Pseudo-random numbers are deterministic sequences that appear to be random. The particular sequence of numbers that is realized depends entirely on the initial state of the pseudo-random number generator. Setting the generator’s initial state to a known value allows a “random” experiment with reproducible results. The command rng(0) initializes the state of the pseudo-random number generator to a known condition of zero, and the MATLAB command rand(a,b) generates an a-by-b matrix of pseudo-random numbers that are uniformly distributed over the interval (0, 1). Radian phases occupy the wider interval (0, 2π ), so the results from rand need to be appropriately scaled.

rng(0); theta_rand0 = 2*pi*rand(N,1);

Next, we recompute and plot m(t) using the randomly chosen θn.

>> m_rand0 = m(theta_rand0,t,omega);

>> plot(t,m_rand0,'k'); axis([-0.01,0.01,-10,10]);

>> xlabel('t [sec]'); ylabel('m(t) [volts]');

set(gca,‘ytick’,[min(m_rand0),max(m_rand0)]); grid on;

For a vector input, the min and max commands return the minimum and maximum values of the vector. Using these values to set y axis tick marks makes it easy to identify the extreme values of the m(t). As seen from Fig. 6.34, the maximum amplitude is now 7.6307, which is significantly smaller than the maximum of 20 when θ*ⁿ* = 0.

Randomly chosen phases suffer a fatal fault: there is little guarantee of optimal performance. For example, repeating the experiment with rng(5) produces a maximum magnitude of 8.2399 volts, as shown in Fig. 6.35. This value is significantly higher than the previous maximum of 7.6307 volts. Clearly, it is better to replace a random solution with an optimal solution.

What constitutes “optimal”? Many choices exist, but desired signal criteria naturally suggest that optimal phases minimize the maximum magnitude of m(t) over all t. To find these optimal phases, MATLAB’s fminsearch command is useful. First, the function to be minimized, called the objective function, is defined.

maxmagm = @(theta,t,omega) max(abs(sum(cos(omega*t+theta*ones(size(t))))));

666 CHAPTER 6 CONTINUOUS-TIME SIGNAL ANALYSIS: THE FOURIER SERIES

Figure 6.34 Test signal m(t) with random θ*ⁿ* found by using rng(0).

Figure 6.35 Test signal m(t) with random θ*ⁿ* found by using rand(‘state’,1).

The anonymous function argument order is important; fminsearch uses the first input argument as the variable of minimization. To minimize over θ, as desired, θ must be the first argument of the objective function maxmagm.

Next, the time vector is shortened to include only one period of m(t).

t = linspace(0,0.01,401);

A full period ensures that all values of m(t) are considered; the short length of t helps ensure that functions execute quickly. An initial value of θ is randomly chosen to begin the search.

>> rng(0); theta_init = 2*pi*rand(N,1);
>> theta_opt = fminsearch(maxmagm,theta_init,[],t,omega);

Notice that fminsearch finds the minimizer to maxmagm over θ by using an initial value theta_init. Most numerical minimization techniques are capable of finding only local minima, and fminsearch is no exception. As a result, fminsearch does not always produce a unique solution. The empty square brackets indicate no special options are requested, and the remaining ordered arguments are secondary inputs for the objective function. Full format details for fminsearch are available from MATLAB’s help facilities.

Figure 6.36 Test signal m(t) with optimized phases.

Figure 6.36 shows the phase-optimized test signal. The maximum magnitude is reduced to a value of 5.3632 volts, which is a significant improvement over the original peak of 20 volts.

Although the signals shown in Figs. 6.33 through 6.36 look different, they all possess the same magnitude spectra. The signals differ only in phase spectra. It is interesting to investigate the similarities and differences of these signals in ways other than graphs and mathematics. For example, is there an audible difference between the signals? For computers equipped with sound capability, the MATLAB sound command can be used to find out.

>> Fs = 8000; t = [0:1/Fs:2]; % Two second records at a sampling rate of 8kHz
>> sound(m(theta,t,omega)/20,Fs); % Play (scaled) m(t) constructed using zero phases

Since the sound command clips magnitudes that exceed 1, the input vector is scaled by 1/20 to avoid clipping and the resulting sound distortion. The signals using other phase assignments are created and played in a similar fashion. How well does the human ear discern the differences in phase spectra? If you are like most people, you will not be able to discern any differences in how these waveforms sound.