Given an input array f(n) = a + ib, an output array g(n) = c + id, and a complex constant (entered from the keyboard) r + is, the available operations are:

T, I      Forward/Inverse Fast Fourier Transform (FFT)

The need to demonstrate the principles and characteristics of the Fourier transform provided the original incentive to develop SIGNALS. The convention for the indices is selected in the MAIN MENU (option "I") to either [-N/2, N/2 - 1] or [0, N - 1]. The normalization factor A‹ in the forward transform also is selected in the MAIN MENU (option "M"). The allowed values for A are 1/N, 1, or SQR(1/N), with the respective normalization factors in the inverse transform of 1, 1/N, and SQR(1/N). The time required for the computation is displayed, to allow comparison with the discrete Fourier transform (see next item).

K      Discrete Fourier Transform (DFT)

This operation is included for demonstrations only. Because the DFT may be computed over "any" number of frequencies within a selected range, the number of samples N in the array would have to be changed whenever the DFT is computed. This would require deleting the other array. The DFT may be computed even at frequencies above/below the Nyquist frequency to demonstrate its periodicity. If computing the DFT of a 256-sample array, select 257 samples in the DFT to produce a set of output frequencies "symmetric" about 0. The values of the displayed transform may be examined by using the VALUES option in the DFT display menu.

C, B     Discrete Cosine Transform (DCT)

This is the fundamental operation for a number of signal/image compression algorithms. The DCT of the entire N-pixel array is computed while assuming that the array is indexed from [0,N-1]. The indices of the array are changed to make this so. The DCT of any real-valued array is real, and the DCT of a complex-valued array is complex. The DCT is computed by FFT on a "symmetric" array of size 2N (with the appropriate phase factor to force symmetry).

H            Hartley Transform

This is a variant of the Fourier transform which sums the cosine and sine functions without the imaginary weight factor. It has the arguable advantage of yielding a real-valued transform for a real-valued input. The transforms of the even and odd real parts of the input are summed, rather than existing in different arrays as is true for the Fourier transform.

This operation, used in image compression, may be viewed as a transform using a bitonal function (values +/- 1). It is most often applied to arrays indexed over the range [0, N - 1], but the generalization to indices [-N/2, N/2 - 1] also is valid.

: ;      Forward/Inverse Complex Cepstrum

The cepstrum (anagram of "spectrum") is the cascade of the operation of the Fourier transform, complex logarithm, and inverse Fourier transform, and is used in "homomorphic" filtering to segment functions that have been multiplied. The logarithm function uses double-precision arithmetic. The inverse cepstrum is the cascade of the inverse Fourier transform, complex exponentiation, and the forward Fourier transform.

R         Convert a Magnitude/Phase array to Real/Imaginary Parts

M        Convert a Real/Imaginary array to Magnitude/Phase

A simple phase unwrapping routine may be used which derives a continuous range of output phases for "well-behaved" functions whose phase varies slowly enough. The phase-unwrapping routine fails where the function is aliased, and may be used to illustrate those locations.

-         Threshold at selected amplitude

The thresholding and clipping operations (below) may be used to demonstrate the effects of nonlinear operations. The range of the output of the threshold may be selected as [0, 1] or [-1,+1]; the latter range is useful when characterizing the spectrum of the operation, as the mean of the output may be approximately zero.

|      Clip at selected amplitude

A less severe distortion than thresholding, the clipping operation retains values of the input function within the range, and may be used to model saturation effects, as in a signal amplifier.

\      Quantize real and imaginary parts of f(n) separately

SIGNALS includes two types of quantization of complex numbers. This is the simpler, where the real and imaginary parts are quantized individually using the same 1-D quantizer. The values of the individual functions may be quantized independently (the quantization error at each sample is discarded), or adaptively by error diffusion (quantization error is subtracted from the next sample to be quantized). This allows the resulting functions to be compared and the total quantization error to be computed.

In the "\" option, several choices are available for quantization:

• R - quantize by rounding to nearest integer (CINT function in BASIC)
• I - truncate to greatest integer, also called two's complement truncation (INT function in BASIC)
• F - truncate absolute value to greatest integer, also called one's complement truncation (FIX function in BASIC)

Other options allow the user to specify the quantization parameters.

M - User selects the maximum and minimum OUTPUT quantization states gmax and gmin and the number of output levels Nlevels. The other N-2 reconstruction levels are computed. The user also chooses whether the quantization error is discarded or diffused, and whether the output is scaled to (usually bipolar) amplitude or displayed as the nonnegative quantum number of the state.

Q - User selects the maximum and minimum INPUT amplitudes fmax and fmin for quantization and the number of output levels Nlevels. The (Nlevels) reconstruction levels gn are computed from the parameters and are located midway between the threshold levels; the maximum and minimum output levels gmaxand gmin have smaller magnitudes than the specified fmax and fmin. This is the scheme described by Pratt in Digital Image Processing, 2nd Edition, p.160.

For example, if quantizing a signal with maximum amplitude of +1 and minimum amplitude of -1 to two levels, both option M and Q place the threshold at 0, but option M places the output levels at +1 and -1, while option Q places them midway between the input levels and the threshold, i.e., at +0.5 and -0.5.

~      Quantize complex amplitudes of f(n) as 2-D vectors.

User selects output states, type of quantization (independent or error-diffused). These operations quantizes the 2-D complex amplitude as a 2-D vector to sets of 2-D quantization states on the complex plane. Again, two choices are available for dealing with the resulting quantization errors at each sample: the errors may be discarded (independent quantization) or propagated to the next sample (error-diffused quantization). This option was included for studies of the effects of quantization of the Fourier transform on the original array, as in Fourier-transform computer-generated holography.

=       Equalize the histogram of f(n)

