Iterative execution of discrete and inverse discrete Fourier transforms with applications for signal denoising via sparsification

H. Robert Frost¹

Abstract

We describe a family of iterative algorithms that involve the repeated execution of discrete and inverse discrete Fourier transforms. One interesting member of this family is motivated by the discrete Fourier transform uncertainty principle and involves the application of a sparsification operation to both the real domain and frequency domain data with convergence obtained when real domain sparsity hits a stable pattern. This sparsification variant has practical utility for signal denoising, in particular the recovery of a periodic spike signal in the presence of Gaussian noise. General convergence properties and denoising performance relative to existing methods are demonstrated using simulation studies. An R package implementing this technique and related resources can be found at https://hrfrost.host.dartmouth.edu/IterativeFT.

¹Department of Biomedical Data Science
Dartmouth College
Hanover, NH 03755, USA
[email protected]

1 Problem statement

We consider a class of iterative discrete Fourier transform [1] techniques described by Algorithm (1). The structure of this algorithm is broadly motivated by iterative methods such as the alternating direction method of multipliers (ADMM) [2] and expectation-maximization (EM) [3]. The family of algorithms we consider take a real-valued vector $\mathbf{x}\in\mathbb{R}^{n}$ as input and then repeatedly perform the following sequence of actions:

•

Execute a function $h():\mathbb{R}^{n}\to\mathbb{R}^{n}$ on $\mathbf{x}$ .
•

Perform a discrete Fourier transform, $dft()$ , on the output of $h()$ . The $j^{th}$ element of the complex-valued vector output by the discrete Fourier transform, $dft(\mathbf{x})[j]$ , is defined as:

$dft(\mathbf{x})[j]=\sum_{k=1}^{n}\mathbf{x}[k]e^{-i2\pi jk/n}$ (1)
•

Execute a function $g():\mathbb{C}^{n}\to\mathbb{C}^{n}$ on the complex vector output by the discrete Fourier transform.
•

Transform the output of $g()$ back into the real domain via the inverse discrete Fourier transform, $dft^{-1}()$ . The $j^{th}$ element of the real-valued vector output by the inverse discrete Fourier transform, $dft^{-1}(\mathbf{c})[j]$ , is defined as:

$dft^{-1}(\mathbf{c})[j]=\frac{1}{n}\sum_{k=1}^{n}\mathbf{c}[k]e^{i2\pi jk/n}$ (2)

This iteration can expressed as:

\mathbf{x}_{k}=\begin{cases}h(\mathbf{x}_{k})\ &k=0\\ h(dft^{-1}(g(dft(\mathbf{x}_{k-1}))))&k\geq 1\\ \end{cases}

(3)

Convergence of the algorithm is determined by a function $c()$ that compares $\mathbf{x}_{k}$ to $\mathbf{x}_{k-1}$ , i.e., the output of $h()$ on the current iteration to the version from the prior iteration. When convergence is obtained, the output from the last execution of $h()$ is returned. See Algorithm (1) for a detailed definition. For this general family of algorithms, a key question relates to what combinations of $h()$ , $g()$ , and $c()$ functions and constraints on $\mathbf{x}$ enable convergence. We are specifically interested in scenarios involving convergence to a value that is relevant for a specific data analysis application, e.g., a denoising or optimization problem.

Algorithm 1 Iterative application of discrete Fourier and inverse Fourier transforms

Inputs:

•

$\mathbf{x}\in\mathbb{R}^{n}$ $\triangleright$ Input data
•

$i_{m}$ $\triangleright$ The maximum number of iterations

Outputs:

•

$\mathbf{y}\in\mathbb{R}^{n}$ $\triangleright$ Output data
•

$i_{c}$ $\triangleright$ Number of iterations completed

i=1

\triangleright

Initialize iteration index

\mathbf{x}_{0}=\mathbf{x}

\triangleright

Initialize

\mathbf{x}_{i}

3:while

i\leq i_{m}

\mathbf{x}^{*}_{i}=h(\mathbf{x}_{i-1})

\triangleright

Apply function

h():\mathbb{R}^{n}\to\mathbb{R}^{n}

\mathbf{x}_{i-1}

5: if

