A guide on spectral methods applied to discrete data -- Part I: One-dimensional signals

Spectral analysis in conjunction with discrete data in one and more dimensions can become a challenging task, because the methods are sometimes difficult to understand. This paper intends to provide an overview about the usage of the Fourier transform, its related methods and focuses on the subtleties to which the users must pay attention. Typical questions, which are often addressed to the data, will be discussed. Such a problem can be the issue of frequency or band limitation of the signal. Or the source of artifacts might be of interest, when a Fourier transform is carried out. Another topic is the issue with fragmented data. Here, the Lomb-Scargle method will be explained with an illustrative example to deal with this special type of signal. Furthermore, a challenge encountered very often is the time-dependent spectral analysis, with which one can evaluate the point in time when a certain frequency appears in the signal. The information to solve such problems and to answer this questions is spread over many disciplines ranging from mathematics, electrical engineering, economic science to astrophysics. The goal of the first part of this paper is to collect the important information about the common methods to give the reader a guide on how to use these for application on one-dimensional data. The second part of this paper will then address the two- and more-dimensional data. The introduced methods are supported by the spectral package, which has been published for the statistical environment R prior this article.


Introduction
Because the field of interest is as wide as the amount of information about spectral methods, it becomes necessary to compose this article out of the most important methods. It is split into three main parts, which point at different levels of abstraction and complexity. The second chapter will give a brief and short overview over the mathematical principles and ideas behind the common methods of spectral analysis. For further reading the appendix contains a selection of derivations. Remember, the intention of this script is not the completeness of mathematical proofs. Its focus lies on the application of methods on discrete data and the mistakes to be avoided. The third chapter will focus on simple one-dimensional methods, their properties, advantages and drawbacks.

Mathematical concepts
This chapter introduces the basic ideas and principles which lead to the methods of spectral analysis for discrete data. Thereby, the focus lays on measurements and time series of physical processes. For now, utilizing the time as variable makes it easier to understand the explained methods. Of course, the time variable t might be replaced by the space variable x, if a distribution of spatial events is of interest. So in this sense space and time are meant to be the location where the signal is defined.
The following definitions and explanations are short and sweet. For further reading it is suggested to refer to "The Scientist and Engineer's Guide to Digital Signal Processing" by Smith (1997) or to "Time-Frequency Analysis" by Cohen (1995).

Definition of a signal
In the context of this article a signal s(t) corresponds to a physical measure p(t) and is therefore real-valued and causal. This means, that with the measurement of the process p(t) the signal s(t) starts to exist at a certain point in time and ends later when the measurement is finished. With that in mind, the signal function s(t) represents a slice of the length T for times t > 0. This can be properly defined as follows: The simple noisy signal in Figure 1 illustrates that. For this example, the underlying physical process p(t) could be the temperature, which is measured in the time range between 0 < t ≤ 1. Evidently the temperature of an object exists before and after the measurement takes place, so the content of s(t) only maps to a time interval of this process.
With respect to a digital representation of information, the measurement of p(t) takes place by acquiring data in the form of sampling the continuous signal at N different instances in time. The resulting sampled signal serieŝ is the data vectorŝ n . Here, the index n addresses the n th element in this N -dimensional data vectorŝ. In equation (2) the basis vector e n ensures thatŝ becomes a vector in a mathematical sense, which simply represents an ordered list of elements. Remember, the absolute value | e n | = 1 is always equal to one, so the sum only helps to iterate along the continuous signal s(t). The arbitrary signal p(t) = 0.4 cos(2π ·1.5t−π/2)+N (0, 0.02)+0.5 with overlapped normal distributed noise is shown here. Only within the rectangular window the measurement takes place and the signal s(t) can exist. Outside this window the physical process might go on, but is not recognized anymore. The points and arrows mark the sampling of that signal.

The fundamental idea of frequency analysis
To explain the idea of frequency analysis right from the beginning, it is necessary to redefine the signal function s(t). From now on, the signal represents a sum of sinusoidal functions in the form: which, in the final consequence, leads to the expansion into orthogonal functions.
The aim of frequency analysis is to calculate the amplitude and phase spectrum of a signal. The result of this procedure is a mapping of amplitudes to their corresponding frequencies. If a measured time series of a physical process is the sum of many individual sinusoidal functions, it could be of interest to estimate each of the amplitudes A i and phases ϕ i for a given frequency f i .
The first approach to do that would be to estimate each A and ϕ by a least square fit against a sinusoidal model for each frequency of interest. This is a legitimate attempt with certain drawbacks and some outstanding advantages, which are discussed below in section 3.8.
Given a signal in the continuous domain s(t) = A · cos(ωt + ϕ), then multiplying this with cos(ω 0 t) as well as sin(ω 0 t) and finally integrating over all times yields for the cos-term: = R (real part or inphase component).
which is going to be named R. Herein, the variable ω 0 denotes the frequency of interest. The sin-term is very similar = I (imaginary part or quadrature component) and is called I. In both equations on the right hand side the right integral vanishes to zero, because each positive area element of a sin-function can be mapped to a negative area element. Only in the case of ω 0 = ω a phase-dependent part is left,otherwise the limit of the integral converges to zero.
With the identity sin 2 α + cos 2 α = 1 the amplitude can be calculated as The phase information is then gained by ϕ = arctan I R The introduced procedure is called quadrature demodulation technique (QDT). Finally, by evaluating the amplitude A at many different frequencies ω 0 , a periodogram emerges. Note, the R and I part can also be interpreted as the real and imaginary part of a signal. The complex exponential formulation then corresponds to with the identity e iφ = cos φ + i sin φ.
One of the main disadvantages arise for finite time ranges 0 < t ≤ T . In that case, this method produces a truncation error ε, which depends on the remainder after the modulus division n = (2πf · T ) mod (2π) with respect to the total time interval T . In other words, if the period 2π/ω of the signal does not fit completely in the integration range, the integrals containing ω + ω 0 will not vanish completely. However, remember the right term of equation (7). Here the integral can be split into an n · 2π periodic part and a residue. The periodic integral of cos ((ω + ω 0 )t + ϕ) vanishes, whereas a residue ε remains in the second integral from the last n · 2π period up to the full length T . Obviously, a minimum of ε is reached if T = n·2π. All this holds if ω = ω 0 . If this is not the case, also the integral cos ((ω − ω 0 )t + ϕ) splits into two parts from 0 to n · 2π and a residue, which leads to a second error ε 2 . Finally, with ω ≈ ω 0 this leakage effect will produce artificial amplitudes in the nearest neighborhood of ω 0 . Any way, note that from inserting the result (15) in equation (11) or (12), mixed products arise, which makes it difficult or even impossible to reject the error ε. All this is far away from being mathematically complete, but it points in the direction towards an error discussion, which is presented in more detail by Jerri (1977, Chap. 6).
At this point, the statistical approach "fitting a sinusoidal model to the data set of interest" becomes attractive. A brief discussion of this concept is given in "Studies in astronomical time series analysis. II" by Scargle (1982). Chapter 3.8 provides an introduction on the Lomb-Scargle method, which is equal to the already mentioned least square fit.

