2. Single-bin sliding discrete Fourier transform

The discrete Fourier transform (DFT) is a numerical approximation of the theoretical Fourier transform (FT) of a continuous and infinite duration signal. It represents the most common tool for engineers to extract the frequency content of a finite and discrete signal sequence, obtained from the periodic sampling of a continuous wave form in time domain.

Let us consider a continuous time signal xðtÞ that is sampled at the rate f <sup>s</sup> ¼ N · f <sup>o</sup> (where f <sup>o</sup> is the fundamental frequency of xðtÞ) to produce the time sequence x½n�. Then the DFT of the sequence x½n� is defined as:

$$X(k) = \sum\_{n=0}^{N-1} x[n] \mathcal{W}\_N^{-k \text{ n}} \tag{1}$$

where <sup>X</sup>ðk<sup>Þ</sup> is the DFT output coefficient, WN <sup>¼</sup> <sup>e</sup><sup>j</sup>2π=<sup>N</sup> is the complex twiddle factor, <sup>N</sup> is the sequence length, k is the frequency domain index ð0 ≤ k ≤ N−1Þ, and n is the time domain index [1].

If Eq. (1) is not properly designed and implemented, the DFT calculation in real-time might represent a considerable bottleneck when developing a DFT-based estimation algorithm, in terms of both measurement reporting latencies and achievable reporting rates. In this respect, in order to improve both latencies and throughput, several efficient techniques to compute the DFT spectrum have been proposed in literature, which can be classified as nonrecursive and recursive algorithms. Among the nonrecursive class, the fast Fourier transform (FFT) algorithm is extensively used for harmonic analysis over an extended portion of the spectrum. When, on the other hand, only a subset of the overall DFT spectrum is necessary to accomplish the desired estimate, the so-called single-bin sliding DFT (Sb-SDFT) turns out to be very effective.

The DFT can also be computed by recursive algorithms which are characterized by a minor number of operations to calculate a single DFT bin. Regardless of this advantage with respect to the class of nonrecursive algorithms, the performances of the two categories usually are not the same. Especially, most of the algorithms in the recursive category suffers of errors due to either the approximations made to perform the recursive update or the accumulation of the quantization errors related to a finite word-length precision [2, 3].

In what follows, four of the most efficient techniques to compute a portion of the DFT spectrum, namely the sliding discrete Fourier transform (SDFT), the sliding Goertzel transform (SGT), the Douglas and Soh algorithm (D&S), and the modulated sliding discrete Fourier transform (mSDFT) will be presented and described.

### 2.1. Sliding discrete Fourier transform

solved using an adaptive coherent sampling mechanism. One of these mechanisms is known as variable sampling period technique (VSPT) and is characterized for the dynamic adjustment of the sampling period to exactly N times the fundamental frequency, thereby avoiding the

The chapter is organized as follows: Section 2 presents a brief review of Sb-SDFT. Section 3 evaluates and compares the four selected Sb-SDFT algorithms in diverse operational conditions, identifying the similarities between them. In order to mitigate the inaccuracies resulting from the spectral leakage effect, a scheme for coherent sampling based on VSPT is introduced in Section 4. Altogether a unified model is also presented to generalize this scheme to all Sb-SDFT along with simulation results. Finally, the conclusions of this chapter are drawn in

The discrete Fourier transform (DFT) is a numerical approximation of the theoretical Fourier transform (FT) of a continuous and infinite duration signal. It represents the most common tool for engineers to extract the frequency content of a finite and discrete signal sequence, obtained

Let us consider a continuous time signal xðtÞ that is sampled at the rate f <sup>s</sup> ¼ N · f <sup>o</sup> (where f <sup>o</sup> is the fundamental frequency of xðtÞ) to produce the time sequence x½n�. Then the DFT of the

where <sup>X</sup>ðk<sup>Þ</sup> is the DFT output coefficient, WN <sup>¼</sup> <sup>e</sup><sup>j</sup>2π=<sup>N</sup> is the complex twiddle factor, <sup>N</sup> is the sequence length, k is the frequency domain index ð0 ≤ k ≤ N−1Þ, and n is the time domain