c(\mathbf{x}^{*}_{i},\mathbf{x}^{*}_{i-1})

then

\triangleright

Check for convergence using indicator function

c()

with domain

\{\mathbb{R}^{n},\mathbb{R}^{n}\}

6: break

\triangleright

If convergence conditions met, stop the iteration

\mathbf{w}_{i}=\textit{dft}(\mathbf{x}^{*}_{i})

\triangleright

Compute discrete Fourier transform of

\mathbf{x}^{*}_{i}

\mathbf{w}^{*}_{i}=g(\mathbf{w}_{i})

\triangleright

Apply function

g():\mathbb{C}^{n}\to\mathbb{C}^{n}

to complex vector

\mathbf{w}_{i}

\mathbf{x}_{i}=\textit{dft}^{-1}(\mathbf{w}^{*}_{i})

\triangleright

Compute inverse discrete Fourier transform of

\mathbf{w}^{*}_{i}

10:

i=i+1

\triangleright

Increment iteration index return

(\mathbf{x}^{*}_{i},i)

\triangleright

Return the output of the final execution of

h()

1.1 Trivial cases

If both $h()$ and $g()$ are the identity function, then convergence occurs after a single iteration and the entire algorithm operates as the identity function, i.e., $h(dft^{-1}(g(dft(\mathbf{x})))=dft^{-1}(dft(\mathbf{x}))=\mathbf{x}$ . Similarly, if just one of $h()$ or $g()$ is the identity function, then the algorithm simplifies to the repeated execution of the non-identity function, e.g., if $g()$ is the identity function then $h(dft^{-1}(g(dft(\mathbf{x})))=h(dft^{-1}(dft(\mathbf{x})))=h(\mathbf{x})$ . In general, we will assume that neither $h()$ nor $g()$ are the identity function and that the number of iterations until convergence, $i_{c}$ , is a function of $\mathbf{x}$ , i.e., if $\mathbf{x}$ is a random vector, then $i_{c}$ is a random variable.

1.2 Generalizations

A number of generalizations of Algorithm (1) are possible:

1.

Matrix-valued input: Instead of accepting a vector $\mathbf{x}\in\mathbb{R}^{n}$ as input, the algorithm could accept a matrix $\mathbf{X}\in\mathbb{R}^{n\times m}$ (or higher-dimensional array) with $\textit{dft}()$ and $\textit{dft}^{-1}()$ replaced by two-dimensional counterparts.
2.

Complex-valued input: Instead of just real values, elements of the input could be allowed to take complex or hypercomplex (e.g., quaternion or octonion) values with a correponding change to the complex, or hypercomplex, discrete Fourier transform.
3.

Alternative invertable discrete transform: The discrete Fourier transform could be replaced by another invertable discrete transform, e.g., discrete wavelet transform [4]. More broadly, a similar approach could be used with any invertable discrete function.

We explore the first generalization in this paper (see Section 4.2) but restrict our interest to the real-valued discrete Fourier transform and leave the other generalizations to future work.

2 Related techniques

Algorithm (1) is related to both standard discrete Fourier analysis techniques and iterative algorithms that alternate between dual problem representations. Although a detailed comparison is not possible without defining $h()$ , $g()$ , and $c()$ , some general observations are possible.

2.1 Relationship to standard discrete Fourier analysis

The discrete Fourier transform is widely used in many areas of data analysis with common applications including signal processing [5], image analysis [6], and efficient evaluation of convolutions [7]. The standard structure for these applications is the transformation of real-valued data (e.g, time-valued signal or image) into the frequency domain via a discrete Fourier transform, computation on the frequency domain representation, and final transformation back into the time/location domain via an inverse discrete Fourier transform. For example, to preserve specific frequencies in a signal, the original data can be transformed into the frequency domain via a discrete Fourier transform, the coefficients for non-desired frequencies set to zero, and a filtered time domain signal generated via an inverse transformation. A similar structure is found in the application of other invertable discrete transforms. A key difference between this type of application and the approach discussed in this paper relates to number of executions of the forward and inverse transforms and the nature of the functions applied to the data in each domain, i.e., $h()$ and $g()$ . Specifically, standard discrete Fourier analysis involves just the single application of forward and inverse transforms separated by execution of a single function. In contrast, the method described by Algorithm (1) involves the repeated execution of forward and inverse transforms interleaved by the execution of $h()$ and $g()$ functions until the desired convergence conditions are obtained.

2.2 Relationship to other iterative algorithms

Algorithm (1) shares features with a range of methods (e.g., ADMM [2], Dykstra’s algorithm [8], and EM [3]) that alternate between coupled representations of a problem on each iteration. For ADMM, an optimization problem is solved by alternating between 1) estimating the value of the primal variable that minimizes the Lagrangian with the Lagrangian variables fixed and 2) updating the Lagrangian or dual variables using the most recent primal variable value. Dykstra’s algorithm is also used to solve optimization problems and, for many scenarios, is exactly equivalent to ADMM [9]. For EM, a likelihood maximization problem is solved by alternatively 1) finding the probability distribution of latent variables that maximizes the expected likelihood using fixed parameter values and 2) finding parameter values that maximize the expected likihood given the most recent latent variable distribution. ADMM, Dykstra’s algorithm and EM can all therefore be viewed as techniques that alternatively update different parameter subsets of a common function. In contrast, Algorithm (1) alternates between two different functions, $h()$ and $g()$ , that are applied to the same data before and after an invertable discrete transform. Thus, while Algorithm (1) is broadly similar to these techniques, it cannot be directly mapped to these methods. We are not aware of existing iterative algorithms that are equivalent to Algorithm (1).