Integral transforms
In contrast to the previously described QDT, which is more or less a statistical approach, the integral transforms are the mathematical backbone of spectral analysis. An integral transform represents the bijective mapping s(t) ↔ S(ω) between the time and frequency domains, which states that all information present in the first will be transformed to the second domain. In the following, capital letters denote always the frequency domain, whereas lowercase letters represent the time domain of the corresponding variable.
Fourier transform First of all the Fourier transform should be introduced. Thereby, the definitions and some of the properties in the continuous and discrete time domain are different. A comparison is given in table 1. Nevertheless, the signal s(t) forms a Fourier pair with S(ω).
The following notation describes a transform of the signal s(t) into the frequency domain, This uses the Fourier operator F, which expresses the calculus according to table 1. The inverse back transform is denoted by F −1 . The Fourier transform gains different features like the frequency shifting which becomes important when amplitude modulation in chapter 3.4 is discussed. All necessary properties are listed in Appendix B.
In addition to the characteristics of the continuous Fourier transform, the discrete Fourier transform (DFT) has a 2π-periodic spectrum. As shown in Appendix C, the result of the DFT is a complex valued data vector, in which the first elementŜ 1 always contains the mean value of the data series. Subsequent bins hold the amplitude values up toŜ N/2 . From this position on the spectrum is mirrored and repeats up to the sampling frequency. In chapter 3 the symmetry of the DFT's and its consequences are discussed in more detail. It will be shown how to take advantage of this to recover the "negative" frequencies. In case of more dimensional data, it turns out that the result of the DFT contains additional information, where also the "negative" frequency components play a role. Further insights on this will be given in part II of this paper.
Additionally, the linear operator H is introduced to execute the Hilbert transform on a function or data set. From the representation in the Fourier space the constant phase shifting feature becomes clear, Hence, equation (19) can be used to easily calculate H (s(t)) in terms of a Fourier transform. Remember, here the sign(ω) function is defined on −∞ < ω < ∞. This must be taken into account, when performing H(ŝ(n)) on discrete data. In this case, the discrete Hilbert transform is calculated by This formalism considers the fact that "negative" frequencies are mirrored into the upper parts of the DFT data vector. However, the result of H (s(t)) is real-valued, with each frequency component phase-shifted by ±π. In conjunction with a real-valued and causal signal s(t), the HT helps to represent the real physical character of measured data, which should only contain positive frequencies, due to real worlds processes.From this, the analytical signal follows, which is a complex representation of s(t) containing a one-sided Fourier spectrum.
One derivation of the Hilbert transform in conjunction with the analytical signal is given in appendix D and its application in terms of data filtering is briefly discussed in section 3.3 and the following.
In addition to that, the concept of instantaneous frequency can be explained with the help of the HT. This topic it discussed in detail in "The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis" by Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung, and Liu (1998). For further reading on that "Amplitude, phase, frequency-fundamental concepts of oscillation theory" by Vakman and Vaȋnshteȋn (1977) should be mentioned.

Spectral methods in one dimension
The application of the introduced Fourier and Hilbert transform will be discussed in detail in this chapter. Each of the following sections includes a theoretical part, which is finally supported by examples. The focus lies on the application of the corresponding methods.
First of all, let the signal function be a set of sinusoidal functions plus an optional noise term. This signal is sampled with an equidistant spacing T s , so the discrete data vector s(n · T s ) · e n of length N represents the sampling series. Using the discrete Fourier transform, see table 1, the spectrum of the data vectorŝŜ can be calculated.
For further reading a brief introduction can be found in "The Shannon sampling theorem-Its various extensions and applications: A tutorial review" by Jerri (1977). A detailed comprehensive mathematical work is given in the book "Fourier Analysis" by Duoandikoetxea Zuazo (2001).