If Eq. (1) is not properly designed and implemented, the DFT calculation in real-time might represent a considerable bottleneck when developing a DFT-based estimation algorithm, in terms of both measurement reporting latencies and achievable reporting rates. In this respect, in order to improve both latencies and throughput, several efficient techniques to compute the DFT spectrum have been proposed in literature, which can be classified as nonrecursive and recursive algorithms. Among the nonrecursive class, the fast Fourier transform (FFT) algorithm is extensively used for harmonic analysis over an extended portion of the spectrum. When, on the other hand, only a subset of the overall DFT spectrum is necessary to accomplish the desired estimate, the so-called single-bin sliding DFT (Sb-SDFT) turns out to be very

The DFT can also be computed by recursive algorithms which are characterized by a minor number of operations to calculate a single DFT bin. Regardless of this advantage with respect to the class of nonrecursive algorithms, the performances of the two categories usually are not

<sup>x</sup>½n�W<sup>−</sup><sup>k</sup> <sup>n</sup> <sup>N</sup> (1)

XðkÞ ¼ ∑ N−1 n¼0

above-mentioned problems.

26 Fourier Transforms - High-tech Application and Current Trends

sequence x½n� is defined as:

2. Single-bin sliding discrete Fourier transform

from the periodic sampling of a continuous wave form in time domain.

Section 5.

index [1].

effective.

A very effective Sb-SDFT method for sample-by-sample DFT bin computation is the so-called sliding discrete Fourier transform (SDFT) technique [4]. Starting from Eq. (1), the DFT can be potentially updated every time-step n, based on the most recent set of samples within a sliding window {x½n�N þ 1�, x½n�N þ 2�, …;x½n�}. The time window is advanced one sample at a time, and a new N-point DFT is calculated. Figure 1(a) illustrates the time domain indexing within the sliding window by showing the input samples used to compute k-bin of an N-points DFT when n ¼ no. The principle used for SDFT is known as the DFT shifting theorem, or the circular shift property [1].

Based on this property, the SDFT can be recursively implemented to calculate Eq. (1) for a desired k-bin, as:

$$X\_k\left[n\right] = \mathcal{W}\_N^k X\_k\left[n-1\right] \text{-x}[n-N] + \text{x}[n] \tag{2}$$

where Xk½n� is calculated by phase shifting the sum of the previous Xk½n�1� with the difference between the current and delayed input sample, x½n� and x½n�N�, respectively [4, 5]. The complex output of the SDFT could be rewritten as:

$$\mathbf{X}\_{k}\begin{bmatrix}n\end{bmatrix} = \mathbf{X}\_{\mathbf{r}k}\begin{bmatrix}n\end{bmatrix} + \mathbf{j}\mathbf{X}\_{\mathbf{ik}}\begin{bmatrix}n\end{bmatrix} \tag{3}$$

where Xr<sup>k</sup> ½n� and Xi<sup>k</sup> ½n� are real and imaginary components of the DFT output coefficient, respectively. The SDFT provides an accurate estimation for the kth component as its amplitude (Ak ½n�) and phase (ϕ<sup>k</sup> ½n�) can be determined by computing the modulus and the argument of the complex result Xk ½n�, as stated by

Figure 1. (a) Samples used to compute Xk½n� within a sliding window, when n ¼ no. (b) Guaranteed-stable SDFT implementation as IIR filter as given by (5).

$$A\_k[n] = \frac{2}{N} \text{ abs}\left(X\_k\left[n\right]\right) \tag{4a}$$

$$\varphi\_k[n] = \arg\left(X\_k\left[n\right]\right) \tag{4b}$$

SDFT is computationally efficient, as it only requires one (complex) multiplication and two additions per time instant. Nevertheless, the implementation of Eq. (2) as an infinite impulse response (IIR) filter in a system with finite word-length precision brings about a rounding error in the implementation of the W<sup>k</sup> <sup>N</sup> coefficient, which may turn the algorithm unstable and/ or increment the estimation error. The first one is a direct consequence of wrong cancellations between singularities and by poles displacement outside the unit circle [2, 3]. Commonly, a damping factor (r, with 0 < r < 1) is used to ensure that all singularities are placed inside the unit circle, hence instability is no longer an issue. Then, the intrinsically stable version of the SDFT is

$$\tilde{X}\_k\left[n\right] = r\mathcal{W}\_N^k \tilde{X}\_k\left[n-1\right] \neg r^N \mathbf{x}[n-N] + \mathbf{x}[n] \tag{5}$$

