Copyright © 1999, 2002, 2007 Kees Krijnen.

The setup of a simple dynamic laser light scattering DSP system (fig.1) consist of:

- A conversion of a light detector analog signal into digital information, in the form of a sequence of binary numbers. This involves two operations: sampling and analog-to-digital (A/D) conversion.
- Perform numerical operations on the digital information by special-purpose digital hardware or PC. In this case the operations are spectrum analysis, "analog" autocorrelation, or both simultaneously.
- Display DSP process on PC monitor continuously.

*Fig. 1*

To facilitate digital processing, a continuous-time signal must be
converted to a sequence of numbers [REF1].
Sampling is essentially a selection of a finite number of data at any finite
time interval as representatives of the infinite amount of data contained in
the continuous-time signal in that interval. Sampling is the fundamental
operation of digital signal processing, and avoiding (or at least minimizing)
*aliasing* is the most important aspect of sampling. Aliasing is the
fundamental distortion introduced by sampling. Aliasing causes the digitized
pattern *not* look like the analog signal digitized. If analog signals
at higher frequencies are present, the sampling process will produce sum and
difference frequencies, in particular, will produce spurious low-frequency
signals - in the signal passband - that cannot be distinguished from the
signal. The *Nyquist* criterion provides a correct *reconstruction*
of the analog signal, stated by *Shannon's* reconstruction theorem.
Reconstruction is the operation of converting a sampled signal back to a
continuous-time signal. Aliasing is caused by sampling at a lower rate than
the Nyquist rate for the given signal.

Nyquist rate: sampling frequency ≥ 2 times highest frequency
signal or *noise* component

Analog noise in data-acquisition systems take three basic forms [REF5]: *transmitted noise*, inherent in
the received signal, *device noise*, generated within the devices used
in data acquisition (preamps, converters, etc.), and *induced noise*,
`picked up' from the outside world, power supplies, logic, or other analog
channels, by magnetic, electrostatic, or galvanic coupling.

Noise is either *random* or *coherent* (i.e., related to some
noise-inducing phenomenon within or outside the system). Random noise is
usually generated within components, such as resistors, semiconductor
junctions, or transformer cores, while coherent noise is either locally
generated by processes, such as modulation/demodulation, or coupled-in.
Coherent noise often takes the form of `spikes', although it may be of any
shape.

The finite resolution of the A/D conversion process introduces `quantization noise', which may be thought of as either a truncation (or roundoff) error, or as a `digitizing' noise, depending on the context.

Noise is characterized in terms of either *root-mean-square* (rms) or
*peak-to-peak* measurements, within a stated bandwidth. Noise leads to
random errors and missing codes in an ADC. Since these errors can lead to a
closed-loop system making inappropriate responses to specific false data
points, high-resolution systems, where random noise tends to limit resolution
and to cause errors, can benefit by digital or analog filtering to smooth the
data.

If higher frequencies signal or noise than the Nyquist rate are part of
the analog signal, a low-pass filter *must* be applied to match the
Nyquist criterion.

The departure point for a discussion of spectral analysis is the Fourier transform equation pair [REF5]:

Fourier Transform

*Inverse* Fourier Transform

where t is time, ω is angular frequency (2πf), x(t) is the signal - a function of time - and X(ω) is its counterpart in the frequency (spectral) domain. These equations give us, at least formally, the mechanics for taking a signal' time-domain representation and resolving it into its spectral weights - called Fourier coefficients. The Fourier series is a trigonometric series F(f) by which we may approximate some arbitrary function f(t) or x(t) - real world signal [REF3]. Specifically, F(f) is the series (another way to describe above frequency-domain function X(ω)):

F(t) = A_{0} + A_{1}Sin(t)+
B_{1}Cos(t) + A_{2}Sin(2t) + B_{2}Cos(2t) + ... +
A_{n}Sin(nt) + B_{n}Cos(nt)

e^{iωt} = Cos(ωt) +
iSin(ωt)

In the limit as n (i.e. the number of terms) approaches infinity, the
trigonometric series is equal to the real world signal. It means that any
real world signal can be considered as a composite waveshape which can be
breaked down into its component sinusoids. A0, A1, B1, ...and A_{n},
B_{n}, are the Fourier coefficients. When the Fourier coefficients
are known the time-domain function x(t) can be easily calculated -
reconstructed - by summation of the component sinusoids. Since the Fourier
equations require continuous integrals, they have only indirect bearing on
digital processors. However, under certain circumstances, a sampled
(digitized) signal can be related faithfully to its Fourier coefficients
through the discrete Fourier transform (DFT):

(*Forward*) DFT

*Inverse* DFT

where: | F(k) = frequency components or transform f(n) = time base data points or inverse transform N = number of data points n = discrete times, 0 to N-1 k = discrete frequencies, 0 to N-1 W _{N} = e^{-i2π/N} = Cos(2π/N) -
iSin(2π/N) |