The essence of band limitation and the Nyquist condition
The question of band limitation is a question of the unique representation of a function s(t) in the frequency space S(ω). In the simplest case, a band limited signal has a spectrum which is only non-zero within the range |ω| ≤ ω max . Outside this interval the spectrum is equal to zero. The sampling theorem by Shannon (1949) and the proof of Whittaker (1935) states, that in such a situation even a time-limited continuous function is completely defined by a discrete and finite set of samples, if all the samples outside the time range are exactly zero (see the exemplary function in Figure 1). Closely related to that is the minimal required sampling frequency (twice the "Nyquist frequency") 1/T s = f s ≥ 2f max , which is necessary to obtain all information about the signal.
The consequence is that there should not be any frequency present in the data set higher than f s /2 to achieve a unique sampling in application. Theoretically, this means a function s(t) is completely determined by its discrete spectrum S(n · 2πf s ) if the signal is bandlimited to |f | < f s /2. Therefore the discrete spectrum of the discrete sampled signalŝ n provides a unique spectral representation of its continuous counterpart s(t).
Application The artificial function y(t) = sin(2π · 4t) + 0.5 · cos(2π · 2t) + 1.5 is given in Figure 2a. It consists of two functions and an offset value. The sampling took place at a rate of T s = 0.05, so that there are theoretically N = 21 possible samples. Note, y(t) is perfectly periodic, so the first sample would equal the last sample. Because the latter marks the first point of the next sampling period, it must be removed from the data set, which now has N * = 20 elements. Ignoring this fact would lead to an error (leakage effect) with the result that the elements to the left and to the right of the frequency peaks are not Mention, that the first sample and the i th + 1 sample at t = 1 would be the same, because of the periodicity of the signal.
zero anymore. The root of this behavior is discussed in section 2.2 and becomes visible in the subsequent figures in this section. The resulting frequency resolution in this example is = 1 so the signal's frequencies fit perfectly into this grid. The absolute amplitude spectrum Ŷ n is given in Figure 2b. Here, one can see that the mean, saved in the first bin, has its original value of A 0 = 1.5, whereas the magnitude of the amplitudes at f n=2 = 2 and f n=4 = 4 are split into the upper and lower frequency domain. To retrieve the correct signal amplitudes the individual values in both half-planes must be added. Remember, the DFT is mirror-symmetric to the Nyquist frequency f s /2 = 10 so that the first amplitude A 1 becomes Figure 2 and Figure 3a display an over-sampled signal, where the required Nyquist condition is fulfilled. The result is a unique mapping of the time domain into frequency domain. In other words, the corresponding frequency vectorf n is only valid for 0 ≤ n < N/2, because it repeats in the inverse order for N/2 ≤ n < N .
However, as Figure 3c indicates, if the requirement of band limitation (f s ≥ 2f max ) is violated, the mapping is not unique anymore. Remember, the frequency resolution in all given examples is δf = 1. In comparison to Figure 3b, where the condition of band limitation is complied with, Figure 3d illustrates that the upper and lower parts of the frequency bands now infiltrate each other. This makes it even more complicated to distinguish between certain frequencies and their physical correspondence. In such a scenario it becomes completely unclear which of the two functions belong to the sampling points.
Example -Under sampling The conditions described above can be extended to the requirement which means that a unique mapping is still possible, yet only if the bandwidth ∆f is not larger than the sampling frequency. In addition to this the frequency range must be known. This is automatically fulfilled in the most common case, whenever the frequency ranges from 0 ≤ f ≤ f s .
The example given in Figure 4 shows the artificial function y(t) = cos(2π · 25 · t) + 0.5 · sin(2π · 27 · t) + 1.5, which is sampled with a sampling frequency of f s = 20. The amplitude spectrum is given in Figure 4b. Obviously y(t) is under-sampled. Remember, the discrete Fourier transform calculates a periodic spectrum, which projects the high-frequency parts into the lower frequency range between 0 and f s . Without any additional information a mapping from spectral domain to time domain will fail. In the present case we assume that the signal frequency is in the range 25 ≤ f ≤ 27, so the whole spectrum can be shifted by an offset of f o = f s . The second f -axis indicates this. The upper and lower frequency bands can be clearly distinguished, so the mapping is unique. In contrast to that, the interprenetation of the upper and lower bands prevents such a mapping in the example in Figure 3d.
Example -The derivative of a function One important property of the Fourier transform is the change of mathematical operations in the spectral domain. According to Appendix B, taking the derivative of a function reduces to a multiplication with iω in the spectral domain. In this manner, differential equations can be reduced to simpler linear algebraic representations.
In contrast to the continuous case, discrete sampled functions must be periodic with a finite bandwidth to fulfill the requirements for the DFT. In case the function is non-periodic within In (b) the analytic (gray line) and estimated numeric (points) derivative dy(t)/ dx is shown.
The crosses indicate the result of the calculation without utilizing the window function. Note, that some crosses are vertically outside the window for x < −2 and x > 2.
the sampled window, the application of a proper window function is necessary to minimize leakage. The following code explains how to estimate the derivative of the non-band-limited function like which is going to be sampled with N = 40 discrete points. require(spectral) x <-seq(-2.5,2.5,length.out = 40) # 40 sample positions y <-win.tukey(x,0.2) * (-x^3+3*x) # window the function # and generate the # sample points Y <-spec.fft(y,x,center = T) # doing the DFT Y$A <-1i * 2*pi * Y$fx * Y$A # calc. deriv. # the realpart of the back transform # contains the signal of interesst dy <-base::Re(spec.fft(Y,inverse = T)$y) It must be mentioned that the function y(x) must be windowed in advance to achieve the periodic property. In this example this is done via the "Tukey"-window, which rises and falls with a quarter cycle of a cosine function (compare dashed curve). Figure 5 shows the output of the listed code above. Pay special attention to the windowed version of y(x) (black curve) and its almost perfect derivative. Here, the effect of the windowing process is to achieve Y (|ω| > 0.5 · ω s ) = 0, which relaxes the function at the right and left side. In consequence, the resulting derivative equals the analytic solution in the mid range, but also deviates at the boundaries. In comparison to that, the unweighted function in (a) (gray line) produces an unusable faulty derivative in (b) (crosses). The reason for that is the discussed missing band limitation.
With a look into line 7 of the above code, notice that the additional parameter center = T is provided to the spec.fft() function. The effect is a zero-symmetric spectrum with the corresponding frequencies ranging from −f s /2 ≤ f < f s , which enables the direct computation of iω · Y (ω) utilizing the elements of the list variable Y. Remember, a pure DFT would map to frequencies in the range of 0 ≤ f < f s . It follows that the negative frequencies −f s /2 ≤ f < 0 were projected into the upper half of the frequency domain f s /2 ≤ f < f s . In case of center = F, the latter frequency range would be generated and used in further calculations, which must be considered to prevent wrong results.