3 Iterative convergence under sparsification

An interesting subclass of the general iterative method detailed in Algorithm 1 involves the use of sparsification functions for both $h()$ and $g()$ with $c()$ identifying convergence when a stable sparsity pattern is achieved in the output of $h()$ . Iterating between real domain and frequency domain sparsification is motivated by the discrete Fourier transform uncertainty principal [10], which constrains the total number of zero values in the real domain data and frequency domain representation generated by a discrete Fourier transform. Attempts to induce sparsity in one domain will reduce sparsity in the other domain with the implication that setting both $h()$ and $g()$ to sparsification functions will not simply result in the generation of a vector of all 0 elements. Instead, for scenarios where the iterative algorithm converages, the solution will represent a stable compromise between real and frequency domain sparsity.

We can explore the general convergence properties of Algorithm 1 using sparsification functions for $h()$ and $g()$ via the following simulation design:

•

Set $\mathbf{x}$ to a length $n$ vector of $\mathcal{N}(0,1)$ random variables.
•

Define $h()$ to generate a sparse version of $\mathbf{x}$ where the proportion $p$ of elements with the smallest absolute values are set to 0.
•

Define $g()$ to generate a sparse version of $\mathbf{w}$ where the proportion $p$ of complex coefficients with the smallest magnitudes are set to $0+0i$ .
•

Define $c()$ to identify convergence when the indices of 0 values in the output of $h()$ are identical on two sequential iterations.
•

The discrete and inverse discrete Fourier transforms are realized using the Fast Fourier Transform.

Refer to caption — Figure 1: Mean iterations until convergence for random length $n$ vectors of $\mathcal{N}(0,1)$ random variables and $h()$ and $g()$ functions that rank order the elements according to absolute value or magnitude and set the bottom 50% to 0. Convergence is based on repeating the same pattern of sparsity after execution of $h()$ on two sequential iterations.

Both general and sparsification versions of Algorithm 1 are supported by the IterativeFT R package available at https://hrfrost.host.dartmouth.edu/IterativeFT. Following this design, we applied the algorithm to 50 simulated $\mathbf{x}$ vectors for each distinct $n$ value in the range from 50 to 1,000 using the sparse proportion of $p=0.5$ and the maximum number of iterations $i_{m}=50$ . Figure 1 below displays the mean number of iterations until convergence as a function of $n$ . Figure 2 illustrates the results from a similar simulation that used a fixed $n$ of 500 and sparse proportion value ranging from $0.1$ to $0.9$ . For all of the tests visualized in Figures 1 and 2, the algorithm converged to a stable pattern of sparsity in $\mathbf{x}$ . Not surprisingly, the number of iterations required to achieve a stable sparsity pattern increased with the growth in either $n$ or $p$ . If $h()$ and $g()$ are changed to set all elements with absolute value or magnitude below the mean to 0 (see simulation design in Section 4), the relationship between mean iterations until convergence and $n$ is similar to that shown in Figure 1. Changing the generative model for $\mathbf{x}$ to include a non-random periodic signal (e.g., sinusodial signal or spike signal) also generates similar convergence results.