Histogram equalization is an array-dependent nonlinear lookup table which maps the amplitudes to a set of levels with approximately equal populations.

Q      Compute Squared Magnitude of the array

Computes the squared magnitude of the data, also called the power or intensity, and is used when modeling optical systems to illustrate the effect of square-law detection.

+      Add complex-valued constant to array:            g(n) = f(n) + (r + is)

The user specifies the real and imaginary parts r and s to be added to the amplitude of all samples.

x      Multiply array by complex-valued constant:     g(n) = f(n) x (r + is)

The user specifies the real and imaginary parts r and s of the complex number to be multiplied by the amplitude of all samples.

^      Raise complex amplitude to complex-valued power:   g(n) = (f(n))(r + is) = e(r+is) LOG(f(n)) = e(r+is) LOG(a+ib)

This operation implements the generalization of de Moivre's theorem in complex analysis, and is performed using double-precision arithmetic. Note that the components of the power may be noninteger, and that the output is generally complex, even if the amplitude is real (i.e., if b = 0). The process is implemented as shown, by raising e to the power computed by multiplying the logarithm of the base by the exponent. The unwrapped phase of the logarithm is used, which may create data overflow when using complex-valued exponents; the maximum allowed value for the real part of the product of the exponent and the logarithm is only +88. This is rarely a problem in most applications, and the amplitudes are tested before performing the calculation. If an overflow is detected, a warning message offers the choice of aborting or continuing the process. The latter case results in "clipped" amplitudes. However, the simplicity of the phase unwrapping algorithm makes it prone to errors when applied to "complicated" functions which are not smooth. Users are warned to use any results from this operation with discretion.

L      Logarithm

Several options are available, including the complex logarithm of the magnitude (suboption L:   g(n) = LOG|f(n)| + 0i = LOG|a + ib|), where the output is the natural logarithm of the magnitude of the complex number; the resulting imaginary part is 0. The output is set to -88 when a = b = 0, which is the limit for single-precision arithmetic. The complete complex logarithm evaluates g(n) = LOG|a + ib| + i TAN-1(b/a), where the imaginary part of the output is the unwrapped phase of the array values. The result is the complex logarithm; the imaginary part of the output is the phase of the array after unwrapping. The magnitude of the logarithm is set to -88 when a = b = 0.

E      Complex Exponential:     g(n) = e(f(n)) = e(a + ib) = ea COS(b) + i ea SIN(b)

Computes the complex power of e = 2.71828..., and may be used to illustrate the Euler relation. The computation is performed using double-precision arithmetic and the result is truncated to single-precision numbers. Again, data arrays are tested for overflow before completing the calculation. A warning message offers the choice of aborting the calculation or continuing and saving "clipped" data.

N      Complex Inverse (reciprocal): g(n) = (a + ib)-1 = (a/(a + ib)2) - i (b/(a + ib)2)

This operation is performed using double-precision numbers and truncated to single precision after testing for overflow. The output is set to zero if a = b = 0.

'      Differentiation

Computes the difference of adjacent amplitudes, except at n = -N/2 (or n = 0 if the index range is [0, N - 1]): g[-N/2] = f[-N/2]. This is tantamount to a boundary condition and ensures that the cascade of the derivative and the integral (below) will return the original function.

\$       Integration:

At each sample, the sum of all samples at all preceding indices is computed.

>      Circularly Translate Array

Allows the array to be translated in a periodic sense, so that values shifted "off the edge" on one side of the array reappear at the other side.

V      "Time"-Reversal

This "flips" the array left to right about the origin. Since the array has an even number of elements, the sample n = -N/2 is not moved. This is another manifestation of the assumed periodicity of the discrete array.

J      Complex Conjugate: g(n) = f*(n) = (a - i b)

Multiplies imaginary part by -1.

&      Autocorrelation (computed via FFTs): g(n) = f(n) & f(n) = SUM( f(n' - n) f*(n') ) over n'

The complex autocorrelation of the array, which is guaranteed to yield a Hermitian result (even real part and odd imaginary part) with peak magnitude at the origin.

*      Convolution with a small kernel h(x): g(n) = f(n) * h(n) = SUM( f(n') h(n - n') ) over n'

Computes the convolution with small kernels of odd size. When the kernel is "off the edge" of the array, the input values are assumed to be zero. A small library of precomputed kernels is available (uniform average, derivatives, etc.), or the user may enter specific values. The time required for the computation is displayed, thus allowing comparison with the same operation via FFTs.

Z      Normalize array: g(n) =   f(n) / f(0)    or   f(n) / fmax

Rescales the array amplitudes so that a specific value is unity. Available choices for the scale factor include the amplitude of the real part and of the imaginary part at the origin, the maximum absolute values of the amplitudes, or a value entered by the user.

!      "Checkerboard" (multiply odd-index pixels by -1): g(n) = f(n) x (-1)n

Used to illustrate the steps in computing a "centered-array" FFT using the operation over the range [0, N - 1].

S      Swap Real/Imaginary Parts of the array

The real and imaginary parts are exchanged; equivalent to multiplying the real part by i and the imaginary part by -i.

)      Nonlinear Filtering over a neighborhood (Median, Maximum, Minimum, Variance)

The most common nonlinear operations may be computed over odd-sized neighborhoods. The array is assumed to "wrap around", so that neighborhood is always the specified size. These operations are used to illustrate the effects of nonlinear processes on the spectrum of the signal.

|      Compute Absolute Value of Real Array (only in Modulation Menu for Real Array)

This (rarely used, even by me) operation is intended to simplify the generation of arrays as magnitude/phase. After generating the funcitons forming the magnitude array, it computes the absolute value to ensure that all values are positive or zero before generating the phase array.