Center a discrete spectrum
As initially explained in section 2.3 the spectrum of a continuous function s(t) is defined in the range of −∞ < ω < ∞. In contrast to that, its counterpart the discrete Fourier transform produces a 2π-periodic spectrum, which is defined for a positive number of samples 0 ≤ n < N and 0 ≤ m frequencies. The negative frequency components are projected into the the range of f s /2 ≤ f < f s , so the spectrum becomes mirror-symmetric with respect to the Nyquist frequency f s /2. Figure 6a shows the result of the continuous Fourier transform applied to the artificial example function (24) y(t) = sin(2π · 4t) + 0.5 · cos(2π · 2t) + 1.5.
It is clear that the mean value y(t) is located at f = 0 and the sinusoidal components are located symmetricly around zero. In contrast to that, the discrete Fourier transform produces a data vectorŶ m , to which only positive discrete frequencies can be assigned. To convert the data vector into a zero-symmetric representation, the function y(t) must be modulated by This method utilizes the shifting or damping property of the Fourier transform as stated in Appendix B Equation 4, assuming a damping factor of a = 1. Doing so the complete spectrum ( Fig. 6b), including the mean value shifts towards the half of the sampling frequency right in the middle, as it is indicated in Figure 6a.
Application In an equally spaced sampling series the discrete time vector t = n·T s changes equation (31) to and simplifies it. From equation (32) follows that the modulation of the data vector can be computed easily by a multiplication of the alternating series (−1) n . Subsequently, all this can be extended into n dimensions, so even 2D spectra can be centered (Gonzalez and Woods 2002, p. 154).  The following R-code which is part of the spec.fft() function explains the method. First, assume a vector x for the spatial location where the measurements y were taken.
The back transform of the centered spectrum is the simple reverse of the described procedure. The function spec.fft() supports both methods to decompose a data vector into its centered or non-centered spectral representation. The usage of the code is explained below.

The analytic signal
The concept of the analytic signal is based on the idea that a real-valued signal, e. g. some measured data can only contain positive frequencies. At least in the 1D case, the interpretation of negative frequencies is almost impossible. Next to that, a real-world signal does only exist for positive times, usually. This makes the signal causal, which also must be taken into account, see Equation (2) in Appendix D. The foundation of the analytic signal is the Hilbert transform. For a short introduction on the Hilbert transformation itself and its properties, refer to section 2.3. Besides that, the attribute "analytic" originates from the Cauchy-Riemann conditions for differentiability, which is fulfilled for analytic functions in this sense (Cohen 1995).

Now, lets introduce the analytic signal
It is formally defined by the sum of a signal s(t) and its Hilbert transform iH (s(t)) multiplied with the imaginary unit. So the real-valued signal s(t) transforms to a(t), which now has an one-sided spectrum. Frequencies above the Nyquist-frequency f s /2 are coerced to be zero and frequencies below that are gained by the factor of two, so that energy is conserved. This can be seen in Figure 7, which shows the result of a DFT applied on the analytic representation of the sampled signal from Figure 2. Starting from that several applications can be derived.
First, all the calculations which utilize the signal's frequencies, e. g. the estimation of the derivative introduced in section 3.1, get simplified with respect to discrete data because the distinction of frequencies above and below f s /2 can be omitted. The example of the derivative in section 3.1 instead utilized the discrete centered spectrum, which provides a frequency vector with positive and negative parts, which can be multiplied directly to the spectrum.
Second, the imaginary part of a(t) equals the Hilbert transform of its real part. The Hilbert transform's phase shifting properties form the basis to calculate, for instance, the envelope function, which is discussed in the next section. The analytic signal also provides a way to estimate the instantaneous frequency of a signal. The latter ends up in the empirical mode decomposition, which gives a time depending spectral decomposition according to Huang et al. (1998), which is out of scope of this guide. In section 3.7 it will be shown how to utilize the Fourier transform and how to overcome its draw backs in conjunction with the analytic signal of data of a non-stationary processes.
Implementation According to Appendix D the function analyticFunction(y) provided by the spectral package calculates the analytic signal. An example is given in Figure 7, which displays the analytic representation of the sampled artificial example function (24) from section 3.1. Here the spectrum is single sided and all the amplitudes have their correct value. Note that all components above f s /2 are zero, because the upper half of the spectrum is projected into the lower half by this method. In consequence, the illustration in Figure  The program code # correct mean value ht <-fft(X, inverse = T) explains how the the analytic signal is calculated. Here the part (1 -sign( f -0.5 )) solves three things. First, this statement shifts the phase with respect to the sign of the virtual frequency vector f, whereby the "negative" frequencies are located in the upper half of the data set. This provides a better solution than Equation (20) because the evenness of the data sets length does not matter anymore. Second, since f is an integer vector, the subtraction f -0.5 circumvents the problem sign(0) = 0, which avoids an error with even length data sets at f = 0. And third, the above code calculates the Hilbert transform and the analytic signal in one step in the spectral domain by combining equations (19) and (31).
The code below produces the single-sided spectrum Y given in Figure 7. Note, that the amplitudes correspond to the real input values of the function (24) and all spectral parts above f s /2 are zero.

Calculating the envelope
Sometimes it becomes necessary to calculate the envelope function of data. Different ap-proaches are possible, for instance finding all maxima and then fiting a spline function to these points. However, in the following the spectral way to do that is introduced.
First, remember the trigonometric identity which becomes the key component of the calculation later. Next to that, assume a function with A(t) as the envelope function, which is modulated with the carrier cos(ω 0 t). It is clear now that calculating the envelope is equivalent to an amplitude demodulation.

Now, the amplitude function is gained by
where y * (t) denotes the signal y(t) phase shifted by π/2, so that y * (t) equals the Hilbert transform Under the condition of sufficient band limitation the statement above can be expressed as which works, because (37) acts like an ideal phase shifter as discussed in section 2.3 and Appendix D.
Application The calculation is straight forward and follows Equation (38). The spectral package includes the function envelope() to perform the calculation of an 1D envelope. The only problem with this is the band limitation, which is necessary to achieve reasonable results. Compare Figure 8b, here the required bandwidth to demodulate the envelope becomes clearly visible in the negative and positive half plane.

Convolution
The Fourier transform does not only map a spatial function into the spectral domain, moreover it also converts several mathematical operations. One example is the derivation of which equals a simple multiplication with the complex variable iω in the spectral domain.
In the following, the convolution ∞ −∞ f (τ )g(t − τ ) dτ will be introduced by the identities and which are valid only if the integrals Example -Polynomial multiplication Suppose two polynomials of the degree n = 2. Then the multiplication of these two would be degree expansion FFT   8  64  209  16  256  500  32  1024  1175  64  4096  2714   Table 2: Computation costs for the polynomial multiplication. The calculation amount for the FFT is estimated by 3 · (2n + 1) log 2 (2n + 1).
which finally gives an expression of degree of 2n = 4. Performing the expansion will end in a convolution of the two coefficient vectors, which must be 2n + 1 elements long to fit the result in the output vector. Doing all that will spend (n + 1) 2 = 9 floating point multiplications on a computer. The statement expresses the convolution and the result, where the symbol * denotes the folding operator. Alternatively, this operation can be done with the help of the Fourier transform. According to (40) requires the transform to be calculated three times.
The great advantage of using the Fourier transform, instead of the straight forward expansion, is the scaling of n · log 2 n for the implementation of the FFT (Cormen, Stein, Leiserson, and Rivest 2001, Chap. 30). As shown in table 2 the FFT accelerates the calculation from a value of n > 32 for this example.
The standard algorithm of the discrete FFT is limited by a vector size of N = 2 n , which can be overcome with the implementation of the DFT method as described by Frigo and Johnson (2005), which even allows prime numbered lengths without a loss of performance.
As a simple example the following code illustrates the method described above. The result is given in the last line of the program listing. The R-command convolve(x, y, conj = F) uses the same mechanism and produces exactly the same results.

