The above Goertzel algorithm computes a single complex DFT spectral bin value for every N input time samples. Here we describe a sliding DFT process whose spectral bin output rate is equal to the input data rate, on a samplebysample basis, with the advantage that it requires fewer computations than the Goertzel algorithm for realtime spectral analysis. In applications where a new DFT output spectrum is desired every sample, or every few samples, the sliding DFT is computationally simpler than the traditional radix2 FFT.
Table 134. SingleBin DFT Computational Comparisons
Method 
Real multiplies 
Real additions 

Singlebin DFT 
4N 
2N 
FFT 
2Nlog2N 
Nlog2N 
Goertzel 
N + 2 
2N + 1 
13.18.1 The Sliding DFT Algorithm
The sliding DFT (SDFT) algorithm computes a single bin result of an Npoint DFT on time samples within a slidingwindow. That is, for the mth bin of an Npoint DFT, the SDFT computes
Let's take care to understand the notation of Xm(q). Typically, as in Chapter 3, the index of a DFT result value was the frequencydomain index m. In Eq. (1385) the index of the DFT result is a timedomain index q = 0, 1, 2, 3, ... , such that our first mthbin SDFT is Xm(0), our second SDFT is Xm(1), and so on.
An example SDFT analysis time window is shown in Figure 1346(a) where Xm(0) is computed for the N = 16 time samples x(0) to x(15). The time window is then advanced one sample, as in Figure 1346(b), and the new Xm(1) is calculated. The value of this process is that each new DFT result is efficiently computed directly from the result of the previous DFT. The incremental advance of the time window for each output computation leads to the name sliding DFT or slidingwindow DFT.
Figure 1346. Analysis window for two 16point DFTs: (a) data samples in the first computation; (b) second computation samples.
We can develop the mathematical expression for the SDFT as follows: the standard Npoint DFT equation, of the mth DFT bin, for the qth DFT of the time sequence x(q), x(q+1), ..., x(q+N–1) is
Equation 1386
(Variable m is the frequencydomain index, where m = 0, 1, 2, ..., N–1.) Likewise, the expression for the next DFT, the (q+1)th DFT performed on time samples x(q+1), x(q+2), ..., x(q+N) , is
Equation 1387
Letting p = n+1 in Eq. (1387), we can write
Equation 1388
Shifting the limits of summation in Eq. (1388), and including the appropriate terms (subtract the p = 0 term and add the p = N term) to compensate for the shifted limits, we write
Equation 1389
Factoring the common exponential term (ej2pm/N), we write
Equation 1390
Recognizing the summation in the brackets being equal to the previous Xm(q) in Eq. (1386), and e–j2pm = 1, we write the desired recursive expression for the sliding Npoint DFT as
Equation 1391
where Xm(q+1) is the new singlebin DFT result and Xm(q) is the previous singlebin DFT value. The superscript m reminds us that the Xm(q) spectral samples are those associated with the mth DFT bin.
Let's plug some numbers into Eq. (1391) to reveal the nature of its time indexing. If N = 20, then 20 time samples (x(0) to x(19)) are needed to compute the first result Xm(0). The computation of Xm(1) is then
Equation 1392
Due to our derivation method's time indexing, Eq. (1392) appears compelled to look into the future for x(20) to compute Xm(1). With no loss in generality, we can modify Eq. (1391)'s time indexing so that the x(n) input samples and the Xm(q) output samples use the same time index n. That modification yields our SDFT timedomain difference equation of
Equation 1393
Equation (1393) reveals the value of this process in computing realtime spectra. We compute Xm(n) by subtracting the x(n–N) sample and adding the current x(n) sample to the previous Xm(n–1), and phase shifting the result. Thus the SDFT requires only two real additions and one complex multiply per output sample. Not bad at all! Equation (1393) leads to the singlebin SDFT filter implementation shown in Figure 1347.
Figure 1347. Singlebin sliding DFT filter structure.
The singlebin SDFT algorithm is implemented as an IIR filter with a comb filter followed by a complex resonator. (If you need to compute all N DFT spectral components, N resonators with m = 0 to N–1 will be needed, all driven by a single comb filter.) The comb filter delay of N samples forces the SDFT filter's transient response to be N samples in length, so the output will not reach steady state until the Xm(N–1) sample. The output will not be valid, or equivalent to Eq. (1386)'s Xm(q), until N input samples have been processed. The ztransform of Eq. (1393) is
Equation 1394
where factors of Xm(z) and X(z) are collected yielding the zdomain transfer function for the mth bin of the SDFT filter as
Equation 1395
This complex filter has N zeros equally spaced around the zdomain's unit circle, due to the Ndelay comb filter, as well as a single pole canceling the zero at z = ej2pm/N. The SDFT filter's complex unit impulse response h(n) and pole/zero locations are shown in Figure 1348 for the example where m = 2 and N = 20.
Figure 1348. Sliding DFT characteristics for m = 2 and N = 20: (a) complex impulse response; (b) pole/zero locations.
Because of the comb subfilter, the SDFT filter's complex sinusoidal unit impulse response is finite in length—truncated in time to N samples—and that property makes the frequency magnitude response of the SDFT filter identical to the sin(Nx)/sin(x) response of a single DFT bin centered at a frequency of 2pm/N.
One of the attributes of the SDFT is that once an Xm(n) is obtained, the number of computations to compute Xm(n+1) is fixed and independent of N. A computational workload comparison between the Goertzel and SDFT filters is provided later in this section. Unlike the radix2 FFT, the SDFT's N can be any positive integer giving us greater flexibility to tune the SDFT's center frequency by defining integer m such that m = Nfi/fs, when fi is a frequency of interest in Hz and fs is the signal sample rate in Hz. In addition, the SDFT requires no bitreversal processing as does the FFT. Like the Goertzel algorithm, the SDFT is especially efficient for narrowband spectrum analysis.
For completeness, we mention that a radix2 sliding FFT technique exists for computing all N bins of Xm(q) in Eq. (1385)[47,48]. That technique is computationally attractive because it requires only N complex multiplies to update the Npoint FFT for all N bins; however, it requires 3N memory locations (2N for data and N for twiddle coefficients). Unlike the SDFT, the radix2 sliding FFT scheme requires address bitreversal processing and restricts N to be an integer power of two.
13.18.2 SDFT Stability
The SDFT filter is only marginally stable because its pole resides on the zdomain's unit circle. If filter coefficient numerical rounding error is not severe, the SDFT is boundedinputboundedoutput stable. Filter instability can be a problem, however, if numerical coefficient rounding causes the filter's pole to move outside the unit circle. We can use a damping factor r to force the pole and zeros in Figure 1348(b) to be at a radius of r just slightly inside the unit circle and guarantee stability using a transfer function of
Equation 1396
with the subscript gs meaning guaranteedstable. (Section 7.1.3 provides the mathematical details of moving a filter's poles and zeros inside the unit circle.) The stabilized feedforward and feedback coefficients become –rN and rej2pm/N, respectively. The difference equation for the stable SDFT filter becomes
Equation 1397
with the stabilizedfilter structure shown in Figure 1349. In this case, we perform five real multiplies and four real additions per output sample.
Figure 1349. Guaranteedstable sliding DFT filter structure.
Using a damping factor as in Figure 1349 guarantees stability, but the Xm(q) output, defined by
Equation 1398
is no longer exactly equal to the mth bin of an Npoint DFT in Eq. (1385). While the error is reduced by making r very close to (but less than) unity, a scheme does exist for eliminating that error completely once every N output samples at the expense of additional conditional logic operations[49]. Determining if the damping factor r is necessary for a particular SDFT application requires careful empirical investigation. As is so often the case in the world of DSP, this means you have test your SDFT implementation very thoroughly and carefully!
Another stabilization method worth consideration is decrementing the largest component (either real or imaginary) of the filter's ej2pm/N feedback coefficient by one least significant bit. This technique can be applied selectively to problematic output bins and is effective in combating instability due to rounding errors which result in finiteprecision ej2pm/N coefficients having magnitudes greater than unity. Like the DFT, the SDFT's output is proportional to N, so in fixedpoint binary implementations the designer must allocate sufficiently wide registers to hold the computed results.
13.18.3 SDFT Leakage Reduction
Being equivalent to the DFT, the SDFT also suffers from spectral leakage effects. As with the DFT, SDFT leakage can be reduced by the standard concept of windowing the x(n) input time samples as discussed in Section 3.9. However, windowing by timedomain multiplication would ruin the realtime computational simplicity of the SDFT. Thanks to the convolution theorem properties of discrete systems, we can implement timedomain windowing by means of frequencydomain convolution, as discussed in Section 13.3.
Spectral leakage reduction performed in the frequency domain is accomplished by convolving adjacent Xm(q) values with the DFT of a window function. For example, the DFT of a Hamming window comprises only three nonzero values, –0.23, 0.54, and –0.23. As such, we can compute a Hammingwindowed Xm(q) with a threepoint convolution using
Equation 1399
Figure 1350 shows this process using three resonators, each tuned to adjacent DFT bins (m–1, m, and m+1). The comb filter stage need only be implemented once.
Figure 1350. Threeresonator structure to compute a single Hammingwindowed Xm(q).
Table 135 provides a computational workload comparison of various spectrum analysis schemes in computing an initial Xm(n) value and computing a subsequent Xm(n+1) value.
To compute the initial windowed Xm(n) values in Table 135, the threeterm frequencydomain convolution need only be performed once, upon arrival of the Nth time sample. However, the convolution needs to be performed for all subsequent computations
We remind the reader that Section 13.3 discusses several implementation issues regarding Hanning windowing in the frequency domain, using binary shifts to eliminate the multiplications in Eq. (1399), as well as the use of other window functions.
Table 135. SingleBin DFT Computation Comparison
Method 
Compute initial Xm(n) 
Compute Xm(n+1) 


Real multiplies 
Real adds 
Real multiplies 
Real adds 

DFT 
4N 
2N 
4N 
2N 
Goertzel algorithm 
N + 2 
2N + 1 
N + 2 
2N + 1 
Sliding DFT (marginally stable) 
4N 
4N 
4 
4 
Sliding DFT (guaranteed stable) 
5N 
4N 
5 
4 
Threeterm windowed sliding DFT (marginally stable) 
12N + 6 
10N + 4 
18 
14 
Threeterm windowed sliding DFT (guaranteed stable) 
13N + 6 
10N + 4 
19 
14 
13.18.4 A LittleKnown SDFT Property
The SDFT has a special property that's not widely known, but is very important. If we change the SDFT's comb filter feedforward coefficient (in Figure 1347) from a –1 to +1, the comb's zeros will be rotated counterclockwise around the unit circle by an angle of p/N radians. This situation, for N = 8, is shown on the right side of Figure 1351(a). The zeros are located at angles of 2p(m + 1/2)/N radians. The m = 0 zeros are shown as solid dots. Figure 1351(b) shows the zeros locations for an N = 9 SDFT under the two conditions of the comb filter's feedforward coefficient being –1 and +1.
Figure 1351. Four possible orientations of comb filter zeros on the unit circle.
This alternate situation is useful: we can now expand our set of spectrum analysis center frequencies to more than just N angular frequency points around the unit circle. The analysis frequencies can be either 2pm/N or 2p(m+1/2)/N, where integer m is in the range 0 m N–1. Thus we can build an SDFT analyzer that resonates at any one of 2N frequencies between 0 and fs Hz. Of course, if the comb filter's feedforward coefficient is set to +1, the resonator's feedforward coefficient must be ej2p(m+1/2)/N to achieve pole/zero cancellation.
URL http://proquest.safaribooksonline.com/0131089897/ch13lev1sec18
Amazon  