4 Detection of spike signals using iterative convergence

To assess the pratical utility of a sparsification version of Algorithm 1 for signal denoising, we applied the technique to data generated as either a one-dimensional vector $\mathbf{x}$ or two-dimensional matrix $\mathbf{X}$ containing a periodic spike signal with Guassian noise. For the one-dimensional cases, $h()$ , $g()$ , and $c()$ were specified as follows:

•

Define $h()$ to generate a sparse $\mathbf{x}$ where all elements $|x_{i}|\leq 1/n\sum_{j=1}^{n}|x_{j}|$ are set to 0.
•

Define $g()$ to generate a sparse $\mathbf{w}$ where all elements $|w_{i}|\leq 1/n\sum_{j=1}^{n}|w_{j}|$ are set to $0+0i$ (here the $||$ operation represents the magnitude of the complex number $w_{j}$ ).
•

Define $c()$ to identify convergence when the indices of 0 values in the output of $h()$ are identical on two sequential iterations.

For the matrix case, $h()$ , $g()$ , and $c()$ can be executed on a vectorized version of $\mathbf{X}$ . In the remainder of this paper, we will refer to the version of Algorithm 1 that uses these $h()$ , $g()$ , and $c()$ functions as the IterativeFT method.

4.1 Detection of spike signal in vector-valued input

For the one-dimensional case, the input vector $\mathbf{x}$ was generated as the combination of a periodic spike signal $\mathbf{s}$ and Gaussian noise $\boldsymbol{\varepsilon}$ , $\mathbf{x}=\mathbf{s}+\boldsymbol{\varepsilon}$ , with:

•

Elements $s_{i},i\in(1,...,n)$ of $\mathbf{s}$ set to 0 for all $i\neq a\lambda$ and generated as $U(\alpha_{min},\alpha_{max})$ random variables for $i=a\lambda,a\in(1,...,b)$ where $b$ is the total number of cycles and $\lambda\geq 1$ is the period.
•

Elements $\varepsilon_{i}$ of $\boldsymbol{\varepsilon}$ are generated as independent random variables with distribution $\mathcal{N}(0,\sigma^{2})$ .

Figure 3 shows an example of $\mathbf{x}$ generated according to this simulation model with $\alpha_{min}=\alpha_{max}=2.5,b=16,\lambda=8$ and $\sigma^{2}=0.5$ . For this specific example, the method converges in seven iterations and perfectly recovers the periodic spike signal.

To more broadly charaterize signal recovery for this simulation design, multiple $\mathbf{x}$ vectors were generated for different values of $\alpha_{min},\alpha_{max},b,\lambda$ , and $\sigma^{2}$ . Figure 4 shows the relationship between the MSE ratio achieved on convergence averaged across 50 simulated $\mathbf{x}$ vectors and the number of cycles, $b$ , captured in $\mathbf{x}$ . The MSE ratio is specifically computed as $MSE_{c}/MSE_{1}$ where $MSE_{c}$ is the mean squared error (MSE) between the output of the method after convergence and the spike signal $\mathbf{s}$ and $MSE_{1}$ is the MSE for the output from the first execution of $h()$ . For this simulation design, the average MSE ratio is very close to 0 for $b\geq 20$ , which reflects near perfect recovery of the input periodic spike signal.

Figure 5 captures the association between the average MSE ratio and the number of iterations completed by the algorithm ( $b$ was fixed at 20 for this simulation). These results demonstrate that signal recovery consistently improves on each iteration of the algorithm with the lowest MSE achieved upon convergence. Figure 6 captures the association between the average MSE ratio achieved on convergence and Gaussian noise variance ( $b$ was fixed at 20 for this simulation). These results demonstrate the expected increase in signal recovery error with increase noise variance and, importantly, show that the method still achieves improved noise recovery relative to just a single execution of $h()$ at high levels of noise.