Spectral filtering in conjunction with the sampling theorem
Every sampled signal has its starting point t = 0, followed by a point in time T where it ends. The sampling procedure itself and the nature of the signal finally define the underlying band limitation. The intention of the first part of this section is to show the relation between the measurement and the consequences which arises if a causal signal is sampled. It turns out that in conjunction with the Fourier transform several convolutions take place in the time and frequency domain by selecting the maximum time T and the maximum frequency.
Regarding that, it is important to understand that the filter process already begins with the sampling of the signal. Figure 9 illustrates that for the function which fits perfectly into the period of the sampling window. In principle, this function is defined on the interval of −∞ < t < ∞. Then, the data acquisition process defines a starting point at t = 0 where the measurement begins and an end point t = T where the measurement stops. This is illustrated by the bold rectangular window in Figure 9a. This window function picks out a range 0 ≤ t < T of y(t) so that the signal is now defined by y(t) multiplied by w(t). According to the convolution properties of the Fourier transform (see section 3.5), this operation corresponds to the convolution of the ideal spectrum Y (ω) = F(y(t)) with the spectral representation of the window function w(t), which is in fact Shannon's sampling function (Shannon 1949). The consequence for the signal's spectrum S(ω) is that the ideal Dirac-spectrum Y (ω) smears out to sin(x)/x like functions, because each arrow is weighted with W (ω). Note, the zeros of the resulting spectrum are in an equidistant spacing of δf = 1/T now. At this moment the original function y(t) is just windowed and not sampled.  The next step is to limit the spectrum in Figure 9b to a maximum frequency, e. g. f max = 1/2T s , by applying an additional rectangular window in the frequency domain The parameter T s becomes the period of the sampling frequency. The important link to the Nyquist frequency would be clarified later. However, the corresponding representation of is displayed by the dashed line in Figure 9a. In the context of the convolution, which now takes place in time domain, the windowed signal is finally convoluted with v(t). But still, the result (52) is not sampled yet which is done in the next step.
Sampling is a procedure where individual values of the function s(t) are selected and stored. With regard to the nature of the sinc-like functions in Figure 9, the sampling should take place at each possible zero of v(t), as well as v(t = 0) which is the only sampling point with a non-zero value of v(t). This behavior guarantees that the convolution (52) rejects everything, but the points of sampling. In the same moment also the distance in between them is set to δt = T s , because the function v(t) provides two zeros in an interval of (2T s ) −1 · t ≤ 2π, which correspond to the next non-zero samples of the neighboring points. In other words, by  choosing the next sample v(t) must be shifted by exactly T s to reject the former sample and to take the current value. Note, violating this restriction, because the location of sampling points jitters, will lead to an error in the discrete spectral representationŜ m .
Finally, this conceptual description gives a second access to the question why the sampling frequency should be at least twice the maximum signal frequency. A detailed mathematical derivation of the discrete Fourier transform and the properties of sampling is given in the book of Debnath and Shah (2015, Chap 3). However, the discrete sampling given in Figure 9a and the finite length of the sampling series causes the convolution (48) to take the only possible peaks in the spectrum at f = {−2, 2}, because all other 2π-periodic values evaluate exactly to zero.
If the signal period does not fit perfectly into the window, which is the case in reality, the kernel function W (ω) never matches the true Dirac, but instead of that it sums up values into the left and right bin of the corresponding discrete Fourier vectorŜ m . To reduce these side effects, other window functions can be used in advance, before the DFT takes place. The intent of the window functions (Hamming-, Blackman-or Tukey-window) is to convolve the spectrum with a steeper decreasing function amplitude compared to W (ω) to minimize the effect of the period's misfit. Moreover w(t) has an infinite spectrum, which is then truncated in the finite and discrete data processing, which leads to additional side effects if the period does not fit the sampling window. However, window functions can help to reduce such interferences, which are generated when local non-periodic events or infinite spectra are present.
Application Convolution filter in the time domain, such as the moving average filter, are easy to calculate with the help of the FFT. Given the signal s(t) = cos(4π · t) + sin(20π · t) + N (0, 0.5)| 0≤t<1 , which is defined in the range of 0 ≤ t < 1. The moving average of the sampled signalŝ n with the filter kernel length of N K = 5ŝ is then calculated for the n th element by taking the weighted sum of N K elements before that element. The filter kernelk for the example in Figure 10 is represented by the weighting coefficientsk With respect to the spectrum ofk, it becomes clear that the length ofk is chosen in a way that the high frequency part sin(20π · t) is rejected completely. But it is also evident that other frequency components still remain in the result or were damped falsely, like the low frequency part cos(4π · t). Figure 10a illustrates that the application of the moving average filter in the time domain will lead to a phase-shift in the signal. Working in the spectral domain, one will gain the possibility to shift the signal back and reconstruct the whole definition range as follows The result of this operation illustrates the dash-dotted line in Figure 10a. The code which calculates this is presented below.
With respect to the spectrum of the kernelK it becomes evident that the moving average filter is not able to suppress higher frequencies completely.  Example -The ideal low pass filter Since the moving average has some drawbacks (because of its infinite frequency response) one should think about an ideal (low pass) filter with a finite frequency response. Such a filter has the advantage that unwanted spectral components can be rejected completely. On the other hand, this results in an infinite long kernel in the time domain, which is nevertheless calculated in the spectral domain.
The basic principle of such a filter is to set each component above a certain frequency |ω max | value to zero S(|ω| > ω max ) = 0. (57) Figure 11 illustrates the procedure. Compared to the previous moving average example in Figure 10, the output signal's shape changed to a more smooth and accurate form. However, there remains a difference compared to the ideal low frequency part (dashed line in Figure 11a), which is the result of the remnant spectral components around the signal frequency f = 2.
A code snippet to calculate the ideal low pass is presented below.
# assuming x and y to hold the time vector # and the sampling points Y <-spec.fft(y, x) # doing the filtering Y$A[abs(Y$fx) > 3] <-0 # reconstruct the filtered signal yf <-spec.fft(Y, inverse = T)$y The spectral package provides the filter.fft(y,x,BW,fc=0,n=3) function to do that in one step. It utilizes the analyticFunction(y) to filter the one-sided spectrum and enables the user to apply a band pass filter with the arbitrary polynomial weight function The parameter n sets the polynomial degree, whereas 3 is close to the moving average. Passing higher values, say n=10, will successively converge to the ideal bandpass solution. All this is because some times it is desirable to have a smooth weight of frequencies and some times not.
The used code then changes to the following example.
yf <-filter.fft(y, x, fc = 0, BP = 3, n = 10) As seen in Figure 11 the bandwidth BP is symmetric around the center frequency fc = 0, so the values must be set in accordance to the desired window which should remain in the spectrum. Utilizing the analytic signal representation avoids mistakes when programming the weight function and keeps the code clean.
Example -Noise reduction via autocorrelation By choosing an appropriate kernel function it becomes possible to extract certain features of the signal. The previous example illustrated how a low pass filter can filter out a noisy signal by zeroing all upper frequency components in the spectrum. In the time domain this corresponds to a folding operation of the input data with an infinite "sinc"-like kernel function.
Now the opposite approach is to ask what the most significant periodic components within the signal are. First of all, the continuous autocorrelation function provides a mechanism to examine a data set or function with respect to its self-similarity. The equations above convolute the signal s(t) with its complex conjugate s * (t). According to the convolution properties, equation (59) can be expressed as the inverse Fourier transform of the product of the signals spectra. Note, for a real-valued signal s(t) = s * (t). Suppose the signal consists of a stationary sinusoidal signal it is quite evident that the acf rejects the noise which overlays s(t). The underlying noise equals to a non-stationary process and therefore it cancels out by utilizing the acf . A generalization of this statement is the Wiener-Khinchin theorem. Its application is discussed in detail by Cohen (1995, Chap 1.9), Marks (2009, Chap. 4) and Wiener (1949).
Assuming the function with some normal distributed noise overlaid. The signalŝ n is then the sampled function s(n · T s ) with T s = 0.01, so it fits perfectly into one sampling period of T = 1. This example is very artificial, real world signals often cause aliasing effects, so attention must be paid.
The discrete spectrum of the acf n can now be calculated by Note, the resulting spectrum is real-valued, so the phase information will be lost. A proper noise reduction can be achieved by defining a weighting vector which sets every spectral componentsmaller than the standard deviation of the ACF m to zero.
The resulting filtered signalŝ is given in Figure 12. For the given example it is important to calculate everything with the analytic function representation to preserve the energy content of the data vectors.

