Abstract

In estimating the frequency spectrum of real-world time series data, we must violate the assumption of infinite-length, orthogonal components in the Fourier basis. While it is widely known that care must be taken with discretely sampled data to avoid aliasing of high frequencies, less attention is given to the influence of low frequencies with period below the sampling time window. Here, we derive an analytic expression for the side-lobe attenuation of signal components in the frequency domain representation. This expression allows us to detail the influence of individual frequency components throughout the spectrum. The first consequence is that the presence of low-frequency components introduces a 1/f component across the power spectrum, with a scaling exponent of . This scaling artifact could be composed of diffuse low-frequency components, which can render it difficult to detect a priori. Further, treatment of the signal with standard digital signal processing techniques cannot easily remove this scaling component. While several theoretical models have been introduced to explain the ubiquitous 1/f scaling component in neuroscientific data, we conjecture here that some experimental observations could be the result of such data analysis procedures.

1  Introduction

The Fourier transform of a function has many applications in physics and engineering, and each is referred to as frequency domain and time-domain representation, respectively (Fourier, 1822). The power spectral density (PSD) is a powerful tool in the analysis of signals and distinguishing systems. However, the main assumption of a Fourier series decomposition, infinitely long data, is violated, since all experimental data are inherently of finite length.

The discrete Fourier transform (DFT) is used to analyze finite-length data (Cooley, Lewis, & Welch, 1969),
formula
1.1
and is related to the continuous-time Fourier transform (Fourier, 1822; Titchmarsh, 1948),
formula
1.2
by considering a discrete time function where is replaced by with and where the sampling frequency is . The sampling theorem of Nyquist (1928) and Küpfmüller (1928) specifies that the sampling rate must be greater than twice the highest-frequency content of the signal. In contrast, however, the effect of low-frequency content has not been studied extensively.
Here, we investigate the case where spectral content exists below the low-frequency bound of the DFT. We find that the presence of these low frequencies has a profound effect on the frequency spectrum. Considering that one period () is necessary for the lowest-frequency content during the sampling time window , where is the number of samples, then the relevant cutoff frequency is
formula
1.3
which is exactly the DFT bin length (Harris, 1978). In his classic 1978 paper, Harris introduced spectral leakage and the effect of different windowing methods to suppress power in the side lobes. Since then, several papers have been published on the design of windows with less side lobe power (see, e.g., Nuttall, 1981, among others). Here an analytical form of these spectral components and the envelope of the DFT spectrum are presented and used to explain some scaling artifacts.
The significant effect of frequency content below the cutoff frequency is illustrated in Figure 1. A signal with a noise-free, 10 Hz oscillating component, observed over a time window of 100 seconds, produces a given spectral estimate (black); however, when an additional, sub-cutoff test frequency (see equation 1.3) component is added to the signal, a roughly linear feature appears across the frequency spectrum in log-log coordinates (see Figure 1a). This feature is surprising in a noise-free signal. Further, when a linear drift (which can be seen as equivalent to an extremely low-frequency component) is added to the signal, a similar linear feature appears (see Figure 1b). While the observations of these effects have been noted previously in the neuroscience literature (Mitra & Bokil, 2007), we present here a systematic analytical derivation of this effect as a result of finite window length. It is important to note that in the numerical spectral estimates, as in the analytic results to follow, we work with the modulus (magnitude) of the Fourier transform. In the literature (Ward & Greenwood, 2007), 1/f noise is referred to as the spectral density of a stochastic process, which has the form .
Figure 1:

Effect of sub-cutoff frequency content (a) and linear drift (b) on spectral scaling. The fast Fourier transform (FFT) magnitude spectral estimate for a signal , where Hz, Hz, and s, is given in black. (a) An additional signal ( Hz), which is an arbitrary frequency smaller than the cutoff (in equation 1.3), is added to . (b) A linear drift is added to . In both cases, the resulting spectral estimate is shown in magenta.

Figure 1:

Effect of sub-cutoff frequency content (a) and linear drift (b) on spectral scaling. The fast Fourier transform (FFT) magnitude spectral estimate for a signal , where Hz, Hz, and s, is given in black. (a) An additional signal ( Hz), which is an arbitrary frequency smaller than the cutoff (in equation 1.3), is added to . (b) A linear drift is added to . In both cases, the resulting spectral estimate is shown in magenta.