4.2 Detection of spike signal in matrix-valued input

To explore the generalization where the input is a matrix rather than a vector, we also tested a simulation model where an input $n\times n$ matrix $\mathbf{X}$ is generated as follows:

•

Generate a length $n$ vector $\mathbf{s}$ that contains a periodic spike signal using the logic from Section 4.1 (see Figure 3 for visualization of a specific example generated according to this model).
•

Create the signal matrix $\mathbf{S}$ as $\mathbf{S}=\mathbf{s}\mathbf{s}^{T}$ .
•

Generate an $n\times n$ noise matrix $\mathcal{E}$ whose elements are independent $\mathcal{N}(0,\sigma^{2})$ random variables.
•

Create the input $n\times n$ matrix $\mathbf{X}$ as $\mathbf{X}=\mathbf{S}+\mathcal{E}$ .

Figure 7 shows an example of $\mathbf{X}$ generated according to this simulation model with $\alpha_{min}=\alpha_{max}=2.5,b=16,\lambda=8$ and $\sigma^{2}=0.5$ . For this specific example, the IterativeFT method converges in 24 iterations and perfectly recovers the periodic spike signal matrix. The convergence properties for this type of matrix input when evaluated across multiple simulations mirror those for the vector input as shown in Figures 4, 5, and 6.

4.3 Analysis of spike signal detection

The simulation-based results displayed in Sections 4.1 and 4.2 can be understood as follows:

•

The elements of the generated input vector $\mathbf{x}$ will be stochastically larger for elements that are non-zero in the periodic signal vector $\mathbf{s}$ given that the expected value of $\mathbf{x}$ is based solely on $\mathbf{s}$ :

$\displaystyle E[\mathbf{x}[j]]$	$\displaystyle=E[\mathbf{s}[j]+\varepsilon[j]]$
$\displaystyle E[\mathbf{x}[j]]$	$\displaystyle=E[\mathbf{s}[j]]+E[\varepsilon[j]]$	E of sum is sum of E
$\displaystyle E[\mathbf{x}[j]]$	$\displaystyle=E[\mathbf{s}[j]]$	$\varepsilon[j]\sim\mathcal{N}(0,1)$ so $E[\varepsilon[j]]=0$

•

These non-zero signal vector elements will therefore be less likely to be set to 0 by the sparsification version of $h()$ , which results in the frequency spectra of $h(\mathbf{x})$ being more strongly based on the frequency spectra of $\mathbf{s}$ than the spectra of $\boldsymbol{\varepsilon}$ .
•

These characteristics of the frequency spectra of the output from $h(\mathbf{x})$ mean that the output of $dft(h(\mathbf{x}))$ will be stochastically larger for elements that correspond to the spectra of $\mathbf{s}$ .
•

The sparsification version of $g()$ will therefore be more likely to retain frequency components associated with $\mathbf{s}$ and set to 0 frequency components corresponding to $\boldsymbol{\varepsilon}$ .
•

Repeated iterations will reinforce these patterns until a stable sparsity structure is obtained in $\mathbf{x}$ , which will tend to correspond to the very sparse frequency spectra of $\mathbf{s}$ .

To understand why the iterative application of both $h()$ and $g()$ on real and frequency representations of the data is needed one can consider scenarios where only one of $h()$ or $g()$ is performed (as outlined in Section 1.1, such scenarios are inherently non-iterative):

•

Only $h()$ is executed: Although a single execution of the sparsification version of $h()$ will tend to reduce the error between $\mathbf{s}$ and $\mathbf{x}$ , noise of a sufficient amplitude will result in the incorrect removal of some parts of $\mathbf{s}$ (i.e., negative noise may lower the value of a spike below the threshold) and incorrect retention of some parts of $\boldsymbol{\varepsilon}$ (i.e., the amplitude of the random noise may be larger than the threshold), e.g. Figure 3. The performance of a single $h()$ execution is assessed more comprehensively in 5.
•

