Sebastian J. Schlecht, Vesa Välimäki, Emanuël A.P. Habets
Analyzing the magnitude response of a finite-length sequence is a ubiquitous task in signal processing. However, the discrete Fourier transform (DFT) provides only discrete sampling points of the response characteristic. This work introduces bounds on the magnitude response, which can be efficiently computed without additional zero padding. The proposed bounds can be used for more informative visualization and inform whether additional frequency resolution or zero padding is required.
The DFT is one of the prime analysis methods for spectral content of a finite-length signal [1]. The DFT, usually computed by applying the fast Fourier transform algorithm, is used for analysis and processing in many applications. Primarily, the magnitude of the DFT is commonly used for signal analysis and visualization purposes, as it describes the strength of the individual frequency components.
Although a DFT of sufficient length carries the complete information of the sequence, direct visualization of the DFT magnitude points can be hard to interpret; see Figure 1 for an example. Many pieces of scientific computing software readily provide plotting features connecting the individual DFT data with line segments. However, such a linear interpolation can be misleading, as the actual shape of the magnitude response is described by the discrete-time Fourier transform (DTFT); see Figure 1. Numerically, the DTFT is usually approximated by a zero-padded higher-order DFT. While sufficient zero padding can lead to higher accuracy, it is often not clear how much zero padding is required.
Figure 1. The magnitude spectrum of a real-valued signal $x(n)$ of length ${N} = {31}$ at nonnegative frequencies. The plot shows the DFT points (blue □) and linear interpolation commonly used for plotting and the corresponding DTFT. The linear interpolation can significantly deviate from the actual magnitude response, occasionally resulting in misleading interpretations.
Some applications require certain bounds on the magnitude spectrum. For instance, a digital filter placed in a feedback loop must have a magnitude response of less than one for stability. Examples for such filters are found in string instrument synthesis [2] and artificial reverberation [3]. Frequency sampling-based filter design is impaired by oscillation between the sampling points, also called the Gibbs phenomenon [4]. The oscillation can be improved by applying appropriate interpolation kernels [5].
Related literature has explored the computation of maximum values of the magnitude response. For example, a theorem by Gronwall [6] gives a bound on the maximum magnitude response of the DTFT if the corresponding DFT points are all bounded by one. It can be shown also, if two trigonometric polynomials of degree N have the same values at ${2}{N} + {2}$ points (or ${2}{N} + {1}$ if zero is repeated), the polynomials are the same, [7, Lemma 3]. Iterative procedures have been devised to compute the maximum magnitude for arbitrary filters [8], [9]. The maximum deviation from piecewise linear approximation is bounded by the maximum value of the second-order derivative [10].
This article proposes an efficient method for computing the upper and lower bounds to the magnitude response. Such bounds can be employed to visualize the continuous magnitude response and make the visual assessment more confident. Further, the bounds can be used for the efficient analysis of extremal response values. A MATLAB implementation of the method and all figures are provided at https://github.com/SebastianJiroSchlecht/bounded_dft.
The DTFT of the sequence $x(n)$ of length N with time ${n} = {0},\ldots,{N}{-}{1}$ is given by \[{X}{(}{f}{)} = \mathop{\sum}\limits_{{n} = {0}}\limits^{{N}{-}{1}}{x}{(}{n}{)}{e}^{\frac{{-}{j}{2}{\pi}}{N}fn} \tag{1} \] where ${f}{/}{N}\in{\mathbb{R}}$ is the continuous frequency in cycles per sample and ${j}^{2} = {-}{1}{.}$ The DFT of $x(n)$ is, then, $X({f}_{k})$ sampled at discrete frequency points ${f}_{k}\in{\mathbb{Z}}{.}$
The magnitude of the DTFT is then given by the pointwise absolute value $\mid{X}{(}{f}{)}\mid{.}$ To avoid the nonlinear absolute value operation, we study the symmetric sequence of length ${2}{N} + {1}$ instead: \[\tilde{x}{(}{n}{)} = {x}{(}{n}{)}\ast{x}{(}{-}{n}{)} \tag{2} \] for ${n} = {-}{N},\ldots,{N},$ where $\ast$ denotes discrete convolution. The corresponding DTFT \[\tilde{X}{(}{f}{)} = {X}{(}{f}{)}{X}{(}{f}{)}^{\ast} = \mid{X}{(}{f}{)}\mid{}^{2} \tag{3} \] is real valued and nonnegative, and it is equal to the power spectrum of $x(n).$ The complex conjugate of $X(f)$ is denoted by ${X}{(}{f}{)}^{\ast}{.}$ For any such symmetric sequence $\tilde{x}(n),$ an x(n) satisfying (2) can be recovered, although the reconstructed phase is not unique; see Edwards [11, Exercise 1.11].
In the following, without loss of generality, we assume that $x(n)$ has a real-valued and nonnegative DTFT and is of odd length N. The DFT is just a sampled DTFT, and, vice versa, the DTFT can be recovered from the DFT by the Dirichlet interpolation [11], i.e., \[\begin{array}{l}{{X}{(}{f}{)} = \mathop{\sum}\limits_{{f}_{k} = {-}{(}{N}{-}{1}{)/}{2}}\limits^{{(}{N}{-}{1}{)/}{2}}{{D}_{N}}{(}{f}{-}{f}_{k}{)}{X}{(}{f}_{k}{)}}\end{array} \tag{4} \] where the (scaled) Dirichlet kernel is \[\begin{array}{l}{{D}_{N}{(}{f}{-}{f}_{k}{)} = \frac{\sin{(}{\pi}{(}{f}{-}{f}_{k}{))}}{{N}\sin\left({\frac{\pi}{N}{(}{f}{-}{f}_{k}{)}}\right)}{.}}\end{array} \tag{5} \]
To establish a bound on X(f), it is then sufficient to establish a bound on the Dirichlet kernel. Ideally, the bound is both tight and fast to compute.
We present lower and upper bounding kernels ${L}_{N}$ and ${U}_{N},$ respectively, such that, for ${f}\in{[}{0},{1}{]}$ and ${f}_{k}\in{[}{-}{N}{/}{2},{N}{/}{2}{],}$ \begin{align*}\begin{array}{l}{\begin{gathered}{\sin{(}{\pi}{f}{)}{L}_{N}{(}{f}_{k}{)}\leq{D}_{N}{(}{f}{-}{f}_{k}{)}} \\{\leq{\sin}{(}{\pi}{f}{\kern0.1em)}{U}_{N}{(}{f}_{k}{).}} \end{gathered}}\end{array} \tag{6} \end{align*}
For the numerator of (5), we have \begin{align*}{\sin}{(}{\pi}{(}{f}{-}{f}_{k}{))} & = \sin{(}{\pi}{f}{)}\cos{(}{\pi}{f}_{k}{)} \\ & \qquad {-}\cos{(}{\pi}{f}{)}\sin{(}{\pi}{f}_{k}{)} \\ & = \sin{(}{\pi}{f}{)}{(}{-}{1}{)}^{{f}_{k}}{.} \tag{7} \end{align*}
Thus, for ${f}\in{[}{0},{1}{],}$ \begin{align*} & \frac{{(}{-}{1}{)}^{{f}_{k}}}{{N}\sin\left(\frac{\pi}{N}{(}{f}{-}{f}_{k}{)}\right)}\leq{U}_{N}{(}{f}_{k}{)} \\ & = \begin{cases}\begin{array}{ll}\frac{1}{{N}\sin\left(\frac{\pi}{N}{(}{0}{-}{f}_{k}{)}\right)}& {\text{for even}}\,{f}_{k} \\ \frac{{-}{1}}{{N}\sin\left(\frac{\pi}{N}{(}{1}{-}{f}_{k}{)}\right)} & {\text{for odd}}\,{f}_{k} \\ \frac{{-}{1}}{{N}\sin\left(\frac{\pi}{N}\left(\frac{1}{2}{-}{f}_{k}\right)\right)} & \begin{array}{l}{\text{for odd}}\\{f}_{k} = {-}{(}{N}{-}{1}{)}{/}{2}\end{array} \\ {\infty} & {\text{for}}\,{f}_{k} = {0},{1}{.} \end{array}\end{cases} \tag{8} \end{align*}
Figure 2(a) shows the bound in (8). The bound is a piecewise constant function tightly placed at the exact curve [the left-hand side of (8)] for each segment.
Figure 2. The Dirichlet interpolation kernels and corresponding upper and lower bound kernels for ${N} = {31}{:}$ the (a) bounded component in (8) and (12) and (b) bounded Dirichlet kernel based on (8) and (12).
The right-hand side of (8) defines, then, the upper bound kernel ${U}_{N}({f}_{k})$ such that \begin{align*} & {D}_{N}{(}{f}{-}{f}_{k}{)}\leq{\sin}{(}{\pi}{f}{)}{U}_{N}{(}{f}_{k}{)} \\ & \qquad {\text{for}}\,{f}_{k}\ne{0},{1}{.} \tag{9} \end{align*}
Figure 2(b) shows the resulting upper bound on the Dirichlet kernel. As expected, the bound is tight near the DFT points but increasingly loose toward the midpoint between DFT points. Also, the closest DFT points have both the largest contribution to the DTFT value and also the largest bounding deviation.
From the upper bound kernel, we derive the DFT upper bound ${X}_{U}(f).$ We split ${f} = {\{}{f}{\}} + \left\lfloor{f}\right\rfloor$ with the fractional part ${\{}{f}{\}} = {f}{-}\left\lfloor{f}\right\rfloor,$ where $\left\lfloor{{\cdot}}\right\rfloor$ denotes the floor function. The DFT is bounded using (8) and computing explicitly the two central points for ${k} = {0},{1},$ i.e., \begin{align*} & {X}{(}{f}{)} = \mathop{\sum}\limits_{{f}_{k} = {-}{(}{N}{-}{1}{)}{/}{2}}^{{(}{N}{-}{1}{)}{/}{2}}{D}_{N}{(}{f}{-}{f}_{k}{)}{X}{(}{f}_{k}{)} \tag{10a} \\ & = \mathop{\sum}\limits_{{f}_{k} = {-}{(}{N}{-}{1}{)}{/}{2}}^{{(}{N}{-}{1}{)}{/}{2}}{D}_{N}\left({\{}{f}{\}} + \left\lfloor{f}\right\rfloor{-}{f}_{k}\right){X}{(}{f}_{k}{)} \tag{10b} \\ & \leq{D}_{N}{(}{\{}{f}{\}}{)}{X}\left(\left\lfloor{f}\right\rfloor\right) + {D}_{N}\left({{\{}{f}{\}}{-}{1}}\right) \\ & {X}\left(\left\lfloor{f}\right\rfloor{-}{1}\right) + \sin{(}{\pi}{\{}{f}{\}}{)}\overline{X}\left(\left\lfloor{f}\right\rfloor\right) = {X}_{U}{(}{f}{)} \tag{10c} \end{align*} with the upper bound DFT \[\begin{array}{l}{\overline{X}\left({\left\lfloor{f}\right\rfloor}\right) = \mathop{\sum}\limits_{{f}_{k} = {-}{(}{N}{-}{1}{)/}{2}}\limits^{{(}{N}{-}{1}{)/}{2}}{{U}_{N}}\left({\left\lfloor{f}\right\rfloor{-}{f}_{k}}\right){X}{(}{f}_{k}{).}}\end{array} \tag{11} \]
Similarly, we can compute a lower bound kernel ${L}_{N}({f}_{k})$ for ${f}\in{[}{0},{1}{]}$ and ${m}\in{\mathbb{N}}{:}$ \begin{align*} & \frac{{(}{-}{1}{)}^{{f}_{k}}}{{N}\sin\left(\frac{\pi}{N}{(}{f}{-}{f}_{k}{)}\right)}\geq{L}_{N}{(}{f}_{k}{)} \\ & = \begin{cases}\begin{array}{ll}\frac{1}{{N}\sin\left(\frac{\pi}{N}{(}{1}{-}{f}_{k}{)}\right)} & {\text{for even}}\,{f}_{k} \\ \frac{{-}{1}}{{N}\sin\left(\frac{\pi}{N}{(}{0}{-}{f}_{k}{)}\right)} & {\text{for odd}}\,{f}_{k} \\ \frac{1}{{N}\sin\left(\frac{\pi}{N}\left(\frac{1}{2}{-}{f}_{k}\right)\right)} & \begin{array}{l}{\text{for even}}\\{f}_{k} = {-}{(}{N}{-}{1}{)}{/}{2}{.}\end{array} \end{array}\end{cases} \tag{12} \end{align*}
Unlike the upper bound ${U}_{N}({f}_{k}),$ the lower bound ${L}_{N}({f}_{k})$ is also finite for ${f}_{k} = {0},{1},$ although only as a poor approximation; see Figure 2(a). Therefore, we use the same neighbor interpolation as for the upper bound; see (10). Figure 2(b) shows the resulting lower bound on the Dirichlet kernel.
Figure 3(a) shows the proposed bounded DFT using the example signal shown in Figure 1. It can be seen that the bounded range (the shaded area) vanishes at the DFT points and becomes larger between them. Note that the lower bound extends into negative values, which is impossible for a magnitude response. The proposed method is generally less suitable for finding good global lower bounds, as the magnitude spectrum can reach zero values. In contrast, the upper bounds are naturally limited because of the finite energy of a finite input signal. The bounds give good estimates on the location of the magnitude response peaks between the sixth and seventh and between the 13th and 14th frequency bins. In particular, the bounds indicate that the filter is not dampening; i.e., there exists ${f}$ with $\kern-0.0em\left|{X(\kern0.3emf\kern0.4em)}\right|\kern0.2em{>}{\kern0.2em1}{\kern0.2em,}$ while this assessment is not readily available from the DFT points $X({f}_{k}).$ Figure 3(b) shows the proposed bounds on a classic brickwall filter, where $x(n)$ is a truncated sinc function. The upper bound indicates the maximum peak magnitude value of around 1.1, while the true maximum peak value is 0.99. Also, in this case, the bounds are close to the DTFT curve, particularly around the cutoff frequency at the eighth bin.
Figure 3. The magnitude response of two example signals $x(n)$ of length ${N} = {31}{.}$ The plot shows the DFT points (blue □) and bounded DFT by the proposed method (the shaded area). The upper and lower bounds clearly indicate the bounded range of the DFT and support sound interpretation of the spectrum. (a) The random sequence $x(n),$ the same as in Figure 1. (b) The truncated Sinc function ${x}{(}{n}{)}{.}$
The proposed DFT bounds can be computed efficiently. Instead of the frequency-domain convolution in (11), the upper bound $\overline{X}\left({\left\lfloor{f}\right\rfloor}\right)$ is computed with a pointwise multiplication in the time domain and a single DFT of size N; i.e., \[{\overline{X}} = {\text{DFT}}_{N}{(}{x}\,{\circ}\,{u}_{N}{)} \tag{13} \] where ${\circ}$ denotes the Hadamard product, and ${u}_{N} = {iDFT}_{N}{(}{U}_{N}{)}$ is the inverse DFT, which can be precomputed if reusing the same transformation lengths repeatedly. The evaluation of specific points of the DFT bounds takes only local computation of the three terms in (10c), which is independent of the overall signal lengths. The overall computation graph for the upper bound is shown in Figure 4.
Figure 4. The proposed method of computing the upper bound of the DFT based on (10). The precomputation requires two DFTs, whereas the evaluation at frequency $f$ takes only computing the three terms on the right-hand side.
We proposed an efficient method to compute a discrete magnitude spectrum’s upper and lower bounds using bounding kernels and two DFTs. The mathematical derivation is given for symmetric signals with real-valued spectra; however, the method is generally applicable to any finite signals. The bounded DFT indicates the possible range in data visualization with a given size analysis order. The bounded DFT can also display whether the given resolution sufficiently resolves the features in the magnitude response or whether additional zero padding is necessary. The bound can naturally be further improved by computing more closest neighbors interpolations in (10), although at the price of more computation.
Sebastian J. Schlecht (sebastian.schlecht@aalto.fi) received his diploma in applied mathematics from the University of Trier, Germany, in 2010; his M.Sc. degree in digital music processing from Queen Mary University of London, U.K., in 2011; and his doctoral degree on the topic of artificial spatial reverberation and reverberation enhancement systems from the International Audio Laboratories Erlangen, Germany, in 2017. He is a professor of practice for sound in virtual reality at the Acoustics Lab, Department of Information and Communications Engineering and Media Labs, Department of Art and Media, Aalto University, 00076 Espoo, Finland. He is a Senior Member of IEEE.
Vesa Välimäki (vesa.valimaki@aalto.fi) is a full professor at the Acoustics Lab, Department of Information and Communications Engineering, Aalto University, 0076 Espoo, Finland. In 2008–2009, he was a visiting scholar at Stanford University. He was the chair of the International Conference on Digital Audio Effects in 2008 and the chair of the Sound and Music Computing Conference in 2017. He is the editor-in-chief of Journal of the Audio Engineering Society. His research interests include signal processing and machine learning applied to audio and music technology. He is a Fellow of IEEE and a fellow of the Audio Engineering Society.
Emanuël A.P. Habets (emanuel.habets@audiolabs-erlangen.de) is a full professor at the International Audio Laboratories Erlangen [a joint institution of the Friedrich-Alexander-Universität Erlangen-Nürnberg and Fraunhofer Institute for Integrated Circuits (IIS)], and head of the Department of Conversational Artificial Intelligence Research at Fraunhofer IIS, 91058 Erlangen, Germany. He served as cochair of the 2013 International Workshop on Applications of Signal Processing to Audio and Acoustics. He is a founding member of the European Association For Signal Processing (EURASIP) Acoustic, Speech, and Music Signal Processing Technical Area Committee, where he served as a vice-chair (2015–2018) and chair (2019–2021). From 2016 until 2018, he was the editor-in-chief of EURASIP Journal on Audio, Speech, and Music Processing. He is a Senior Member of IEEE.
[1] J. R. Deller, “Tom, Dick, and Mary discover the DFT,” IEEE Signal Process. Mag., vol. 11, no. 2, pp. 36–50, Apr. 1994, doi: 10.1109/79.273893.
[2] B. Bank and V. Välimäki, “Robust loss filter design for digital waveguide synthesis of string tones,” IEEE Signal Process. Lett., vol. 10, no. 1, pp. 18–20, Jan. 2003, doi: 10.1109/LSP.2002.806707.
[3] S. J. Schlecht and E. A. P. Habets, “Accurate reverberation time control in feedback delay networks,” in Proc. Int. Conf. Digit. Audio Effects (DAFx), Edinburgh, U.K., Aug. 2017, pp. 337–344.
[4] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed. Englewood Cliffs, NJ, USA: Prentice-Hall, 2010.
[5] X. Huang, S. Jing, Z. Wang, Y. Xu, and Y. Zheng, “Closed-form FIR filter design based on convolution window spectrum interpolation,” IEEE Trans. Signal Process., vol. 64, no. 5, pp. 1173–1186, Mar. 2016, doi: 10.1109/TSP.2015.2494869.
[6] T. H. Gronwall, “A sequence of polynomials connected with the nth roots of unity,” Bull. Amer. Math. Soc., vol. 27, no. 6, pp. 275–279, Mar. 1921, doi: 10.1090/S0002-9904-1921-03411-2.
[7] R. P. Boas, “Inequalities for the derivatives of polynomials,” Math. Mag., vol. 42, no. 4, pp. 165–174, Sep. 1969, doi: 10.2307/2688534.
[8] J. J. Green, “Calculating the maximum modulus of a polynomial using Steckin’s lemma,” SIAM J. Numer. Anal., vol. 36, no. 4, pp. 1022–1029, Jan. 1999, doi: 10.1137/S0036142997331335.
[9] G. De La Chevrotière, “Finding the maximum modulus of a polynomial on the polydisk using a generalization of Steckin’s lemma,” SIAM Undergraduate Res. Online, vol. 2, no. 2, pp. 71–80, Jan. 2009, doi: 10.1137/09S010460.
[10] C. L. Frenzen, T. Sasao, and J. T. Butler, “On the number of segments needed in a piecewise linear approximation,” J. Comput. Appl. Math., vol. 234, no. 2, pp. 437–446, May 2010, doi: 10.1016/j.cam.2009.12.035.
[11] R. E. Edwards, Fourier Series, vol. 64. New York, NY, USA: Springer-Verlag, 1979.
Digital Object Identifier 10.1109/MSP.2022.3228526