where <sup>X</sup><sup>~</sup> <sup>k</sup> <sup>½</sup>n� is the estimated DFT output coefficient. While Eq. (5) is numerically stable, it no longer computes the exact value of XðkÞ in Eq. (1), since a small error is induced by the damping factor. The z domain transfer function for the estimated kth bin of the SDFT is

$$H\_{\rm S\,DFT}(z) = \frac{1 - r^N z^{-N}}{1 - rW\_N^k z^{-1}}\tag{6}$$

The stable SDFT algorithm given by Eq. (5) leads to the filter structure shown in Figure 1(b). This structure is basically an IIR filter that comprises a comb filter followed by a complex resonator. The comb filter makes the transient response N−1 samples in length; therefore, the output will reach steady state when the stored waveform equals the input signal.

### 2.2. Sliding Goertzel transform

The number of multiplications required in the SDFT can be reduced by creating a new pole/ zero pair in its HSDFTðzÞ system function. This is achieved by multiplying the numerator and denominator of <sup>H</sup>SDFTðz<sup>Þ</sup> in Eq. (6) by the factor <sup>ð</sup>1−rW<sup>−</sup><sup>k</sup> <sup>N</sup> <sup>z</sup><sup>−</sup><sup>1</sup><sup>Þ</sup> yielding:

$$H\_{\rm SGT}(z) = \frac{(1 - rW\_N^{-k}z^{-1})(1 - r^Nz^{-N})}{1 - 2r\cos\left(2\pi k/N\right)z^{-1} + r^2z^{-2}}\tag{7}$$

The transfer function represented by Eq. (7) is commonly known as the sliding Goertzel transform (SGT). Because the poles are placed on the z domain unit circle, the SGT implementation is also potentially unstable. Once more a damping factor r can be used in Eq. (7), to move the singularities inside the unit circle and to ensure the system stability.

This method can be implemented by the following pair of finite difference equations:

$$\boldsymbol{\sigma} \cdot \boldsymbol{\mathfrak{v}} \, [\boldsymbol{n}] = \mathbb{C}\_1 \, \boldsymbol{\mathfrak{v}} \, [\boldsymbol{n} - 1] - \mathbb{C}\_2 \, \boldsymbol{\mathfrak{v}} \, [\boldsymbol{n} - 2] + \boldsymbol{\mathfrak{x}}[\boldsymbol{n}] - \boldsymbol{r}^N \boldsymbol{\mathfrak{x}}[\boldsymbol{n} - \boldsymbol{N}] \, \tag{8a}$$

$$
\dot{X}\_k \begin{bmatrix} n \end{bmatrix} = \boldsymbol{\upsilon} \begin{bmatrix} n \end{bmatrix} - r \boldsymbol{\mathcal{W}}\_N^{-k} \boldsymbol{\upsilon} \begin{bmatrix} n-1 \end{bmatrix} \tag{8b}
$$

where <sup>C</sup><sup>1</sup> <sup>¼</sup> <sup>2</sup><sup>r</sup> cos <sup>ð</sup>2πk=N<sup>Þ</sup> and <sup>C</sup><sup>2</sup> <sup>¼</sup> <sup>r</sup>2, with 0 <sup>&</sup>lt; <sup>r</sup> <sup>&</sup>lt; 1. The SGT is implemented as an IIR filter that consists of a comb filter followed by the standard Goertzel filter, as depicted in Figure 2(a). The resulting system only has real coefficients so its computational complexity is decreased in relation to that of the SDFT [6, 7].

### 2.3. Douglas and Soh algorithm

Ak½n� ¼ <sup>2</sup>

<sup>X</sup><sup>~</sup> <sup>k</sup> <sup>½</sup>n� ¼ rW<sup>k</sup>

error in the implementation of the W<sup>k</sup>

28 Fourier Transforms - High-tech Application and Current Trends

2.2. Sliding Goertzel transform

denominator of <sup>H</sup>SDFTðz<sup>Þ</sup> in Eq. (6) by the factor <sup>ð</sup>1−rW<sup>−</sup><sup>k</sup>

SDFT is

SDFT is computationally efficient, as it only requires one (complex) multiplication and two additions per time instant. Nevertheless, the implementation of Eq. (2) as an infinite impulse response (IIR) filter in a system with finite word-length precision brings about a rounding