Only $g()$ is executed: In this scenario, only a single execution of $dft()$ is performed on $\mathbf{x}$ , followed by execution of $g()$ before finally transforming back via $dft^{-1}()$ , i.e., the output is given by $dft^{-1}(g(dft(\mathbf{x})))$ . If one had prior knowledge of the spectra of $\mathbf{s}$ , then it would be possible specifically target the desired frequencies via $g()$ , e.g., a bandpass filter (though as shown in Section 5 below, even prior knowledge of an appropriate frequency band may not yield good denoising performance). Lacking this prior knowledge, however, sparsification via $g()$ can only be based the magnitude of frequency domain coefficients and, while it will tend to preserve the portions of the spectra of $\mathbf{x}$ associated with $\mathbf{s}$ , random noise will cause the incorrect removal of spectral components that are due to $\mathbf{s}$ and the incorrect retention of spectral components due to $\boldsymbol{\varepsilon}$ . The performance of a single $g()$ execution is assessed more comprehensively in 5.

5 Comparative evaluation of denoising performance

To evaluate the performance of the IterativeFT method relative to existing signal denoising techniques, a simulation study was undertaken on data generated according to several spike signal patterns with varying signal-to-noise ratios and spike frequencies.

5.1 Simulation design

The input data was generated as a length $n$ vector $\mathbf{x}=\mathbf{s}+\boldsymbol{\varepsilon}$ where the elements of $\boldsymbol{\varepsilon}$ were generated as independent $\mathcal{N}(0,\sigma^{2})$ random variables and the signal $\mathbf{s}$ followed one of three different periodic spike signal patterns:

1.

Uniform spike: Elements $s_{i},i\in(1,...,n)$ of $\mathbf{s}$ are set to 0 for all $i\neq a\lambda$ and set to the constant $\delta$ for $i=a\lambda,a\in(1,...,b)$ where $b=floor(n/\lambda)$ is the total number of cycles and $\lambda\geq 2$ is the period. The first panel in Figure 8 illustrates this signal pattern, which was also used for the example in Section 4.1 above.
2.

Alternating sign spike: The sparsity pattern is similar to the uniform spike model above but the non-zero values alternate between $\delta$ and $-\delta$ . The first panel in Figure 11 below illustrates this signal pattern.
3.

Alternating size spike: The sparsity pattern is similar to the uniform spike model above but the non-zero values alternate between $\delta/2$ and $\sqrt{7/4}\delta$ (these values yield a signal power equivalent to that for a signal with $\delta$ sized spikes). The first panel in Figure 14 below illustrates this signal pattern.

For the simulation results presented in Sections 5.3-5.5 below, $n=1024$ and $\delta=2$ . The noise variance, $\sigma^{2}$ , was set to yield a signal-to-noise ratio ranging between 0.05 and 2. The $n$ value of 1024 corresponds to a Nyquist frequency of 512 Hz and the spike period $\lambda$ was set to yield a ratio of spike frequency to Nyquist frequency ranging between 0.0625 (for $\lambda=32$ ) and 1 (for $\lambda=2$ ). For the single examples in Figures 8, 11 and 14, $n=128$ , $\lambda=8$ (i.e., a ratio of spike frequency to Nyquist frequency of 0.25) and $\sigma^{2}$ was set to yield a signal-to-noise ratio of 1.

5.2 Comparison methods

We compared the denoising performance of the IterativeFT method (Algorithm 1 using the $h()$ , $g()$ and $c()$ functions defined in Section 4) against four other techniques:

•

Real domain thresholding: This method performs a single thresholding operation on $\mathbf{x}$ using $h()$ .
•

Frequency domain thresholding: This method performs a single frequency domain thresholding operation using $g()$ , i.e., the denoised output is generated as $dft^{-1}(g(dft(x)))$ .
•

Butterworth bandpass filtering: This method is realized by a Butterworth [11] 3-order bandpass filter (as implemented by the butter() function in v0.3-5 of the gsignal R package [12]) where the band is set to the spike signal frequency (as a fraction of the Nyquist frequency) $\pm 0.05$ , e.g, if the spike signal is 0.25 of the Nyquist frequency, the band is set to 0.2 to 0.3.
•