We have investigated this low-frequency-induced artifact present in spectral analysis of finite-length signals. This study demonstrates that although the effect becomes dominant for test frequencies below (see equation 1.3), small artifacts are present throughout the frequency spectrum for all test frequencies. Finally, we demonstrate that these effects are robust to removal with standard filtering techniques (see the appendix). While theoretical models have been proposed to explain the ubiquitous scaling observed in neuroscientific data, for example as a result of neuronal shot noise (Milstein, Mormann, Fried, & Koch, 2009), we propose that the effect studied here is important for the analysis of neuroscientific data and that care should be taken with interpreting power spectral estimates.

In summary, this letter has three main points:

  1. Fourier analysis is based on the assumption of infinite window length. This assumption is always necessarily violated in real-world data analysis.

  2. Consequences of this violation include nonlinearities resulting from multiple frequencies in the data, spectral leakage, and scaling artifacts induced by low-frequency components in the data (below ).

  3. The interplay of frequencies and phase is often neglected in the analysis of real-world data.

The letter is organized as follows. In section 2, the analytical forms of the continuous-time Fourier transform for infinite and finite windows are derived. From the finite-length form, we derive the analytical form of the spectrum and the analytical forms of its piecewise envelopes. We make the connection to spectral leakage in section 2.1.1, the observed low-frequency induced scaling artifact in section 2.3, and nonlinearities in the spectrum when dealing with multiple frequencies in the data in section 2.1.3. We perform the analyses following Goertzel (1958) and an alternative approach. In section 2.3, we present an analytical derivation of the 1/f scaling artifact and briefly address the interplay between frequency and phase, which will be the subject of future research. Section 3 presents concluding remarks and provides perspectives on how this spectral scaling artifact could be an important consideration for both theoreticians and experimentalists using spectral analysis. The appendix presents rigorous numerical evaluations of the robustness of this artifact to common signal processing techniques, such as windowing and prefiltering.

2  Analytical Results

Previously, we introduced a mathematical approach to time-domain spectral analysis (Lainscsek & Sejnowski, 2015), which is a generalization of the Goertzel algorithm (Goertzel, 1958), a well-known technique for efficient numerical evaluation of individual terms in the DFT.

2.1  Goertzel Algorithm

The Goertzel algorithm uses equation 1.2 in the one-sided form,
formula
2.1
to evaluate the signal spectrum at given individual frequencies . Starting from a signal , in the standard Goertzel approach,
formula
2.2
is then computed from the real part and the imaginary part using the expectation operator and by
formula
2.3

2.1.1  Single-frequency signal

For a single-frequency signal with infinitely dense sampling, the real part and the imaginary part for a finite-length window are then
formula
2.4
where the spectral leakage terms are
formula
2.5
For , this is
formula
2.6
but for finite window length, there are additional error terms,
formula
2.7
with the expressions for and in equations 2.4. The expressions for using equations 2.5 are
formula
2.8

2.1.2  Time Shift Theorem for a Single Frequency

Concerning the Fourier transform of a single-frequency signal with a time shift , it is well known (see, e.g., Erdélyi, 1954) that
formula
2.9
Thus, for an arbitrary time shift, only a phase shift is added to the spectral representation. The substitution in equation 2.9 is possible only for infinitely long data. This theorem therefore holds only in the limit of infinite window length; however, with finite windows, the effect of a time shift on the error terms (see equations 2.7) must be considered. For a single frequency signal and finite window length, equation 2.9 changes to
formula
2.10

2.1.3  Multiple Frequency Signal

For a signal with frequencies and phases, the Fourier transform for finite-length windows can be composed of the real and imaginary parts,
formula
2.11
where and are the expressions from equations 2.4 and 2.5, respectively. The Fourier transform
formula
2.12
will then have cross-terms resulting from the square of the expressions in equation 2.11, which is a nonlinear function of the frequencies and phases. This means that in addition to the linear superposition of Fourier transformations of the signal components, there are error terms that result from nonlinear spectral leakage in the case of finite window length. Finally, as already noted, these nonlinear error terms must be considered in the context of the time-shift theorem in the case of multiple frequency signals.

2.2  Alternative Representation