Provided that the signal is sampled frequently enough (Nyquist rate), and assuming that the signal is periodic, the above DFT equations hold exactly. What is most interesting from the standpoint of DSP is that the DFT equation provides us with a means to estimate spectral content by digitizing an incoming signal and simply performing a series of multiply/accumulate operations.

A process is required that can isolate, or `detect' any single
*harmonic* (sine wave) component within a complex waveshape [REF3]. To be useful quantitatively, it will
have to be capable of measuring the amplitude and phase of each harmonic
component. While measuring any harmonic component of a composite wave, all of
the other harmonic components must be ignored.

According to Fourier the way to detect and measure a sine wave component
within a complex waveshape is to multiply through by a unit amplitude sine
wave of the identical frequency and then find the *average value* of
the resultant. The two functions must have the same number of data points,
corresponding domains, etc. The reason this works is that both sine waves are
both positive at the same time, and both negative at the same time, so that
the products of these two functions will have a positive average value. Since
the average value of a sinusoid is normally zero, the original digitized sine
wave has been `rectified' or `detected'. This scheme will work for any
harmonic component, changing the frequency of the *unit amplitude
sine* wave to match the frequency of the harmonic component being
detected. The same approach will work for the cosine components if the unit
sine function is replaced with a unit cosine function.

Orthogonality is the method used to extract the sinusoid being analyzed while ignoring all other sinusoidal components [REF4]. Orthogonality implies that if one sinusoidal is multiplied by another, an average value of zero will obtained (averaged over any number of full cycles), if either:

- the sinusoids are dephased by 90° or,
- the sinusoids are integer multiple frequencies.

It is necessary to deal with *full cycles* of a sinusoid for the
average value to be zero. This requirement can be satisfied by scaling the
frequency term so that exactly one full sinusoid fits into the domain of the
function working with.

The composite waveform is generated by summing in harmonic components:

F(t) = A_{0} + A_{1}Sin(t)+ B_{1}Cos(t) +
A_{2}Sin(2t) + B_{2}Cos(2t) + ... +
A_{n}Sin(nt) + B_{n}Cos(nt) |

If this composite function is multiplied by Sin(Kt) (or, alternatively, Cos(Kt)), where K is an integer, the following terms on the right hand side of the equation are created:

A_{0}Sin(Kt) + A_{1}Sin(t)Sin(Kt) +
B_{1}Cos(t)Sin(Kt) + ... + A_{k}Sin2(Kt) +
B_{k}Cos(Kt)Sin(Kt) + ... + A_{n}Sin(nt)Sin(Kt) +
B_{n}Cos(nt)Sin(Kt) |

Treating each of these terms as individual functions, if argument (Kt) equals the argument of the sinusoid it multiplies, that component will be `rectified'. Otherwise, the component will not be rectified. From what is shown above, two sinusoids of ± 90° phase relationship (i.e. Sine/Cosine), or integer multiple frequency, represent orthogonal functions. As such, when summed over all values within the domain of definition, they will all yield a zero resultant (regardless of whether they are handled as individual terms of combined into a composite waveform). This is precisely what is demanded from a procedure to isolate the harmonic components of an arbitrary waveform.

The DFT algorithm consist of stepping trough the digitized data points of the input time-domain function, multiplying each point by sine and cosine functions, and summing the resulting products into accumulators (one for the sine component and one for the cosine). When every point has been processed in this manner, the accumulators (i.e. the sum-totals of the preceding process) is divided by the number of data points. The resulting quantities are the average values for the sine and cosine components at the frequency being investigated as described above. This process must be repeated for all integer multiple frequencies up to the frequency that is equal to the sampling rate minus 1 (i.e. twice the Nyquist frequency minus 1).

If the digitized matrix f(t) is considered to be a row matrix, the process can be described as a matrix multiplication [REF6,7]:

W_{N} = e^{-i2π/N} =
Cos(2π/N) - iSin(2π/N)

The DFT W_{NN} square coefficient matrix is called the array of
*Twiddle Factors.* One is needed for the Sine function and one for the
Cosine function.

Unfortunately, the large number of multiplications required by the DFT
limits its use in real-time signal processing. The computational complexity
of the DFT grows with the square of the number of input points; to resolve a
signal of length N into N spectral components F(k) requires N^{2}
*complex *multiplications (4N^{2} *real* multiplications).
Given the large number of input points needed to provide acceptable spectral
resolution, the computational requirements of the DFT are prohibitive for
most applications. For this reason, tremendous effort was devoted to
developing more efficient ways of computing the DFT. This effort produced the
fast Fourier transform (FFT) invented by *Cooley and Tukey* (1965). The
FFT uses mathematical shortcuts to reduce the number of calculations the DFT
requires.

The DFT (and FFT) requirement to deal with full cycles of a sinusoid rises inevitably the need for a high-pass filter. The frame length and timing of data acquisition determines the low cut-off frequency. The Nyquist criterion requires a low-pass filter, under experimental conditions this means a band-pass filter has to be applied to the input signal.

Suppose an N-point DFT is accomplished by performing two N/2 -point DFTs
and combining the outputs of the smaller DFTs to give the same output as the
original DFT. The original DFT requires N^{2} complex multiplications
and N^{2} complex additions. Each DFT of N/2 input samples requires
(N/2)^{2} = N^{2}/4 multiplications and additions, a total of
N^{2}/2 calculations for the complete DFT. Dividing the DFT into two
smaller DFTs reduces the number of computations by 50 percent. Each of these
smaller DFTs can be divided in half, yielding N/4 -point DFTs. If we continue
dividing the N-point DFT calculation into smaller DFTs until we have only
two-point DFTs, the number of complex multiplications and additions is
reduced to Nlog_{2}N. For example, a 1024-point DFT requires over a
million complex additions and multiplications. A 1024-DFT divided down into
two-point DFTs needs fewer than ten thousand complex additions and
multiplications, a reduction of over 99 percent!

Dividing the DFT into smaller DFTs is the basis of the FFT. A radix-2 FFT divides the FFT into two smaller DFTs, each of which is divided into two smaller DFTs, and so on, resulting in a combination of two-point DFTs. The radix-2 Decimation-In-Time (DIT) FFT divides the input time-domain sequence into two groups, one of even samples and the other of odd samples. N/2-point DFTs are performed on these sub-sequences, and their outputs are combined to form the N-point DFT.

The input time-domain sequence in equation (1), is divided into even and odd sub-sequences to become the equation:

(2)

By substitution equation (2) becomes:

(3)

By substitution equation (3) can also be expressed as two equations:

(4.1)

(4.2)

Together these equations (4.1) and (4.2) form an N-point FFT. Each of
these two equations can be divided to form two more. The final decimation
occurs when each pair of equations together computes a two-point DFT
(*one* point per equation). The pair of equations that make up the
two-point DFT is called a radix-2 `butterfly' (fig. 2). The butterfly is the
core calculation of the FFT. The entire FFT is performed by combining
butterflies in patterns determined by the FFT algorithm based on the
*Similarity*, *Shifting*, *Addition* and *Stretching*
transform theorems. An array of 2^{n} data points yields *n stages
of computation* in FFT.