Wavelet filtering: This method is based on the usage example for the threshold.wd() function in v4.7.4 of the wavethresh [13] R package. Specifically, a discrete wavelet transform is first applied to $\mathbf{x}$ using the Daubechies least-asymmetric orthonormal compactly supported wavelet with 10 vanishing moments and interval boundary conditions [14]. This wavelet transform is realized using the wd() function in the wavethresh R package with bc=”interval” and the default wavelet family and smoothness settings. The wavelet coefficients at all levels are then thresholded using a threshold value computed according to a universal policy and madmad deviance estimate on the finest coefficients followed by application of an inverse wavelet transformation.

5.3 Results for uniform spike model

Figures 8-10 illustrate the comparative denoising performance for the evaluated techniques on data simulated according to the uniform spike model. Figure 8 shows the outputs for a single example vector with a signal-to-noise (SNR) ratio of 1 and ratio of spike signal frequency to Nyquist frequency of 0.25. For this example, the IterativeFT method is able to perfectly recover the spike signal. In contrast, the other evaluated techniques have only marginal signal recovery performance.

Figure 9 shows the relationship between denoising performance and spike signal frequency at a signal-to-noise (SNR) value of 1. Denoising performance is plotted on the y axis as the mean MSE across 25 simulated inputs between the denoised data and the underlying signal and is plotted relative to the MSE measured on the undenoised data. A relative MSE value equal to 1 indicates null performance, i.e., the denoising method is equivalent to leaving the data unchanged; values above 1 indicate that the method further corrupts the signal. For this simulation model, the IterativeFT technique offers dramatically better denoising performance than the other four techniques with the exception of the high frequency domain where the Butterworth bandpass filter has slightly better performance. The simple real and frequency domain thresholding methods, i.e., a single application of either the $h()$ or $g()$ sparsification functions, are slightly better than null and relatively insensitive to signal frequency. The wavelet filter offers the good performance at low frequency values with performance decreasing towards the null level as the signal frequency increases. By contrast, the Butterworth filter is worse than null at low frequency values with performance steadily improving as frequency increases. The three notable dips in relative MSE for the IterativeFT method correspond to spike periods of 32, 16 and 8, which are all factors of 1024, the length of the input vector.

Figure 10 illustrates the relationship between denoising performance and the SNR value at a spike frequency that is 0.25 of the Nyquist frequency. Similar to Figure 9, both of the simple thresholding methods have performance that is only marginally better than null but, as expected, does improve slightly increasing SNR. The IterativeFT method is substantially better than all comparison techniques at SNR values above 0.5. At low SNR values, performance of the IterativeFT method decreases and falls below both the wavelet and Butterworth filtering techniques at SNR values below $\sim$ 0.25. The relative performance of the wavelet and Butterworth methods shows a steady decline with increasing SNR, which is surprising the error should decrease as SNR increases. This unexpected trend is due to the fact that the relative rather than absolute MSE is shown. In particular, the absolute MSE for these methods does decrease as the SNR increases, however, the MSE relative to no modifications increases with both filtering methods giving worse than null performance at SNR values greater than $\sim$ 1.1.

5.4 Results for alternating sign spike model

Figures 11-13 illustrate comparative denoising performance on data simulated according to the alternating sign spike model. Figure 11 shows the outputs for a single example vector with a signal-to-noise (SNR) ratio of 1 and ratio of spike signal frequency to Nyquist frequency of 0.25. Similar to the uniform spike example in Figure 8, the IterativeFT method perfectly recovers the spike signal with poor performance by the comparative techniques.

Figures 12 and 13 visualize the performance of the evaluated methods on the alternating sign model relative to signal frequency and SNR. The results are generally similar to those for the uniform spike model with the exception that the Butterworth bandpass filter has worse than null performance across all tested signal frequencies.

5.5 Results for alternating size spike model

Figures 14-16 illustrate comparative denoising performance on data simulated according to the alternating size spike model. Figure 14 shows the outputs for a single example vector with a signal-to-noise (SNR) ratio of 1 and ratio of spike signal frequency to Nyquist frequency of 0.25. Similar to the uniform spike and alternating sign examples, the IterativeFT method provides excellent signal recovery though it does underestimate the magnitude of the smaller spike; the other evaluated methods are worse by comparison.