or increment the estimation error. The first one is a direct consequence of wrong cancellations between singularities and by poles displacement outside the unit circle [2, 3]. Commonly, a damping factor (r, with 0 < r < 1) is used to ensure that all singularities are placed inside the unit circle, hence instability is no longer an issue. Then, the intrinsically stable version of the

<sup>N</sup>X<sup>~</sup> <sup>k</sup> <sup>½</sup>n�1�−<sup>r</sup>

where <sup>X</sup><sup>~</sup> <sup>k</sup> <sup>½</sup>n� is the estimated DFT output coefficient. While Eq. (5) is numerically stable, it no longer computes the exact value of XðkÞ in Eq. (1), since a small error is induced by the damping factor. The z domain transfer function for the estimated kth bin of the SDFT is

<sup>H</sup>S DFTðzÞ ¼ <sup>1</sup>−rNz<sup>−</sup><sup>N</sup>

The stable SDFT algorithm given by Eq. (5) leads to the filter structure shown in Figure 1(b). This structure is basically an IIR filter that comprises a comb filter followed by a complex resonator. The comb filter makes the transient response N−1 samples in length; therefore, the

The number of multiplications required in the SDFT can be reduced by creating a new pole/ zero pair in its HSDFTðzÞ system function. This is achieved by multiplying the numerator and

The transfer function represented by Eq. (7) is commonly known as the sliding Goertzel transform (SGT). Because the poles are placed on the z domain unit circle, the SGT implementation is also potentially unstable. Once more a damping factor r can be used in Eq. (7), to move

This method can be implemented by the following pair of finite difference equations:

<sup>N</sup> <sup>z</sup><sup>−</sup><sup>1</sup><sup>Þ</sup> yielding:

<sup>1</sup>−2<sup>r</sup> cos <sup>ð</sup>2πk=NÞz<sup>−</sup><sup>1</sup> <sup>þ</sup> <sup>r</sup><sup>2</sup>z<sup>−</sup><sup>2</sup> (7)

<sup>N</sup> <sup>z</sup><sup>−</sup><sup>1</sup>Þð1−rNz<sup>−</sup><sup>N</sup><sup>Þ</sup>

output will reach steady state when the stored waveform equals the input signal.

<sup>H</sup>SGTðzÞ ¼ <sup>ð</sup>1−rW<sup>−</sup><sup>k</sup>

the singularities inside the unit circle and to ensure the system stability.

1−rW<sup>k</sup>

<sup>N</sup> abs <sup>ð</sup>Xk <sup>½</sup>n�Þ (4a)

ϕk½n� ¼ arg ðXk ½n�Þ (4b)

<sup>N</sup> coefficient, which may turn the algorithm unstable and/

Nx½n�N� þ <sup>x</sup>½n� (5)

Nz<sup>−</sup><sup>1</sup> (6)

The implementation of a SDFT or SGT requires a damping factor to guarantee the algorithm stability. The trade-off for the system stability is that the calculated value is no longer exactly equal to the kth-bin of an N-point DFT in Eq. (1). In Ref. [8], a technique that significantly reduces this error, without compromising the stability, is developed. This method is a periodically time-varying system designed to generate an <sup>X</sup><sup>~</sup> <sup>k</sup> <sup>½</sup>n� output that is mathematically equal to XðkÞ in Eq. (1) at every Nth time instant.

This technique is implemented by the following pair of finite difference equations:

$$\boldsymbol{\tilde{X}}\_{k}[\boldsymbol{n}] = \begin{cases} r\mathcal{W}\_{N}^{k} \boldsymbol{\tilde{X}}\_{k}[\boldsymbol{n}-1] \text{-} r\mathbf{x}[\boldsymbol{n}-N] + \mathbf{x}[\boldsymbol{n}], & (\boldsymbol{n} \bmod \text{N}) = \mathbf{0} \quad (\text{a}) \\\mathcal{W}\_{N}^{k} \boldsymbol{\tilde{X}}\_{k}[\boldsymbol{n}-1] \text{-} r\mathbf{x}[\boldsymbol{n}-N] + \mathbf{x}[\boldsymbol{n}], & \text{else} \end{cases} \tag{9}$$

The algorithm described by Eq. (9) will be referred to as the Douglas and Soh algorithm (D&S). The filter implementation of Eq. (9), shown in Figure 2(b), requires two multiplications and two additions as well as the control logics to determine when n mod N ¼ 0. In the figure, the

