- 1 Standing Wave Discrete Fourier Transform
- 2 Standing wave fast Fourier transform
- 3 Method and apparatus
- 4 Usage
- 5 Inverse standing wave discrete Fourier transform
- 6 Magnitude and phase
- 7 Variations spectral leakage and windowing
- 8 Multi-dimensional transforms
- 9 Standing wave discrete cosine transform
- 10 References
Standing Wave Discrete Fourier Transform
In applied mathematics, the standing wave discrete Fourier transform (SWDFT) is an apparatus and method to efficiently compute the frequency spectrum of a signal with identical magnitude to that of the discrete Fourier transform (DFT). The SWDFT is log2(N) times more efficient than a fast Fourier transform (FFT) when calculating N harmonics of a streaming data signal, such as bandlimited samples coming from an analog-to-digital converter (ADC). The standing wave fast Fourier transform (SWFFT), the standing wave discrete cosine transform (SWDCT), the standing wave fast cosine transform (SWFCT), and their inverses are also discussed.
The sequence of M complex numbers xm, ..., xM+m−1 is transformed into the sequence of N complex numbers X0[m], ..., XN−1[m] by the SWDFT according to the formula
where i is the imaginary unit and is a primitive M'th root of unity.
The inverse SWDFT is given by
If the sequence xm is a real bandlimited signal, then N=M/2 is the Nyquist component. An orthonormal transform can be obtained by setting N=M and using a scaling factor in both transform equations. For the SWDFT equation above Xk is the 1st spectral estimate of samples x0 to xM-1, the spectral estimate consisting of complex harmonics k=0 to N-1, Xk corresponds to SNUM=M in Figure 1. Likewise Xk is the 2nd spectral estimate of samples x1 to xM, Xk corresponds to SNUM=M+1 in Figure 1. Each spectral estimate Xk[m] of N complex harmonics k=0 to N-1 is calculated in Figure 1 in N operations, there is no computational latency, blocking or shuffling of data, and a specific number of harmonics can be calculated with no impact on computational efficiency. For the SWFFT the number of multiplies per spectral estimate Xk[m] can be reduced to an average of approximately log2(N). There is one spectral estimate calculated for each ADC sample, k is HCTR in Figure 1, and N and M are the same variable between the SWDFT equations and the method and apparatus described in Figure 1.
Signals and natural data usually have portions created by resonances amenable to periodic modeling, these periodic portions change in instantaneous frequency and amplitude over time. By continuously shifting the signal by one sample and performing spectral decomposition, the continuous stream of spectral components Xk[m] can be used for analysis, for example when the signal's periodic characteristics change, for example musical pitch. The intention of this analysis is to take a portion of the signal x0, ..., xM-1 and perform a basis function decomposition (eg harmonically related complex exponentials, orthogonal basis functions, or other model of the data being analyzed), then take a shifted version of the sequence and x1, ..., xM and perform a phase shifted basis function decomposition, etc. Analysis of the continuous stream of coefficients Xk[m] has application in data analysis of natural formations, video and audio processing, DNA sequencing, cell modeling, genetic research, bioengineering, spectral analysis, image and pattern analysis, etc.
The SWDFT uses a filter bank structure to continuously output spectral estimates with identical magnitude to that of a discrete Fourier transform (DFT) in N operations or less per spectral estimate. A new spectral estimate is computed for each ADC sample, and after initialization the SWDFT requires a factor of log2(N) fewer operations than an FFT. The SWDFT can be implemented in an FPGA, ASIC, ASSP or other hardware, or in software. A method to calculate N harmonics of a signal is shown in Figure 1.
// Standing-Wave Discrete Fourier Transform
// Program to calculate N harmonics of a signal using the
// previous M data samples
// Set N = M/2 to calculate all harmonics from DC up to the Nyquist component
// Bandlimited samples arrive from the Analog to Digital Converter (ADC):
Done = false // Set true when done
Har(N,M) = zeros // Initialize harmonic memory
Acc(N) = zeros // Initialize harmonic accumulators
Snum = 1 // Initialize spectrum counter
Ptr = 0 // Initialize circular memory pointer
j = (-1)^.5 // Math is complex
pi = 3.141592654 // Define pi
Print("Spectrum",Snum) // Current spectral estimate number
Data = ADC // Get next time-domain sample
for Hctr = 0 to N-1 // Loop for each harmonic - M/2 is the Nyquist component
Dataexp = Data * exp(-j*2pi*Hctr*Ptr/M) // Redundancy
Acc(Hctr) = Acc(Hctr)-Har(Hctr,Ptr) + Dataexp // Remove old & add new
Har(Hctr,Ptr) = Dataexp // Update harmonic memory
if Snum>=M then Print("Harmonic ",Hctr," = ",Acc(Hctr)/M)
Ptr = (Ptr + 1) modulo M // Increment circular memory pointer
Snum = Snum + 1 // Next spectral estimate
until Done = true // Loop forever
Figure 1. SWDFT - A method for Real-Time Spectral Analysis.
May be used freely for benevolent scientific research, education, and 503(c)(3) non-profit public benefit purposes if implemented in software using a serial method. Any other use including military or for-profit purposes requires a license.
While the FFT typically processes data in blocks of size N, requires N log2(N) operations per spectral estimate, and has a latency of log2(N), the SWDFT continuously processes data from an ADC, and computes a new spectral estimate for each newly arriving ADC sample in N operations. Unlike an FFT, the SWDFT requires no shuffling of blocks of data, has no interlace decomposition so harmonics are not reordered, and has minimal latency. For the SWDFT, an operation is defined as a complex accumulate, complex subtract, and a complex multiply.
The SWDFT algorithm requires an initialization of M × N operations, where M is the number of time-domain samples per SWDFT, N is the number of harmonics, with M > N. A memory buffer of size M × N complex data words is required. The ability to perform a DFT (with modified phase data) in N operations is accomplished by tuning each channel of a digital heterodyne filter bank to a harmonic of the DFT, and then replacing each of that channel's two low-pass filters with time-averages over the previous M "partially demodulated" samples. Using this method the real and imaginary parts of the harmonic for each filter bank channel are demodulated from the sampled time-domain signal.
Phase data is also computed; however, it is referenced to an absolute time origin. The key to the algorithm is the use of phase shifting harmonically related complex-exponential basis functions for each new spectral estimate. The SWDFT can return N estimates of the spectrum in the same number of operations as is required to calculate one DFT, an operation for the DFT being defined as a complex multiply and complex accumulate. The SWDFT was developed for the estimation of patient-sensitivity, baseline-pressure, and infusion-delay during surgery, and was first published in 1990 as part of the Master’s Thesis “Estimation of Patient-Sensitivity, Baseline-Pressure, and Infusion-Delay using Fast Spectral Estimates for Closed-Loop Regulation of Mean Arterial Pressure” at the University of California, San Diego..
A paper reviewing the closed-loop controller, not the actual method and apparatus used, whose solution led to development of the Standing-wave discrete Fourier transform and its derivations and inverses, was presented at the Engineering in Medicine and Biology Society, 1990., Proceedings of the Twelfth Annual International Conference of the IEEE.
After initialization, as each new time-domain data sample arrives, an estimate of the spectrum is calculated in N operations, based on the previous M time-domain data samples. If Nyquist's criterion is satisfied, the SWDFT after initialization will theoretically return a unique eigenvalue decomposition of the signal in N operations, assuming the previous M samples are properly updated into a circulant matrix. If the time domain samples are real, then only harmonics up to the Nyquist component (M/2) need be calculated in Figure 1.
Standing wave fast Fourier transform
The number of multiplications per operation can be reduced by identifying redundancy in the orthogonal basis functions (complex exponentials) of the SWDFT, and can be eliminated from the filter bank calculation. A description of multiplication redundancy for can be found in the book "Musical Applications of Microprocessors" by Hal Chamberlin. Systematically identifying and removing the redundancy from the SWDFT yields an even more efficient method of real-time harmonic spectral analysis, called the Standing Wave Fast Fourier Transform, or SWFFT.
While the FFT typically processes data in blocks of size N, requires Nlog2(N) operations per spectral estimate, and has a latency of log2(N), the SWFFT continuously processes data from an ADC, and computes a new spectral estimate for each newly arriving ADC sample in < N operations. Unlike an FFT, if only multiplication redundancy is eliminated, the SWFFT requires no shuffling of blocks of data, has no interlace decomposition so SWFFT harmonics are not reordered, and has minimal latency.
In the statement in Figure 1 with the comment "Redundancy", all cases for redundant multiples can be identified, and the term "Dataexp = Data * exp(-j*2*pi*Hctr*Ptr/M)" need only be calculated for the first occurrence. For example Hctr=3 and Ptr=7 will give the same Dataexp term as Hctr=7 and Ptr=3. M can be selected to maximize the redundancy, for example setting M = 1 × 2 × 3 × 4 × 5 = 120 requires an average of 14 complex multiplies per 60 component spectral estimate calculated for each new ADC sample. Even with the time averages replaced with traditional heterodyne filter bank low-pass filters, for example to minimize spectral leakage, the multiplication redundancy in the quadrature carriers of each harmonic still exists and can still be eliminated. With N = M = 256, a 256 point SWFFT requires an average of 20.344 complex multiples per spectral estimate, while the FFT requires (N/2)log2(N) = 1024. For this case, the SWFFT requires around a factor of 50 fewer multiples than the FFT, and there is a computational latency advantage of the SWFFT over the FFT of log2(N) = 8 clocks.
Method and apparatus
Let M equal the number of bandlimited time-domain data samples, as required by Nyquist's criterion, used to calculate the spectrum. The structure of the SWDFT is such that not all of the harmonics need be calculated for the algorithm to be efficiently utilized. For example, to compute seventh harmonic of a signal over the previous 65536 samples, set N = 1 and M = 65536. Initialization requires N × M = 65536 operations, and after initialization the filter bank channel requires only one (N = 1) operation per newly arriving data sample to calculate the 7th harmonic using a rectangular window applied to the previous 65536 data samples. For this example Hctr is held to 7 in Figure 1.
In order to calculate the DC component of a signal over the previous M time-domain data samples, for each data sample which arrives from the ADC, store it in a memory buffer of depth M, and add it to an accumulator. After M samples have arrived, the accumulator holds the DC component of the signal over the previous M samples. When the (M + 1)st sample arrives, lookup the memory buffer with address modulo M, subtract this memory buffer value from the accumulator, then store the new ADC sample in memory to replace to old value, and also add the new ADC sample the accumulator. The accumulator again holds the DC value of the previous M data samples, and this process continues for each newly arriving time-domain data sample. This is the moving average of the signal over the previous M data samples.
In a similar fashion, in order to calculate the Hth harmonic of a signal over the previous M band-limited time-domain data samples, for each data sample which arrives, multiply it times a complex-exponential of the form , store this product in a memory buffer (of size M), and then add it to an accumulator. After M samples have arrived, the accumulator contains the Hth harmonic of the signal using a rectangular window applied to the previous M data samples. When the (M + 1)st sample arrives, lookup the memory buffer with address modulo M, subtract this memory buffer value from the accumulator, store the new ADC sample times the complex-exponential in memory to replace the old value, and also add this new result to the accumulator. The accumulator again contains the Hth harmonic of the previous M data samples, and this process continues for each newly arriving time-domain data sample. This is a quadrature harmonically weighted moving average of the signal over the previous M data samples. Regardless of the phase of the complex-exponential, due to its periodicity, the quadrature components remain orthogonal. A program to calculate N harmonics of a signal is shown in Figure 1. For a real time-domain signal, the complex conjugate symmetry of the Fourier transform halves the number of operations required, only the 1st M/2 harmonics need be calculated.
The standing wave discrete Fourier transform (SWDFT) and variations, derivations, and inverses such as the ISWDFT, SWFFT, SWDCT, SWFCT, may be used freely for benevolent scientific research, education, and 503(c)(3) non-profit public benefit purposes if implemented in software using a serial method (not implemented in parallel for example in an ASIC, ASSP, FPGA, DSP chip, hardware etc). Any other use including military or for-profit purposes requires a license from the author of “Estimation of Patient-Sensitivity, Baseline-Pressure, and Infusion-Delay using Fast Spectral Estimates for Closed-Loop Regulation of Mean Arterial Pressure”, David Todd Johnson, published at the University of California, San Diego in 1990.
Inverse standing wave discrete Fourier transform
Changing the sign in the complex exponential exponent in Figure 1 calculates the inverse SWDFT. For the inverse SWDFT, the ADC is replaced with streaming frequency domain content, and Hctr is now a time domain sample estimate pointer. After initialization, each spectrum of M SWDFT harmonics are transformed into N time domain sample estimates in N operations. The inverse FFT requires N log2(N) operations. If the time-domain signal is real, then the number of operations in the inverse SWDFT is halved by using the symmetry properties of the inverse Fourier transform, reducing the loop to M/2 in Figure 1. The first half of the even part of the time-domain signal is contained in the real part of Acc(Hctr), and the first half of the odd part of the time-domain signal is contained in the imaginary part of Acc(Hctr). Normalization factors such as are the same between the SWDFT and DFT for both forward and inverse transforms, they are omitted in Figure 1 for computational efficiency. The inverse SWDFT provides maximal overlapping of time-domain signal estimates, while the SWDFT provides maximal overlapping of frequency-domain spectral estimates.
Magnitude and phase
The magnitude is identical to performing an DFT over the previous M data samples. The phase will appear referenced to a single time origin. Magnitude and phase for each spectral component Hctr are given by the following equations:
For example, if the signal is many cycles of a harmonic sine wave, the SWDFT's corresponding harmonic will have a constant magnitude and phase for all spectral calculations. An FFT over the previous M samples for each newly arriving time-domain sample will result in a different phase angle for each successive spectral estimate. For the SWDFT, each harmonic’s phase appears as a constant, referenced to the quadrature components of a complex exponential beginning at the time origin of when the algorithm is started. The additional phase shift term will exactly cancel when the inverse SWDFT is performed, and does not affect the magnitude of the frequency-domain spectral estimate. If one multiplies an SWDFT harmonic by the complex exponential , then the phase and magnitude of the DFT or FFT and the SWDFT are identical.
Variations spectral leakage and windowing
Variations include methods for shaping the spectrum and controlling spectral leakage using windowing via spectral convolution (see Harris, Fredric j. (January 1978). ""On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform"". Proceedings of the IEEE 66 (1): 51–83. doi:10.1109/PROC.1978.10837.), and sweeping, tuning or real-time adapting each heterodyne quadrature carrier waveform. For example, by adjusting M for each Hctr value in the statement in Figure 1 with “Dataexp = Data * exp(-j*2*pi*Hctr*Ptr/M[Hctr])”, additional frequency resolution can be achieved, because each heterodyne can be precisely tuned with an arbitrarily accurate resolution, making the signal periodic within the sampling window. Modifying the moving time averages back to traditional heterodyne low-pass filters should be considered over time domain windowing, as time domain windowing requires the mathematical equivalent of an additional window multiply for each spectral estimate and additional quadrature multiplies. Windowing via spectral convolution using the Hanning window is efficient due to using binary shifts instead of multiplies. Methods to account for the additional phase shift term of the SWDFT should be considered.
Taking the example above of tuning each quadrature carrier even further, one could replace each quadrature heterodyne carrier with a quadrature arbitrary function generator, for example replacing the sine carrier with a sum of independent amplitude harmonically related sine waves, and the cosine carrier with a sum of independent amplitude harmonically related cosine waves. Each arbitrary waveform quadrature modulator could be precisely tuned to the particular frequencies relevant to the signal being analyzed. As the basis function decomposition methods become more exotic allowing for better analysis and processing of certain signal types, mathematical properties such basis functional orthogonality and unique signal decomposition may be lost.
Multi-dimensional transforms are obtained by cascading SWDFTs, one cascade per dimension. Harmonics from one SWDFT feed the sample data input of the next SWDFT. Additional harmonic memory Har(Hctr,Ptr) in Figure 1 is required for each additional cascade. If all additional N spectral estimates are used in the cascade, the number of spectral estimates grows by a factor N for each new dimension.
Hybrids transforms are possible. For example, a 2-dimensional FFT could be streamed into an SWDFT in the time dimension. This transform has maximal overlapping of spectral estimates in the time dimension only.
Standing wave discrete cosine transform
The SWDFT method can be applied to other transforms, for example a generalization of the popular type-II form of the discrete cosine transform (DCT-II) with M samples and N components is given by:
The sequence of M numbers xm, ..., xM+m−1 is transformed into the sequence of N numbers X0[m], ..., XN−1[m] by the SWDCT-II according to the formula
The inverse SWDCT-II is given by
If the sequence xm is a real bandlimited signal, then N=M/2 is the Nyquist component. An orthonormal transform can be obtained by setting N=M and using a scaling factor in both transform equations.
For the DCT-II form of the standing wave discrete cosine transform (SWDCT-II), replace the statement with the comment "Redundancy" in Figure 1 with
- Dataexp = Data * cos((pi*Hctr*(Ptr + 0.5))/M).
A standing wave fast cosine transform is accomplished by identify redundancy in the "Redundancy" statement in Figure 1 and removing it from the filterbank calculation. The inverse standing wave discrete cosine transform is obtained in Figure 1 by instead using the negative cosine and replacing ADC samples with streaming spectral samples calculated from the forward transform.
A two-dimensional SWDCT is created by cascading two SWDCTs. A steaming two-dimensional SWDCT encoder has the advantage of maximal overlapping of spectral estimates. If used for image compression, the additional spectral estimates can be used to minimize compression artifacts. The SWDCT can streamed in the time dimension as well, additional spectral estimates are calculated.
- Johnson, David T. Estimation of Patient-Sensitivity, Baseline-Pressure, and Infusion-Delay using Fast Spectral Estimates for Closed-Loop Regulation of Mean Arterial Pressure. University of California, San Diego, S&E Books, Call Number TK3.7 .J64 1990
- Identification Of Patient-sensitivity, Baseline-pressure, And Infusion-delay Using Fast Spectral Estimates During Closed-loop Control Of Mean Arterial Pressure. Johnson, D.T.; Schneider, A.M.; Martin, J.F.; Quinn, M.L. University of California This paper appears in: Engineering in Medicine and Biology Society, 1990., Proceedings of the Twelfth Annual International Conference of the IEEE Issue Date: 1-4 Nov 1990 On page(s): 1878 - 1879 Print ISBN: 0-87942-559-8 Digital Object Identifier: 10.1109/IEMBS.1990.692063 Date of Current Version: 06 August 2002
- Chamberlin, Hal Musical Applications of Microprocessors. Hayden Books, Nov 1985