Non-stationary processes -spatially dependent spectral analysis
Many signals in reality do not fulfill the requirement of stationarity. It is very often the case that the measured signal is overlaid by a trend or an "temporal-local" event that takes place only once. All these in-stationary and non-periodic processes will spread out into all  frequencies in the corresponding Fourier decomposition of the data. If a proper band limitation cannot be achieved the mapping between physical frequencies and the spectral representation might fail too, as stated in section 3.1.
The equation · sin (2π · 40 · t) (65) describes a very artificial example of a non-stationary signal, which is illustrated in Figure 13a.
Here two different events occur. The first one has a low frequency at f 1 = 20 Hz, whereas the second one oscillates faster with f 2 = 40 Hz. Both parts are modulated with an Gaussian envelope to switch each of them on and off. Finally, the right panel (b) shows the resulting time variant decomposition of the signal y(t). The waterfall() function from the spectral package can be used to calculate such kinds of waterfall diagrams.
Note, this type of analysis produces a two-dimensional time-frequency map like a shifting FFT, which selects only a window of the data before doing the spectral decomposition for one time step. The amplitude functions A 1 (t) and A 2 (t) play the role of arbitrary window functions which mask the time series. The interpretation is that A i (t) contains the information in time where the finite process is located. All this relies on the assumption that the signal of interest can be modeled as In contrast to the shifting FFT the introduced approach is to do an amplitude demodulation by calculating the signal's envelope. The result is the amplitude function A(t) for a given center frequency f c , which corresponds to the frequency of the physical process. Finally, a sufficient bandwidth around f c guarantees that A(t) will be reconstructed correctly.
The first advantage of this method is that for each single frequency of interest the whole data set is taken into account. Note, the shifting FFT would select only a certain range of the data, so information about low frequencies can be lost if the selected window becomes too small. However, the overall frequency resolution decreases to the length of the window size of the shifting FFT. In the end, two signals that are very close in the frequency space might not be distinguishable anymore. The second point is the necessary band pass filter at f c , which is changed according to the actual frequency. This draws attention to the uncertainty principle, which states that with increasing frequency the locating of an event becomes sharper but the exact frequency gets more incorrect.
Implementation As stated above, the key component is a band pass filter BP, which must be applied twice to calculate the two envelopes in Figure 13. The following code shows how this can be done in R with the appropriate functions provided by the spectral package.
A1 <-Re( envelope(Re(filter.fft(y,x,fc = 20,BW = 10,n = 10))) ) A2 <-Re( envelope(Re(filter.fft(y,x,fc = 40,BW = 10,n = 10))) ) The next step is to do this calculation for many frequencies. A reasonable range is to start from f = 0 up to f s /2. To accelerate the code the waterfall() function uses a fast version of the envelope() function, which combines the filtering and the Hilbert transform as follows.
Finally, the bandwidth calculation is done by an empirical approach, like Here δf = 1/∆x denotes the frequency step and wd is the normalized width of the resulting band pass. The task of BW (f c ) is to widen the frequency band for higher frequencies, whereas low frequencies take a very small band width. This somehow pays attention to the uncertainty principle, which states that "a narrow waveform yields a wide spectrum and a wide waveform yields a narrow spectrum and both the time waveform and frequency spectrum cannot be made arbitrarily small simultaneously" (Skolnik 1980).
is going to be analyzed. Figure 14 illustrates the results. The waterfall diagram in panel (b) displays all the features of Equation (68). Note that even the |x| function and the sin x 2 term can be distinguished, whereas the normal time-invariant Fourier spectrum in Figure 14c hides this property completely. Here the 10 Hz carrier is the only "correctly" visible signal component.
A second example is given in Figure 15a. Here the given data represents the outflow temperature of a cooling system with about 300 kW power. It consists of a heat exchanger outside the building and a secondary water loop inside from which to draw the heat. A sprinkler system improves the performance of the outside heat exchanger, so that outflow temperatures below the ambient temperature become possible. Panel (b) illustrates the ambient temperature near the external air heat exchanger. The impact of the suns radiation becomes visible in a signal amplitude with a very long period of about T sun ≈ 80 min. In terms of period lengths the  high eigenfrequencies of the hydraulic system become visible at the bottom, approximately within a range of 10 min to 20 min, when the control loop reaches its limit.
Finally, the following code shows how to invert the frequencies to periods. Note, that a resampling of the matrix must be performed, because the unequally spaced vector 1/f now maps to the corresponding spectrum. This task is done by the plot.fft() function, which also plots objects with the attribute mode = "waterfall".