Figure 2. (a) Guaranteed-stable SGT implementation as IIR filter as given by (8). (b) Guaranteed-stable D&S algorithm implemented as IIR filter as given by (9).

change between Eqs. (9a) and (9b) is performed by switch S<sup>1</sup> Therefore, the switching period of S<sup>1</sup> in Figure 2(b) is equal to N · Ts, where T<sup>s</sup> is the sampling period, and its duty cycle is equal to one sample. It is worth mentioning that the effect of the nonlinear operation of D&S algorithm in the dynamic response is negligible as it only changes its structure every N samples.

### 2.4. Modulated sliding discrete Fourier transform

There is an alternative way of avoiding the reduction in accuracy generated by the damping factor, without compromising stability. SDFT implementation in Eq. (2) is marginally stable, however, for the particular case of k ¼ 0 (DC component estimation). It takes the following form:

$$X\_0[n] = X\_0[n-1] \text{-} \mathbf{x}[n-N] + \mathbf{x}[n] \tag{10}$$

The absence of the W<sup>k</sup> <sup>N</sup> coefficient, which typically leads to stability issues when it is represented with finite precision, allows to implement the recursive expression without the damping factor r. Therefore, the recurrence in Eq. (10) is unconditionally stable and does not accumulate errors. The modulated sliding discrete Fourier transform (mSDFT) algorithm uses the Fourier modulation property to effectively shift the DFT bin of interest to the position k ¼ 0 and then use Eq. (10) for computing that DFT bin output. This is accomplished by the multiplication of the input signal <sup>x</sup>½n� by the modulation sequence <sup>W</sup><sup>−</sup><sup>k</sup> <sup>n</sup> <sup>N</sup> . This approach allows to exclude the complex twiddle factor from the resonator and avoids accumulated errors and potential instabilities [9]. The recursive realization of the mSDFT is:

$$X\_k^0 \left[ n \right] = X\_k^0 \left[ n - 1 \right] - \mathbf{x} \left[ n - \mathbf{N} \right] \mathcal{W}\_N^{-k(n-N)} + \mathbf{x} \left[ n \right] \mathcal{W}\_N^{-k \; n} \tag{11a}$$

$$X\_k\left[n\right] = \mathcal{W}\_N^{k\,n} X\_k^0\left[n\right] \tag{11b}$$

where X<sup>0</sup> <sup>k</sup> ½n� is a complex constant related to the phase of the complex twiddle factor, since the modulation moves the desired kth-bin to k ¼ 0 (0 Hz). The relation between the desired Xk ½n� and the computed X<sup>0</sup> <sup>k</sup> ½n� is given by Eq. (11b). It is worth noticing that if the application only requires DFT magnitude estimation, the complex multiplication in Eq. (11b) is unnecessary because <sup>j</sup>X<sup>0</sup> <sup>k</sup> j is equal to jXðkÞj. The filter structure of the mSDFT algorithm in Eq. (11) is depicted in Figure 3(a). In contrast of traditional recursive DFT algorithms, the mSDFT method

Figure 3. (a) Guaranteed-stable mSDFT implementation as IIR filter as given by (11). (b) Guaranteed-stable mSDFT implementation as IIR filter as given by (12).

is unconditionally stable and does not accumulate errors because its singularities are exactly placed on the unit circle, regardless of the finite precision used. These advantages are possible due to the removal of the complex twiddle factor from the resonator loop.

If multiple DFT frequency bins are to be computed, the mSDFT in Eq. (11) requires a comb filter for each frequency bin. On the other hand, given the periodicity of <sup>W</sup><sup>−</sup><sup>k</sup> <sup>n</sup> <sup>N</sup> , as shown in Ref. [9], Eq. (11) can be rewritten as

$$X\_k^0 \left[ n \right] = X\_k^0 \left[ n - 1 \right] + W\_N^{-k} \left( -\mathbf{x}[n - N] + \mathbf{x}[n] \right) \tag{12a}$$

$$X\_k\left[\mathfrak{n}\right] = \mathcal{W}\_N^{k\,\,\,n} X\_k^0\left[\mathfrak{n}\right] \tag{12b}$$

Whenever multiple DFT frequency bins are to be computed, Eq. (12) becomes a more efficient approach as only one comb filter is needed (Figure 3(b)).