*fig. 2*

In general, the input of FFTs are assumed to be complex. When input is purely real - for example, in case of digitized analog data - their symmetric properties compute the DFT very efficiently [REF9].

One such optimized real FFT algorithm is the packing algorithm. The original 2N-point real input sequence is packed into an N-point complex sequence. Next, an N-point FFT is performed on the complex sequence. Finally the resulting N-point complex output is unpacked into the 2N-point complex sequence, which corresponds to the DFT of the original 2N-point real input.

Using this strategy, the FFT size can be reduced by half, at the FFT cost function of N operations to pack the input and unpack the output. Thus, the real FFT algorithm computes the DFT of a real input sequence almost twice as fast as the general FFT algorithm.

The TMS320C542 (see *PALS Experimental Setup*
chapter) real FFT algorithm is a radix-2, DIT DFT algorithm in four
phases:

In phase 1 - *packing and bit-reversal of real input* - the input is
bit-reversed so that the output at the end of the entire algorithm is in
natural order. First, the original 2N-point real input sequence is copied
into contiguous sections of memory and interpreted as an N-point complex
sequence. The even-indexed real inputs form the real part of the complex
sequence and the odd-indexed real inputs form the imaginary part. This
process is called packing. Next, this complex sequence is bit-reversed and
stored into the data processing buffer.

In phase 2 - *N-point complex FFT* - an N-point complex DIT FFT is
performed in place in the data-processing buffer. The twiddle factors are in
Q15 format and are stored in two separate tables, pointed to by sine and
cosine. Each table contains 512 values, corresponding to angles ranging from
0 to almost 180 degrees. The indexing scheme used in this algorithm permits
the same twiddle tables for inputs of different sizes.

Phase 3 - *separation of odd and even parts* - separates the FFT
output to compute four independent sequences, which are the even real, odd
real, even imaginary, and the odd imaginary parts.

Phase 4 - *generation of final output *- performs one more set of
butterflies to generate the 2N-point complex output, which corresponds to the
DFT of the original 2N-point real input sequence.

Often, the last stage in a FFT spectral analysis is the computation of the
*power spectrum* from the frequency-domain transform of the input
time-domain sequence.

The power spectrum is a meaningful way of representing random signals in
the frequency domain by their *mean-square value* or *energy
content*. Now from Parseval's theorem [REF10]:

(5)

where ω = 2πν and ν is expressed in Hertz.
Equation (5) states that the energy content of f(t) is given by 1/(2π)
times the area under the |F(ω)|^{2} curve. For this reason, the
quantity |F(ω)|^{2} is called the power spectrum or energy
spectral density of f(t).