Figures 15 and 16 visualize the performance of the evaluated methods on the alternating size model relative to signal frequency and SNR. The results are generally similar to those for the uniform spike model with the exception that performance for the Butterworth bandpass filter is never better than the IterativeFT method.

6 Conclusion

In this paper we explored a family of iterative algorithms based on the repeated execution of discrete and inverse discrete Fourier transforms on real valued vector or matrix data. One interesting member of this family, which we refer to as the IterativeFT method, is motivated by the discrete Fourier transform uncertainty principle and involves the application of a thresholding operation to both the real domain and frequency domain data with convergence obtained when real domain sparsity hits a stable pattern. As we demonstrated through simulation studies, the IterativeFT method can effectively recover periodic spike signals across a wide range of spike signal frequencies and signal-to-noise ratios. Importantly, the performance of the IterativeFT method is significantly better than standard non-iterative denoising techniques including real and frequency domain thresholding, wavelet filtering and Butterworth bandpass filtering. An R package implementing this technique and related resources can be found at https://hrfrost.host.dartmouth.edu/IterativeFT. Areas for future work include on this method include:

•

Expanding on the discussion in Section 4.3 to more thoroughly explore the theoretical basis for the denoising performance shown in Sections 4 and 5.
•

Expanding the comparative evaluation to include other denoising techniques, e.g., basis pursuit [15], and a broader range of periodic signal and noise models, e.g., composition of spike signals at multiple frequencies/amplitudes, mixtures of harmonic and spike signals, etc.
•

Exploring other classes of $h()$ , $g()$ and $c()$ functions and associated analysis problems, e.g., soft thresholding.
•

Exploring the generalization of the IterativeFT method to complex or hypercomplex-valued inputs.
•

Exploring the generalization of the IterativeFT method to other invertable discrete transforms, e.g., discrete wavelet transform [4].

Acknowledgments

This work was funded by National Institutes of Health grants R35GM146586 and R21CA253408. We would like to thank Anne Gelb and Aditya Viswanathan for the helpful discussion.

References

[1] Ronald N Bracewell. The Fourier transform and its applications. McGraw Hill, Boston, 3rd ed edition, 2000.
[2] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn., 3(1):1–122, jan 2011.
[3] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39(1):1–38, 1977.
[4] Stphane Mallat. A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way. Academic Press, Inc., USA, 3rd edition, 2008.
[5] D Sundararajan. Fourier Analysis-A Signal Processing Approach. 1st ed. 2018 edition.
[6] Ann Maria John, Kiran Khanna, Ritika R Prasad, and Lakshmi G Pillai. A review on application of fourier transform in image restoration. In 2020 Fourth International Conference on I-SMAC (IoT in Social, Mobile, Analytics and Cloud) (I-SMAC), pages 389–397, 2020.
[7] Lu Chi, Borui Jiang, and Yadong Mu. Fast fourier convolution. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 4479–4488. Curran Associates, Inc., 2020.
[8] Richard L. Dykstra. An algorithm for restricted least squares regression. Journal of the American Statistical Association, 78(384):837–842, 1983.
[9] Ryan J. Tibshirani. Dykstra’s algorithm, admm, and coordinate descent: Connections, insights, and extensions. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, pages 517–528, Red Hook, NY, USA, 2017. Curran Associates Inc.
[10] David L. Donoho and Philip B. Stark. Uncertainty principles and signal recovery. SIAM Journal on Applied Mathematics, 49(3):906–931, 1989.
[11] S. Butterworth. On the Theory of Filter Amplifiers. Experimental Wireless & the Wireless Engineer, 7:536–541, October 1930.
[12] Van Boxtel, G.J.M., et al. gsignal: Signal processing, r package version 0.3-5 edition, 2021.
[13] Guy Nason. wavethresh: Wavelets Statistics and Transforms, 2022. R package version 4.7.2.
[14] Albert Cohen, Ingrid Daubechies, and Pierre Vial. Wavelets on the interval and fast wavelet transforms. Applied and Computational Harmonic Analysis, 1(1):54–81, 1993.
[15] Shaobing Chen and D. Donoho. Basis pursuit. In Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers, volume 1, pages 41–44 vol.1, 1994.