Fragmented and irregularly sampled data
Sometimes the nature of measurement prevents an equally spaced sampling. Observations of astrophysical bodies like the sun or the moon are examples of objects only visible at a certain time of day. Concerning a definite property of these bodies, the corresponding time series would be fragmented, what makes a spectral analysis with a standard Fourier transform almost impossible. A nice application example can be found in "The Lick Planet Search: Detectability and Mass Thresholds" by Cumming, Marcy, and Butler (1999). Another issue would be the measurement with randomly sampled data, for instance the occurrence of an event like the appearance of a malfunctioning product in a production lane, or even a measurement with a high jitter.
The time instances t n when a sample is taken now depend on the sample number n. This means, the time interval ∆t = t n−1 − t n = const.
is not constant anymore, so the sampling time becomes an additional data vector. This enables one advantage in case of randomly sampled data. Given a stationary process, it becomes possible to estimate signal amplitudes for frequencies which are in the order of f max ≈ O min(∆t) −1 , with the minimal distance between two points. In fact, this is related to an average Nyquist frequency f c , where which can be defined by the time range of the time series and the total number of samples.
Lets introduce the Lomb-Scargle periodogram as a consistent estimator for spectral properties of discrete data. The estimation in that sense is not a mathematical transform like the Fourier transform. The frequencies of interest can be selected freely, so the whole method analyses a data set against amplitude and phase. The associated mathematical details of this statistical approach can be found in the papers by Mathias, Grond, Guardans, Seese, Canela, Diebner, and Baiocchi (2004), Hocke andKämpfer (2009), Scargle (1982) and Lomb (1976), who originally developed the method.
In principle, the concept of the procedure is about a least square optimization of the parameters A i and B i for the model function which is fitted to the data. Remember, the introductory chapter 2.2 shows a slightly different approach in conjunction with the trigonometric identities (4) and (5), which also become the key component here. While the traditional calculation of a Newton algorithm takes several iteration steps until the result converges, the optimal parameters can alternatively be estimated by one single matrix inversion, see Mathias et al. (2004). This can be done only if the orthogonal property of the trigonometric model function is taken into account. The mathematical proof can be found in "A generalized inverse for matrices" by Penrose and Todd (1955). As a necessary condition to obtain an optimal result for A i and B i , the relation N n=1 cos (ω i · (t n − τ i )) · sin (ω i · (t n − τ i )) = 0 (71) must be fulfilled. Given that, it can be shown that the error of the resulting fit against the normal quadrature approach from chapter 2.2 is minimized. From (71) the calculation of τ i can be derived tan (2ω i · τ i ) = n sin(2ω i · t n ) n cos(2ω i · t n ) .
Next to that, the parameters R, I, C and S will be defined as follows: and Finally, the power spectral density P , the absolute amplitude A and the phase ϕ can be calculated from these coefficients. Hereby σ describes the standard deviation of the discrete data vector y n . In comparison to equation (15) from chapter 2.2, the equations above work quite similar. But now the modified function argument -extended by τ i , the correction terms C(ω i ) and S(ω i ) -leads to the optimal least square fit solution.
In the limit of an infinite number of equally spaced samples, the Lomb-Scargle estimator converges to the Fourier transform, so it becomes a consistent estimator (Mathias et al. 2004). The last statement is very important, because with an increasing number of samples the error reduces until the result converges to the true amplitude and phase.
To value the significance of the estimated amplitude the "false alarm probability" (FAP) is defined. Here the probability p, that there is no other larger amplitude P than P 0 , is expressed in terms of the exponential function above. For small values of p the approximation can be used. The free parameter M counts the independent frequencies in the data. These are difficult to measure a priori, but it turns out that with M = N/2 sufficient results can be achieved. A brief discussion on this issue can be read in the work of Townsend (2010, Cap. 6.2), Zechmeister and Kürster (2009)and Cumming et al. (1999).
Implementation The calculation of the phase is tricky. Invoking the arctan2() function is mandatory. Nevertheless, to prevent errors the phase ϕ is calculated by which take the normalized components I and R.
Next to that, if the suggestions made by Townsend (2010) are taken into account then the above algorithm can be shortened. The problem is that the straight forward implementation runs two times over the whole data set, while calculating τ prior the rest of the parameters. A keen refactoring of the equations above will minimize the computational cost. For a certain frequency ω i the amplitudes can be calculated out of the parameters in one single loop. A vectorized R-code example is given below. Here the data and the corresponding frequencies form the matrix omega.x, which is processed successively.
Application To show the power of the method let us first assume a simple example. The function is going to be sampled with N = 101 equally spaced samples. Remember the last point of the data vector equals the first point of the period of the signal. As discussed in the application section of chapter 3.1, the violation of the periodicity would lead to an error in the result of the Fourier transform. In addition to that, approximately 30% of the data where deleted.
Subsequently, Figure 16 shows the fragmented signal and the corresponding periodogram, which is calculated with the spec.lomb(x, y, f) function.
x <-seq(0,1,by=0.01) x.new <-seq(0,1,by=1e-3) yorg <-function(x) return(sin(2*pi*7*x)) l <-spec.lomb(x=x, y=y, f=seq(0,25,by=0.1)) lf <-filter.lomb(l,newx=x.new,phase="lin",threshold=3) The code above shows how Figure 16 can be created. Note that the frequency vector f contains 250 different frequencies to analyze the data. Compared to the single sided spectrum of the unfragmented data's analytic signal representation (black dots in fig. 16b) the Lomb-Scargle periodogram produces large side band amplitudes to the left and right of the main amplitude. This is quite typical if the data is non-uniformly sampled or even patchy. The dashed line represents the false alarm probability, which tends to zero if the corresponding amplitude is significant.
The spectral package also provides a filter.lomb() function, with which the most significant amplitudes can be extracted for reconstruction. Provided the continuous sampling vector x.new, the result is a new data set in which the remaining gaps are filled. A more complex example is given in Figure 17. Here the function y (x) = sin (2π · (x + N (0, ∆x/4))) + sin (2π · 20 · (x + N (0, ∆x/4))) + N (0, 1) consists of two sin-terms, the location of sampling jitters as well as the amplitude itself. In the simple test case here, normal distributed noise is assumed and the resulting data vectors x n and y n are N = 101 elements in length. However, for reconstruction the most significant amplitude values are selected again, so that the dashed line in panel (a) will fill the gaps correspondingly.
For further reading about non-uniformly sampled data the adaptive approach introduced by Stoica, Jian Li, and Hao He (2009) is recommended. This method has a better signal to noise ratio, but will have much more computational cost.