If we replace the two real and imaginary probe signals in the above representation with a single probe signal with frequency and phase , we get
formula
2.13
where is the value of the probing phase that maximizes (Lainscsek & Sejnowski, 2015).

2.2.1  Single Frequency

For a single frequency, this calculation then yields
formula
2.14
which is the time domain spectrum for signals of infinite length. Generalizing this next to finite-length signals with a single frequency ,
formula
2.15a
formula
2.15b
Again, it is apparent that each individual frequency has an effect throughout the spectrum, as shown above and in Figure 2. The red lines in Figure 2 depict all possible phases in the probe signal . Although these two approaches are not completely equivalent, the functional representation contained in equations 2.15 is much simpler. Figure 3 illustrates the error made in both approaches and the small but nonnegligible difference between them. Light blue lines indicate the error made at each phase relative to the window and dark blue is the sinc function. Importantly, the error made in both approaches is very similar at all window lengths.
Figure 2:

Analytical form of for test frequencies with several random phases and (red; see equation 2.15b) and derived envelopes (see equations 2.16). The width of and together is exactly the DFT bin length and is centered in the main lobe.

Figure 2:

Analytical form of for test frequencies with several random phases and (red; see equation 2.15b) and derived envelopes (see equations 2.16). The width of and together is exactly the DFT bin length and is centered in the main lobe.

Figure 3:

Error made for a single sinusoidal signal with a frequency of 10 Hz and sampling rate of 1000 Hz in the alternative representation (top) and for the Goertzel algorithm (middle). The difference (bottom) between both approaches, in percent of the signal amplitude, is small but nonvanishing. The light blue lines are the errors made for the different phases (; ) in equation 2.15a, and the dark blue lines represent the sinc function (; ).

Figure 3:

Error made for a single sinusoidal signal with a frequency of 10 Hz and sampling rate of 1000 Hz in the alternative representation (top) and for the Goertzel algorithm (middle). The difference (bottom) between both approaches, in percent of the signal amplitude, is small but nonvanishing. The light blue lines are the errors made for the different phases (; ) in equation 2.15a, and the dark blue lines represent the sinc function (; ).

2.2.2  Envelopes

Outside the frequency , the envelope of with test frequency is
formula
2.16a
formula
2.16b
formula
2.16c
formula
2.16d
where and . The main lobe spans the range between and and is surrounded by the side lobes (see Figure 2). The lengths of and are each half the DFT bin width. Equations 2.16 hold only if is an integer multiple of . If , the side lobes are not equidistant and the envelopes cannot be expressed in the form of equations 2.16. It is not possible to express the envelope in general if there are multiple frequencies in the signal. Then a linear superposition of the envelopes of equations 2.16 for each of the frequencies in the signal will lie above the spectrum. Further analytical investigations will be the subject of future research.

2.3  Analytical Derivation of the 1/f Scaling Artifact

For a frequency lower than , only the term contributes (for ):
formula
2.17
The magnitude of the spectrum then scales as with a small correction term. It is important to note again that beause we are working in magnitude of the Fourier transform, the power spectrum will then scale with .

As a direct consequence of equations 2.15b, 2.16b, and 2.16c for two frequencies greater than , where is smaller than , the signal phases (both absolute and relative) determine their discriminability in the resulting spectrum. This effect will be the subject of future study.

We can see in equation 2.15b that there is a phase dependence that results in slopes ranging between 1 and 2. The slope is 2 only in the case . In Figure 4, we tested the numerical Goertzel algorithm against the analytical expression in equation 2.15b and found a nearly perfect fit.
Figure 4:

Effect of the phase of the low-frequency component ( Hz, Hz, and s) relative to the finite-length window. is shown for phases between 0 and as separate linear plots (top) and as polar plots where the logarithm of the frequency is plotted along the radius (bottom). The analytical form in equation 2.15b is shown on the left side and the numerical results for the Goertzel algorithm on the right side. The slope is 1 nearly everywhere. The slope is 2 only for . The dashed black lines (top plots) are the analytical form in equations 2.16.

Figure 4:

Effect of the phase of the low-frequency component ( Hz, Hz, and s) relative to the finite-length window. is shown for phases between 0 and as separate linear plots (top) and as polar plots where the logarithm of the frequency is plotted along the radius (bottom). The analytical form in equation 2.15b is shown on the left side and the numerical results for the Goertzel algorithm on the right side. The slope is 1 nearly everywhere. The slope is 2 only for . The dashed black lines (top plots) are the analytical form in equations 2.16.

3  Conclusion

We have presented analytical derivations illustrating the nonlinear effects present in spectral analysis. In the appendix, we support our results with numerical simulations. These aspects are a direct consequence of violating the assumption of infinite-length sinusoidal components in Fourier analysis. These nonlinear effects result from couplings of the spectral leakage contributions from individual frequency components and can become profound in certain instances. We applied the theoretical results to the case of low-frequency components in a noise-free signal, where a specific low-frequency-induced spectral scaling artifact can occur.

To summarize, following are the main points explored in this work:

  • Consequences of analyzing finite-length signals in Fourier analysis

    • Nonlinearities resulting from multiple frequencies in the data. For real-world, finite-length data, the limit of is not realized, which causes nonlinear interactions between single-frequency components.

    • Spectral leakage. We analyzed the spectral leakage effects introduced by Harris (1978) in an analytical manner.

    • Scaling artifacts induced by sub-cutoff frequencies (below ). Scaling artifacts can be caused by special spectral leakage effects.

  • Interplay of frequencies and phase is often neglected in analysis of real-world data. Due to nonlinearities caused by finite-length windows, frequencies and phases are interconnected in a nonlinear way. Consequences of this are the subject of ongoing research.

Recognition of the possibility of this artifact and nonlinear interconnections between phase and frequency is important for neuroscientists, physicists, and engineers who use spectral analysis in their work. As illustrated in the appendix, this artifact is robust to removal through common linear filtering and detrending techniques, and thus may be an important practical consideration for analysis of experimental data. More robust filtering approaches and further consequences of nonlinearities in the signal processing technique, will be the subject of future work.

Appendix A: Numerical Robustness Tests

We numerically tested the robustness of this spectral scaling artifact under various common spectral analysis situations. We began with the dependence of the spectral scaling artifact on the amplitude of the added sub-cutoff frequency component. Starting with a sinusoid of unit amplitude, we added low-frequency components with amplitudes between 0 and 1; we show the resulting frequency spectra in Figure 5. The scaling of the resulting spectra is invariant to the amplitude of the sub-cutoff component, consistent with the slope in equation 2.17. It is important to note that we verified the robustness of this spectral scaling bias using many different spectral estimation techniques (Welch power spectral estimate, Goertzel algorithm, Hann and Hamming windows, multitaper spectral estimate; data not shown but provided as online supplementary code).
Figure 5:

Effect of the amplitude of the low-frequency component on the spectral scaling. FFT estimates of the frequency spectrum for varying amplitude of the 0.009 Hz (0.9) sub-cutoff frequency component. Amplitude of the sub-cutoff frequency added to the original 10 Hz sinusoid ranges from 0 (dark red) to 1 (blue); the original sinusoid is of unit amplitude. Note that the specific scaling component analyzed in the main text appears immediately on the addition of the low-frequency component at all amplitudes.

Figure 5:

Effect of the amplitude of the low-frequency component on the spectral scaling. FFT estimates of the frequency spectrum for varying amplitude of the 0.009 Hz (0.9) sub-cutoff frequency component. Amplitude of the sub-cutoff frequency added to the original 10 Hz sinusoid ranges from 0 (dark red) to 1 (blue); the original sinusoid is of unit amplitude. Note that the specific scaling component analyzed in the main text appears immediately on the addition of the low-frequency component at all amplitudes.

Further, these spectral artifacts are difficult to remove using standard filtering techniques. We tested the efficacy of both finite impulse response (FIR) and infinite impulse response (IIR) high-pass digital filters in recovering the original signal spectrum. Filtering was carried out in the forward direction only. Linear-phase FIR filters ranging from 100th to 5000th order were constructed using a Hamming window (cf. Matlab fir1). The numerical stability of the FIR approach allows construction of a filter at high order, which is necessary to achieve sufficient magnitude response in the low-frequency stopband. Butterworth IIR filters (cf. Matlab butter) were constructed at fourth order, displaying a slightly stronger magnitude response. In neither case were the filters generally able to recover the original signal spectrum above the cutoff frequency (0.5 Hz); the spectral scaling component remained (see Figure 6 and supplementary Matlab code for replication). The high-order FIR filters were successful, however, in one specific case. Interestingly, the success of the FIR filtering operation depends critically on the phase angle of the low-frequency component, and in a specific manner: the original signal spectrum is recovered only if the imaginary part of the Fourier transform of the low-frequency component is nearly zero, that is, if the phase angle of the low-frequency component is zero or (see Figure 7). Such a nontrivial dependence of the filtering operation on the phase of signal components is hitherto unknown, to the best of our knowledge, and will be the subject of future study.
Figure 6:

Effect of high-pass prefiltering on low-frequency intrusions. As in Figure 1, the FFT spectral estimates are given for an example 10 Hz sinusoid (black) and for the original signal with an added 0.009 Hz low-frequency component with random phase angle (magenta). A version of this composite signal filtered with a high-pass, 5000th-order, linear-phase FIR filter is plotted in turquoise.

Figure 6:

Effect of high-pass prefiltering on low-frequency intrusions. As in Figure 1, the FFT spectral estimates are given for an example 10 Hz sinusoid (black) and for the original signal with an added 0.009 Hz low-frequency component with random phase angle (magenta). A version of this composite signal filtered with a high-pass, 5000th-order, linear-phase FIR filter is plotted in turquoise.

Figure 7:

Dependence of filtering on the phase of the sub-cutoff frequency component. The magnitude of slope for the spectral scaling is plotted as a polar function of the low-frequency phase for the original signal with an added 0.009 Hz low-frequency component (magenta) and for the signal filtered as in Figure 6 (turquoise). For this calculation, mean magnitude spectra were obtained by averaging over 100 trials at each phase using the Welch power spectral estimate, and fits were conducted on the smoothed results. The unit circle is plotted in black. Slope is determined by linear fit in log-log coordinates between 0.7 and 3 Hz, where the spectral scaling artifact is most prominent in the filtered version (see Figure 6). While the slope remains near 1 for the unfiltered signal, the slope of the filtered signal returns to zero only for angles close to .

Figure 7:

Dependence of filtering on the phase of the sub-cutoff frequency component. The magnitude of slope for the spectral scaling is plotted as a polar function of the low-frequency phase for the original signal with an added 0.009 Hz low-frequency component (magenta) and for the signal filtered as in Figure 6 (turquoise). For this calculation, mean magnitude spectra were obtained by averaging over 100 trials at each phase using the Welch power spectral estimate, and fits were conducted on the smoothed results. The unit circle is plotted in black. Slope is determined by linear fit in log-log coordinates between 0.7 and 3 Hz, where the spectral scaling artifact is most prominent in the filtered version (see Figure 6). While the slope remains near 1 for the unfiltered signal, the slope of the filtered signal returns to zero only for angles close to .

Finally, we inspected preprocessing by signal differencing, a common digital signal processing technique, to further test the robustness of the observed scaling artifact (see Figure 8). While signal differencing does remove much of the scaling artifact, a scaling component remains apparent in the resulting spectrum (see Figure 8 from 10 to 500 Hz). We further tested numerical detrending of signals. In Figure 9 we added frequency noise and in Figure 10 we added frequency and amplitude noise. This, along with the filtering tests presented above, demonstrates the difficulty in removing spectral artifacts resulting from signals outside the DFT cutoff and for combinations of signals more generally. Thus, while signal differencing is indeed useful for detrending experimental data in general, a signal preprocessed in these ways may still contain the spectral scaling artifacts observed here.
Figure 8:

Test with signal differencing. FFT spectral estimates for an example 10 Hz sinusoid (black), the same signal with added low-frequency component (blue), and this second signal preprocessed by differencing (red). The observed spectral scaling artifact remains in the differenced signal, illustrating the difficulty in removing this artifact through standard signal processing techniques.

Figure 8:

Test with signal differencing. FFT spectral estimates for an example 10 Hz sinusoid (black), the same signal with added low-frequency component (blue), and this second signal preprocessed by differencing (red). The observed spectral scaling artifact remains in the differenced signal, illustrating the difficulty in removing this artifact through standard signal processing techniques.

Figure 9:

