The SWDFT
The SWDFT extends the discrete Fourier transform (DFT) from the frequency domain to the time-frequency domain. Whereas the DFT makes a single large computation on the entire signal, the SWDFT takes a sequence of smaller DFTs on temporally ordered subsets of the signal. To be clear, this vignette refers to data as a signal, but you can easily substitute signal with time-series, data, input, etc. In any case, if x is a length N signal, then the DFT of x computes N complex-valued coefficients, and is defined by:
\[\begin{eqnarray} a_k &=& \sum_{j=0}^{N-1} x_j \omega_n^{-jk} \nonumber \\ k &=& 0, 1, \ldots, N - 1, \end{eqnarray}\]where
\[\begin{eqnarray} \omega_N &=& \cos(\frac{2 \pi}{N}) + i \sin(\frac{2 \pi}{N}). \end{eqnarray}\]Each DFT coefficient measures the correlation between x and a sine and cosine function with a particular frequency, where frequency specifies how fast the sine and cosine functions oscillate. This explains why researchers are interested in large DFT coefficients: large DFT coefficients indicate high correlation between our signal a sine or cosine function with a particular frequency.
A pictorial example of the DFT for a length N = 8 signal is given below, where the blue squares are the input signal, and the red circles are the DFT coefficients:
In contrast, the SWDFT of takes a sequence of smaller DFTs, and outputs a 2D complex array coefficients. In this 2D array, the first dimension represents frequency, and the second dimension represents time, which is why the SWDFT is a time-frequency method. For example, the SWDFT for the same length N = 8 signal used in the previous signal is:
The size of the smaller DFTs is known as the window size, so in this example, the window size is four. To see how the SWDFT works, let’s walk through the computations in detail. Since our window size is length four, step one computes a length four DFT, and stores the output in the first window position of the 2D array:
Step two moves the window one position to the right, computes the DFT of the second through fifth data-points, and stores the results in the second window position:
This process continues until the window position reaches the final data-point:
And that’s it. Formally, the SWDFT for a length N signal with a length n window size is defined by:
\[\begin{eqnarray} a_{k, p} &=& \sum_{j=0}^{n-1} x_{p-n+1+j} \omega_{n}^{-jk} \nonumber \\ k &=& 0, 1, \ldots, N - 1 \nonumber \\ p &=& n - 1, n, \ldots, N - 1. \end{eqnarray}\]In this notation, The k subscript in \(a_{k, p}\) indexes frequency, and the p subscript indexes window position, where window position often represents time. In this package, the swdft function implements the SWDFT. For example, you can take the SWDFT of a length 100 random Gaussian noise series with:
N <- 100
window_size <- 32
x_random <- rnorm(n=N, mean=0, sd=1)
a <- swdft::swdft(x=x_random, n=window_size)But the SWDFT is most useful for signals with frequency domain characteristics that change over time. For this reason, this vignette demonstrates the how this package works using local cosine signals as data. The next section defines a local cosine signal, and shows how you can generate them with the local_signal function.
Local Cosine Signals
The sliding window discrete Fourier transform (SWDFT) is useful for signals frequency characteristics that change over time. Therefore, the signals we use to demonstrate how the swdft package works should have changing frequency characteristics. To accomplish this, we use local cosine signals as data. A local cosine signal is simply a cosine function multiplied by an indicator function, where the indicator function is only nonzero for a contiguous subset of the signal. Formally, a length N local cosine signal is defined by:
The parameters of a local cosine signal are:
- A: Amplitude of the cosine function
- f: Frequency of the cosine function
- \(\phi\): Phase of the cosine function
- S: Start of the local cosine signal
- L: Length of the local cosine signal
You can generate a local cosine signal with the local_signal function. This example also adds a small Gaussian noise factor to the local cosine signal:
# --- Generate a local cosine signal + Gaussian noise ---
set.seed(999)
## Specify the parameters for the local cosine signal
signal_length <- 96 # The length of the local cosine signal
freq <- 3 / 32 # Frequency, interpreted as "3 complete cycles every 32 data-points"
amplitude <- 1 # The amplitude of the local cosine signal
phase <- 0 # Phase of the local cosine signal
periodic_start <- 20 # When the local cosine signal starts
periodic_length <- 50 # How long the local cosine signal lasts
signal <- swdft::local_signal(N=signal_length, A=amplitude,
Fr=freq, phase=phase, S=periodic_start,
L=periodic_length)
## Generate the Gaussian noise factor, and add to the local cosine signal
sigma <- .3 # Standard deviation of the Gaussian noise factor
noise <- rnorm(n=signal_length, mean=0, sd=sigma)
x <- signal + noise
## Plot what the local cosine signal + noise looks like
plot(x, pch=19, cex=1.4, xlab="", ylab="Signal",
main="Local Cosine Signal plus Gaussian Noise")
lines(signal, lwd=2, col="red")
In this figure, the red line is the local cosine signal, and the black dots are the local cosine signal plus noise.
The swdft function
The cornerstone of this package is the swdft function. The swdft function is analogous to base R’s stats::fft function for the time-frequency case. Just as the fft function returns a 1D complex array of DFT coefficients, the swdft function returns a 2D complex array of SWDFT coefficients. You can take the SWDFT of the local cosine signal (x) generated in the previous section with the following code:
window_size <- 32
a <- swdft::swdft(x=x, n=window_size)An important thing to k now is that by default, we pad the input signal (x) with zeros, so that the dimensions of the SWDFT match the dimensions of the input signal:
length(x)
#> [1] 96
dim(a$a)
#> [1] 32 96If you don’t want to pad the signal, you can use the pad=FALSE option:
a_nopad <- swdft::swdft(x=x, n=window_size, pad=FALSE)
dim(a_nopad$a)
#> [1] 32 65To visualize the SWDFT, because SWDFT coefficients are complex-valued, it is common to display the squared modulus of the SWDFT coefficients. For example, if \(a_{k, p}\) is the SWDFT coefficient at frequency k and window position p, the squared modulus of \(a_{k, p}\) is given by:
\[\begin{eqnarray} Mod(a_{k, p})^2 &=& Re(a_{k, p})^2 + Im(a_{k, p})^2. \end{eqnarray}\]In R, you can compute the squared modulus of a complex number with the Mod function. For example, you can compute the squared modulus of the SWDFT coefficients with:
amod <- Mod(a$a)^2The implementation of the swdft function simply takes a fast Fourier transform (FFT) in each window position, and can be easily understood by looking at the code of the internal swdft_fft function:
swdft_fft <- function(x, n, taper) {
N <- length(x)
P <- N - n + 1
a <- array(data = NA, dim = c(n, P))
for (p in n:N) {
a[, p - n + 1] <- stats::fft(z = x[(p - n + 1):p] * taper)
}
return(a)
}The actual implementation of the swdft function is a wrapper, which calls either swdft_fftw or swdft_fft, depending on whether the fftwtools library is installed. The fftwtools package is an R-wrapper to the Fastest Fourier Transform in the West (FFTW) library (Frigo and Johnson 2005), so swdft_fftw is faster than swdft_fft.
The swdft function returns an S3 object of class swdft. This class includes the 2D complex array of SWDFT coefficients (a$a), along with basic information about the parameters used (see ?new_swdft for more information). You can visualize the squared modulus of the SWDFT coefficients with the plot.swdft method:
plot(a)
The y-axis of this plot represents frequency, and the x-axis represents window position, and the x-axis typically represents time. In the next section, we describe several options to customize plot.swdft visualizations.
The plot.swdft method
The plot.swdft method visualizes objects created with the swdft function. This section walks through the most frequently used options to customize these visualizations.
The default color used by plot.swdft is grayscale. However, I often use the blue-to-red colorscale tim.colors, which comes from the fields package. You can use the tim.colors colorscale with the option col="tim.colors":
plot(a, col="tim.colors", title="SWDFT using the 'tim.colors' colorscale")
Often, we need to display the frequency dimension (i.e. the y-axis) in different units, such as Hertz. For a SWDFT with length n windows, the \(k^{th}\) frequency can be interpreted as “k complete cycles in a length n period of time”. Hence, the default frequency unit option we use is freq_type="cycles", which indexes the frequency dimension with \(k = 0, 1, \ldots, n - 1\). Another popular frequency unit is Hertz, which can be used with the freq_type="hertz" option. If you use Hertz, you also need to specify the sampling frequency with the fs= option, for example:
sampling_frequency <- 100
plot(a, freq_type="hertz", fs=sampling_frequency)
The final frequency unit option we provide is freq_type="fraction". In many fields, it is common to represent frequency as the decimal value \(\frac{k}{n}\), for \(k = 0, 1, \ldots, n-1\), rather than the integer k used in the freq_type="cycles" option. The freq_type="fraction" option is used in base R’s plot.spec method for the spectrum function. An example of this option is given by:
plot(a, freq_type="fraction")
When freq_type="fraction", we follow convention by only displaying the frequency range 0.0 to 0.5. This is a convention because for real-valued signals, the Fourier coefficients larger than 0.5 are equivalent to the complex conjugate of a Fourier coefficient less than 0.5. This implies that that squared modulus values are identical. This fact explains why, in the previous plots, the frequencies flipped across the middle of the y-axis are identical. Technically speaking, this phenomenon is called “aliasing”, and a good reference to understand aliasing is Section 2.5 of (Bloomfield 2004).
Finally, you sometimes need to change either the axis labels or values, for instance, you may want the x-axis to represent time. To change the x-axis values, you can use the custom_xaxis= option (there is a corresponding custom_yaxis as well). If you use this option, you need to make sure that the length of the custom axis is the same as the dimension of the swdft object, which you can check in R with dim(a$a). For example, if I wanted the x-axis of the local cosine signal to represent years from 1990 to 1995, I would use:
years <- 1900:1995
plot(a, freq_type="fraction", col="tim.colors", custom_xaxis=years, xlab="Years")
Those are the most frequently used options. For details and more options, see the ?plot.swdft documentation.
Tapering and Smoothing
The SWDFT is essentially a sequence of discrete Fourier transforms (DFTs), and DFT coefficients are notoriously noisy. By noisy, I mean that as the signal gets larger, the variance stays the same, which implies that DFT coefficients are not a consistent estimator of the spectral density function. If the last sentence didn’t make sense, don’t worry, it’s a complicated topic, and a good reference to learn more is Chapter 18 of (Kass, Eden, and Brown 2014). For our purposes, the takeaway is that DFT coefficients aren’t perfect from a statistical perspective, and it’s possible to improve them. Statisticians have improved DFT coefficients in two ways: tapering and smoothing. Tapering applies a window to each end of the signal before taking the DFT, and smoothing smooths the final DFT coefficients, for example with a moving average. The swdft function provides options for both tapering (Section 18.3.6 of (Kass, Eden, and Brown 2014)) and smoothing (Section 18.3.3 of (Kass, Eden, and Brown 2014)).
To taper the SWDFT coefficients, we provide the widely used cosine bell taper, based on R’s spec.taper function. You can use this taper with the taper_type="cosine" option. In addition, you can specify the proportion of each end of the signal to taper with the p= option. An example of is given by:
a_taper <- swdft::swdft(x=x, n=window_size, taper_type="cosine", p=.2)
plot(a_taper)
To smooth the SWDFT coefficients, we provide the smooth= option. This option builds a smoothing kernel using base R’s stats::kernel function. The additional option m= controls the width of the kernel, and the option num_convs= specifies the number of times the kernel is convolved with itself. See ?stats::kernel for more details. An example of smoothing the SWDFT coefficients is given by
a_smooth <- swdft::swdft(x=x, n=window_size, smooth="daniell", m=1, num_convs=0)
plot(a_smooth)
It is important to note that the most widely used tapering approach combines multiple tapers, which is often referred to as the Multitaper method (Thomson 1982). The swdft package does not provide an implementation of the multitaper method, but I would be willing to if there’s significant interest. For an implementation, see the excellent multitaper R-package.