Conclusion
Today, spectral analysis is one of the key methods in data processing. It helps to filter signals, to accelerate calculations or to identify processes. While the complex mathematics is already documented in the given references and the common text books, the present work focuses on the application of the basic spectral methods in one dimension. Hereby working principles of the introduced methods and their possible misunderstandings are of interest. Small examples support the explanations about the work of operation and show the "do's and don'ts". The idea of the present work is to show how the complex underlying theory of spectral methods can be applied to solve real world tasks. A dedicated tool supporting this intent is the spectral package in the statistic environment R, which was developed contemporaneous to this paper. Further work will compose part II of this paper, which will cover the application of spectral methods in two and more dimensions.

Appendix
A.
The origin of f s /2 To derive the band limitation we have to assume a function y(t) for which the Fourier transform Y (ω) exist. There might be a natural maximum frequency ω max ≤ ω B in the signal, so the integration boundaries of the back transform can be limited to |ω B |. In the next step this function is sampled at discrete time instances t = n · T s , with a sampling frequency f s = 1/T s . It follows y(n · T s ) = 1 2π By substituting dω ↔ 2π df the modified equation (2) can be written as The pointer e iϕ is periodic with respect to 2π, so the span of ∆f /f s must not exceed unity to reject duplication. It follows ∆f f s ≤ 1 so that the maximum integration boundaries can be concluded in the first instance.
Due to the requirement that the frequency span must be less than the sampling frequency, the under-sampling of functions becomes possible. This changes the meaning of a typical Nyquist limit to a limitation of bandwidth which is necessary to obtain uniqueness.

B. Properties of the Fourier operator
In the following the properties of the Fourier operator are listed.

C. Derivation of the discrete Fourier-transform
The derivation of the discrete Fourier transform (DFT) is described below in detail. The already introduced and sketched out properties of the DFT will become clear with the discrete formulation of the Fourier transform. We will see that the spectral information becomes periodic and in the 1D case redundant.
First, let y be a real-valued function of time, with which is going to be sampled later on. For that the Dirac-function must be defined as a distribution as follows: The Dirac pulse is defined as δ(t = 0) = 0 for all values of t except t = 0. At this point the Dirac function returns an infinite value, δ(t = 0) → ∞. However, the behavior (3) turns the Dirac function into a distribution. The finite value of this integral enables subsequent calculations to sample a point from the continuous function y = x(t). This is called "sifting" property and can be written as In conjunction with the Fourier transform (4) leads to a necessary simplification for the derivation of the DFT.
Next to that the continuous function x(t) is going to be sampled with the help of the sampling function ∞ n=−∞ δ(t − nT s ). The sampling series evaluates x(t) at n reference points in time by multiplying the sampling function to x(t).
Because of the properties of the Dirac pulse the sum (5) converges to the value A(t) → ∞.
In the next step, the sampling series A(t) is going to be inserted in the Fourier integral This transforms the function x(t) from time or spatial domain into frequency domain, expressed by the complex valued spectrum X(ω). The capital letter X(ω) supports the result, which is a bijective mapping between x(t) ↔ X(ω). At this point both functions x(t) and its x (t) e −im·ω 0 ·t dt (7) with m ∈ N. Inserting the sampling series A(t) in (7), now discrete points enter the discrete spectrum: The distribution character (3) and the sifting property of the Dirac function help to rearrange the equation into the DFT form Bronštejn and Semendjaev (2001); Hoffmann (1998), x (nT s ) e −im·ω 0 ·nTs .
For a finite data set of N equidistant data points, ω 0 = 2π/T and a total length of T = N · T s , this can be written to: Here T s denotes the inverse sampling frequency T s = f −1 s and m ∈ N measures the normalized discrete frequencies on which the spectrum is being evaluated.
This discrete form of the Fourier transform provides several properties which define the spectrum X(m). 2. X(ω) can be interpreted as a line spectrum because it is only evaluated at ω = 2π N m discrete frequencies.
The main feature, which arises from (10), is the periodic spectrum. In consequence of that, frequencies which are larger than half of the sampling frequency are projected into the lower half of the spectrum. In case of a signal which is not band limited, an explicit mapping between a certain frequency and the underlying physical process is not possible. Although, an additional information is provided which enables a fold back. An example is given in chapter 3.1.
is not so obvious to resolve. The following derivation will explain this in detail. First, the exponent of the left hand side is extended with −ε · ω, which is then taken in the limit to zero.