Test with numerical detrending. (Left) A 10 Hz sinusoid with an added noisy oscillation (, where represents frequency noise). The original signal is plotted in the top panel, and the detrended version is plotted at bottom. (Right) The spectra for the 10 Hz sinusoid (black), 10 Hz with added noisy, low-frequency component (magenta), and the detrended version of the composite signal (blue).

Figure 9:

Test with numerical detrending. (Left) A 10 Hz sinusoid with an added noisy oscillation (, where represents frequency noise). The original signal is plotted in the top panel, and the detrended version is plotted at bottom. (Right) The spectra for the 10 Hz sinusoid (black), 10 Hz with added noisy, low-frequency component (magenta), and the detrended version of the composite signal (blue).

Figure 10:

Test with numerical detrending under the influence of added white noise (SNR = 10 dB). (Left) A 10 Hz sinusoid with an added noisy oscillation (, where and represent amplitude and frequency noise, respectively). The original signal is plotted in the top panel, and the detrended version is plotted at bottom. (Right) The spectra for the 10 Hz sinusoid (black), 10 Hz with added noisy, low-frequency component (magenta), and the detrended version of the composite signal (blue).

Figure 10:

Test with numerical detrending under the influence of added white noise (SNR = 10 dB). (Left) A 10 Hz sinusoid with an added noisy oscillation (, where and represent amplitude and frequency noise, respectively). The original signal is plotted in the top panel, and the detrended version is plotted at bottom. (Right) The spectra for the 10 Hz sinusoid (black), 10 Hz with added noisy, low-frequency component (magenta), and the detrended version of the composite signal (blue).

Acknowledgments

This work was supported by the Howard Hughes Medical Institute and the Crick-Jacobs Center for Theoretical and Computational Biology; the U.S. Office of Naval Research under grants N00014-10-1-0072 and N00014-12-1-0299; and NIH under grants R01-EB009282 and T32-EY20503-5. We thank Andrew Viterbi for helpful discussions.

References

Cooley
,
J.
,
Lewis
,
P.
, &
Welch
,
P.
(
1969
).
The finite Fourier transform
.
IEEE Transactions on Audio and Electroacoustics
,
17
(
2
),
77
85
.
Erdélyi
,
A.
(Ed.). (
1954
).
Tables of integral transforms
, vol.
1
.
New York
:
McGraw-Hill
.
Fourier
,
J.
(
1822
).
Thïéorie analytique de la chaleur
.
Chez Firmin Didot
.
Goertzel
,
G.
(
1958
).
An algorithm for the evaluation of finite trigonometric series
.
American Mathematical Monthly
,
65
(
1
),
34
35
.
Harris
,
F. J.
(
1978
).
On the use of windows for harmonic analysis with the discrete Fourier transform
. In
Proceedings of the IEEE
,
66
(
1
),
51
83
.
Küpfmüller
,
K.
(
1928
).
Über die Dynamik der selbsttätigen Verstärkungsregler
.
Elektrische Nachrichtentechnik
,
5
(
11
),
459
467
.
Lainscsek
,
C.
, &
Sejnowski
,
T.
(
2015
).
Delay differential analysis of time series
.
Neural Computation
,
27
(
3
),
594
614
.
Milstein
,
J.
,
Mormann
,
F.
,
Fried
,
I.
, &
Koch
,
C.
(
2009
).
Neuronal shot noise and Brownian 1/f behavior in the local field potential
.
PloS One
,
4
(
2
), e4338.
Mitra
,
P.
, &
Bokil
,
H.
(
2007
).
Observed brain dynamics
.
New York
:
Oxford University Press
.
Nuttall
,
A. H.
(
1981
).
Some windows with very good sidelobe behavior
.
IEEE Transactions on Acoustics, Speech and Signal Processing
,
29
(
1
),
84
91
.
Nyquist
,
H.
(
1928
).
Certain topics in telegraph transmission theory
.
Trans. AIEE
,
47
,
617
644
.
Titchmarsh
,
E.
(
1948
).
Introduction to the theory of Fourier integrals
(2nd ed.).
Oxford University
:
Clarendon Press
.
Ward
,
L.
, &
Greenwood
,
P.
(
2007
).
1/f noise
.
Scholarpedia
,
2
(
12
),
1537
.
Revision 90924
.