Interpolation Algorithms of DFT for Parameters Estimation of Sinusoidal and Damped Sinusoidal Signals

Discrete Fourier Transform (DFT) is probably the most popular signal processing tool. Wide DFT use is partly dedicated to fast Fourier Transform (FFT) algorithms (Cooley & Tukey, 1965, Oppenheim et al., 1999, Lyons, 2004). DFT may also be efficiently computed by recursive algorithms in the window sliding by one sample (Jacobsen & Lyons, 2003, Duda, 2010). Unfortunately, DFT has two main drawbacks that deteriorate signal analysis which are (Harris, 1978, Oppenheim et al., 1999): 1) spectral leakage, and 2) sampling of the continuous spectrum of the discrete signal. Spectral leakage is reduced by proper time windows, and the frequency bins between DFT bins are computed by interpolated DFT (IpDFT) algorithms, thoroughly presented in this chapter.


Introduction
Discrete Fourier Transform (DFT) is probably the most popular signal processing tool. Wide DFT use is partly dedicated to fast Fourier Transform (FFT) algorithms (Cooley & Tukey, 1965, Oppenheim et al., 1999, Lyons, 2004. DFT may also be efficiently computed by recursive algorithms in the window sliding by one sample (Jacobsen & Lyons, 2003, Duda, 2010. Unfortunately, DFT has two main drawbacks that deteriorate signal analysis which are (Harris, 1978, Oppenheim et al., 1999: 1) spectral leakage, and 2) sampling of the continuous spectrum of the discrete signal. Spectral leakage is reduced by proper time windows, and the frequency bins between DFT bins are computed by interpolated DFT (IpDFT) algorithms, thoroughly presented in this chapter.

Basic theory 2.1 Signal model
IpDFT algorithms may be derived for discrete sinusoidal or damped sinusoidal signals. Discrete sinusoidal signal is defined as and discrete damped sinusoidal signal is defined as where A>0 is signal's amplitude, 0<ω 0 <π is signal's frequency in radians or radians per sample also referred as angular frequency or pulsation, and ω 0 =π rad corresponds to the half of the sampling rate F s in hertz, -π<φ≤π is the phase angle in radians, n is the index of the sample, N is the number of samples, and d is damping factor. If discrete signals (1) and (2) result from sampling analog counterparts then 00 2( / ) www.intechopen.com where F 0 is the frequency of analog signal, v(t)=Acos(2πF 0 t+φ) or v(t)=Acos(2πF 0 t+φ)e dFst , in hertz, F s is sampling frequency in hertz, and t is continuous time in seconds.
In section 3 it is shown how to estimate parameters of (1) and (2) i.e. A, ω 0 , φ and d with the use of DFT. If the investigation refers to analog counterpart signal than parameters of the discrete signal should be rescaled adequately, for example F 0 =F s ·ω 0 /(2π).

DFT analysis
Equations in this section are taken from the textbook (Oppenheim et al., 1999). Fourier transform (FT) of infinite length discrete time signal x n is defined as () jj n n n Xe xe where n is integer sample index that goes from minus to plus infinity and ω is continuous frequency in radians (angular frequency, pulsation). Continuous spectrum X(e jω ) defined by (4) is periodic with the period 2π. The notation X(e jω ), instead of X(ω), stresses up the connection between FT and Z transform.
For finite length discrete time signal v n containing N samples DFT is defined as Signal models (1) and (2) were given with rectangular window defined as Discrete signal is always analyzed with the time window as only finite number of samples may be read into computer. If no other window is purposely used then the signal is analyzed with rectangular window (7).
Based on FT convolution property signal windowing (6)  where X(e jΘ ) and W(e jΘ ) are the FT spectra of the infinite length signal x n and the time window w n . Thus, according to (8) we observe the convolution of the signal spectrum with the window spectrum, and not the signal spectrum alone.
FT of infinite length signal x n =cos(ω 0 n+φ) is a pair of impulses at frequencies ±ω 0 +2πk thus the spectrum of the windowed sinusoidal signal (1) Equation (9) is used as starting point in derivation of IpDFT algorithms. It is also used in leakage correction algorithms e.g. in (Radil et al., 2009).
According to (9) the spectrum V(e jω ) of the discrete, windowed, sinusoidal signal is the sum of two periodic replicas of window spectrum W(e jω ) shifted to the frequency ±ω 0 and rescaled by complex amplitude (A 0 /2)e ±jφ .  (1) with rectangular window and the modulus of its spectrum: a) signal (1), b) spectrum components for positive and negative frequencies (9), c) continuous FT spectrum (4)(9), d) sampling the continuous spectrum by DFT bins (5); ω E denotes estimated value Fig. 1 illustrates equations given above for sinusoidal signal (1) analyzed with rectangular window for A=1, ω 0 =1 rad, φ= 1.3 rad, and N =8. Fig. 1b depicts spectrum components for positive and negative frequencies (9). The sum of those components gives continuous spectrum shown in Fig.1c that may also be computed from the FT definition (4). As seen from Fig.1c the energy is not concentrated in the single frequency bin ω 0 (as would be for infinite length observation), but spills over to all neighboring frequencies. This phenomenon is called spectral leakage and may be the reason of significant estimation errors. In the given example spectral components for positive and negative frequencies (9) influence each other and maximum shown in Fig. 1c is moved from ω 0 =1 rad to ω E =1.04 rad (E stands for estimated value), that is estimation error equals 4%. Amplitude estimation error equals approx -15%. Estimation errors depend also on signal's phase φ as the sum (9) is complex. Because of the spectral leakage, the signal in the example disturbs its own spectrum. The impact of the spectral leakage would be stronger if the signal was a sum of sinusoidal signals, as every sinusoidal component would disturb its own spectrum and the spectrum of all others sinusoidal components. In DFT analysis spectral leakage is reduced by application of time windows other then rectangular and the error caused by the sampling of continuous spectrum only in frequencies ω k =(2π/N)k is practically eliminated by IpDFT algorithms.
DFT (5) is derived for periodic signals; as a consequence frequency analysis is correct only for the signals containing integer number of periods. The signal with integer number of periods is called coherently or synchronously sampled. In field measurements, for sinusoidal signals (1) close to coherent sampling is obtained with PLL (Phase Locked Loop) that keeps integer ratio of signal's frequency F 0 to sampling frequency F s , see (3). The frequency of coherently sampled signal is ω 0 =(2π/N)k rad and equals the frequency of DFT bin with index k. Damped sinusoidal signals (2) Fig. 2b shows continuous Fourier spectrum (4) and DFT bins of this signal. Upper OX axis is scaled in DFT index k, and lower OX axis is scaled in frequency in radians.
The range of frequencies intentionally exceeds 2π (one period) to stress up the periodic nature of the spectrum of discrete signals and the fact that spectral leakage for small frequencies originates from neighboring period (i.e. the part of the spectrum for negative frequencies). For coherent sampling of sinusoidal signal only one DFT bin is nonzero and the analysis is practically not affected by spectral leakage and sampling of the continuous spectrum. Fig. 3a depicts sinusoidal signal which does not contain integer numbers of periods. The frequency of this signal is ω 0 =2.2(2π/N)≈0.86 rad and lays between DFT bins with index k=2 and k=3. Fig. 3b shows FT spectrum and DFT spectrum of this signal. Estimation of signal's frequency based on the highest magnitude DFT bin is biased by the error of 0.79-0.86=-0.07 rad or -8%. High estimation errors would also be obtained for amplitude and phase estimation based on the highest magnitude DFT bin. Those errors caused by the sampling of continuous spectrum can be significantly reduced by IpDFT algorithms.

Time windows
The application of time windows in DFT analysis is thoroughly reviewed in (Harris, 1978), and also described in signal processing textbooks, e.g. (Oppenheim et al., 1999, Lyons, 2004. Other time windows with interesting properties are described in (Nuttall, 1981).
Window properties are determined by its spectrum W(e jω ) that consists from the main lobe which is the highest peak in the spectrum and side lobes. The shape of the spectrum of rectangular window (7) may be observed in Fig. 1b. The main lobe of the window should be as narrow as possible, and side lobes should be as low as possible. Narrow main lobe improves frequency resolution of DFT analysis, while low side lobes reduce spectral www.intechopen.com leakage. Rectangular window is the one with the narrowest main lobe, which is an advantage and the highest side lobes which is disadvantage. All the other time windows reduce side lobes, and thus spectral leakage, by the cost of widening main lobe i.e. reducing frequency resolution. It is also known that rectangular window has the best noise immunity although systematic errors caused by leakage may be dominant for signal containing small number of cycles.
Time windows are defined as cosine windows or non cosine windows. The cosine windows may be written in the form where w in A m w is introduced to distinguish from signal's amplitude in (1)(2).
IpDFT algorithms may only be analytically derived for Rife-Vincent class I (RVCI) windows. Coefficients A m w for RVCI windows are given in Tab.1. For M=0 RVCI window is rectangular window, and for M=1 RVCI window is Hanning (Hann) window. RVCI windows have the advantage of the fastest decay of the side lobes but they also have wide main lobe, which may be observed in Fig. 4. RVCI windows are also referred as cos α (X), =0,2,4,6... defined in (Harris, 1978 (10) Other examples of popular cosine windows are Hamming window (M=1, A 0 w =0.54, A 1 w =0.46), and Blackman window (M=2, A 0 w =0.42, A 1 w =0.5, A 2 w =0.08).
Cosine windows (10) are the sum of frequency modulated rectangular windows, thus the spectrum of order M cosine window is where ω m =(2π/N)m and W R (e jω ) is the spectrum of rectangular window Optimal due to main lobe width and side lobes level are non cosine Kaiser-Bessel and Dolph-Chebyshev windows (Harris, 1978). Those windows generally perform well as highlighted in (Harris, 1978 Fig. 12). Fig. 4 shows comparison of windows spectra for RVCI, Kaiser-Bessel and Dolph-Chebyshev windows for similar attenuation of the first side lobe approx -31.5 dB (Fig. 4a) and approx -101 dB (Fig. 4b). RVCI window has the fastest decay of side lobes and the widest main lobe. Kaiser-Bessel window has narrow main lobe and slowly decay side lobes, and Dolph-Chebyshev window has the narrowest main lobe and the side lobes on the same level. The markers in Fig. 4 denote DFT bins. Cosine RVCI window has M+1 nonzero DFT coefficients as stated by (11). For non cosine Kaiser-Bessel and Dolph-Chebyshev windows all DFT coefficients are non zero.

Interpolated DFT algorithms
IpDFT algorithm for sinusoidal signal analyzed with rectangular window was introduced in (Jain et al., 1979). In (Grandke, 1983) similar derivation was presented for sinusoidal signal analyzed with Hanning window. IpDFT algorithms for higher order RVCI windows are given in (Andria et al., 1989). In (Offelli & Petri, 1990) IpDFT algorithm for arbitrary cosine window was proposed based on polynomial approximation. In (Agrež, 2002) multipoint IpDFT was introduced with the feature of reducing long range leakage and thus reducing systematic estimation errors. IpDFT algorithm for the signal analyzed with arbitrary, even non cosine window was given in (Duda, 2011a).
IpDFT algorithms for damped sinusoidal signal analyzed with rectangular window are described in (Yoshida et al., 1981, Bertocco et al., 1994. In (Duda et al., 2011b) those algorithms were put in the same framework, and new algorithms for damped signal with rectangular window were proposed. Application of RVCI windows for damped signals in IpDFT algorithms was discussed in (Agrež, 2009, Duda et al., 2011b.

10
The IpDFT problem for sinusoidal and damped sinusoidal signals is depicted in Fig. 5 and may be formulated as follows. Based on the DFT spectrum V k (5) of the signal x n analyzed with the known window w n (6) find the frequency correction δ so to satisfy the equation where ω 0 is signal's frequency, N is the number of samples and k is the index of DFT bin with the highest magnitude. If |V k +1|>|V k -1|, as in Fig. 5, then there is '+' in (13).
For coherent sampling depicted in Fig. 2 frequency correction δ=0, and there is no need to use IpDFT algorithm. For the example of non coherent sampling shown in Fig. 3 frequency correction is δ=0.2 and it would be δ=-0.2 if the signal's frequency was ω 0 =1.8(2π/N)≈0.71 rad.
where it was used ω k =ω 0 -δ2π/N, ω k+1 =ω 0 -δ2π/N+2π/N that goes from Fig. 5 and (13). As stated by (9) the spectrum of the sinusoidal signal consists from the spectrum of the window moved from the position ω=0 to ω=ω 0 , thus the ratio of signal DFT (14) may be substituted by the ratio of window DFT The approximation sign '≈' is used in (15) instead of equality, because the ratios may slightly differ due to spectral leakage. Solving (15) for frequency correction δ we obtain two-point (2p) IpDFT formulas. Analytic solution of (15) is only possible for RVCI windows.
The ratio of DFT bins may also be defined with three bins as www.intechopen.com

Interpolation Algorithms of DFT for Parameters Estimation of Sinusoidal and Damped Sinusoidal Signals
which may be substituted by Solving (17) for δ we obtain three-point (3p) IpDFT algorithm.
Let us define the following summation where M is the order of cosine window (10) and A m w is the vector of window coefficients.
Then (15) and (17) In the next subsections IpDFT algorithms will be derived and explained. In all derivations the spectrum of the sinusoidal signal (9) is approximated by

Rectangular window (RVCI M=0)
Rectangular window w n R is defined by (7) and the spectrum of this window is given by (12).
By approximating sine functions by their arguments in (22) we have From (23) frequency correction is Signal's frequency in next computed from (13).
The amplitude of the frequency bin ω 0 may be found from the following proportion From (25) signal's amplitude is The phase of the frequency bin ω 0 may be found from the following equation From (27) signal's phase is The sign '+' or '-' in (28) is selected the same way as in (13).
In similar way, for the three-point interpolation we get from (17) and (12) and finally

Hanning (Hann, RVCI M=1) window
Periodic Hanning window is defined as 2 0.5 0.5cos , 0 , Hanning window (31) may be interpreted as the sum of rectangular window and frequency modulated rectangular window, thus based on FT properties, the spectrum of the Hanning window is the following sum of the spectra of rectangular windows Inserting (12) into (32), taking the approximation e j(π/N)(N-1) ≈-1+j π/N, and assuming π/N<<1 the modulus of Hanning window spectrum is www.intechopen.com By calculating modulus of DFT bins in (15) which may be rewritten as As seen in (35) coefficients by the powers of δ equal zero, and we get simple equation with the solution for frequency correction Signal's frequency is next computed from (13).
Similarly to (26), signal's amplitude is Signal's phase, computed as by (28) with the difference that in (28) frequency correction (24) is used and in (39) frequency correction (37) is used. For higher order cosine windows (M>1) we may write down (19) and demand that the coefficients by the powers of δ equal zero, as in (35) for Hanning window, and next solve the equations for δ. Described procedure allows us to find cosine window coefficients A m w that give analytic IpDFT solutions, and it turns out that by this procedure RVCI windows are found. For example, for M=2 grouping coefficient by the powers of δ gives

Higher order RVCI windows
which are RVC1 M=2 coefficients listed in Tab. I.
For the signal analyzed with RVCI window for two-point interpolation we get For M=0 and M=1 (41) agrees with previously derived formulas (24) and (37) for rectangular and Hanning window. Signal's frequency is next computed from (13). Signal's amplitude is For three-point interpolation from (20) For three-point interpolation and rectangular window frequency correction δ is computed from (30) and not (43). Equation (43) does not hold for rectangular window, i.e. M=0, because the spectrum of rectangular window contains only one nonzero DFT bin.
Signal's amplitude for three-point interpolation may be computed from the proportion 00 0 The phase of the frequency bin ω 0 for two-point and three-point interpolation is estimated the same way as for rectangular window (28) For known frequency ω 0 amplitude and phase of the signal may also be computed from the definition of Fourier transform (4) or in the least squares (LS) sense.

Sinusoidal signal -arbitrary windows
IpDFT formulas (41), (43) are only valid for RVCI windows, which have wide main lobe that deteriorates frequency resolution and noise performance. In practice it is often desired to analyze the signal with the window having better properties. It is known from literature, e.g. (Harris, 1978), that optimal parametric non cosine Kaiser-Bessel and Dolph-Chebyshev windows often perform superior over other windows including RVCI windows. In the following part IpDFT algorithm for arbitrary, even non cosine, windows is described for two-and three-point interpolation. The algorithm is based on polynomial approximation of the ratio R(δ) (15), (17) for the window of interest.
First, the ratio R(δ) (15), (17) is computed numerically for the selected window based on window spectrum, and then the dependence δ=f δ (R) is approximated by polynomial. During analysis the ratio R(δ) (15), (17) is evaluated from DFT bins, and frequency correction is estimated from previously computed δ=f δ (R).
where P d denotes L degree polynomial, and R is computed from window's spectrum. Signal's frequency is given by but this time R is computed from DFT bins of the analyzed signal.
Signal's amplitude is determined from the dependence X N =f X (δ), where for two-point interpolation and for three-point interpolation www.intechopen.com Signal's amplitude for two-point interpolation is and signal's amplitude for three-point interpolation is Signal's phase is computed the same way for two-point and three-point interpolation based on the dependence P N =f P (δ) Signal's phase is where P p is polynomial approximating P N =f P (δ), i.e. P N ≈P p (δ). Fig. 6 shows the dependence δ=f δ (R) for RVCI windows for two-point and three-point interpolation. Fig. 7 presents systematic errors of frequency estimation in dependence of the order of approximation polynomial for selected RVCI, Kaiser-Bessel and Dolph-Chebyshev windows for signal's frequencies ω 0 =0.05 rad and ω 0 =1 rad and signal's length N=512. Approximation polynomial was fitted to 64 points of window spectrum computed numerically via Fourier transform (4) in the frequency range from 0 to 2π/N rad, i.e. W(e -jω ) was computed by (4) for the set of frequencies from ω=0 to ω=2π/N rad with the increment (2π/N)/63 rad. Systematic errors were defined similar to (Schoukens et al., 1992), for each frequency ω 0 , test signals were generated with the phase from the interval <-π/2, π/2> changed with the step π/20 and the maximum absolute value of differences between estimated and true frequency was selected. It is seen from Fig. 7 that approximation polynomial of order 5 may give acceptable small systematic errors, nevertheless, in results shown in section 4 approximation polynomial of order 10 is used. For that order systematic errors, in Matlab 64 bit precision, for RVCI windows are practically the same for analytic IpDFT formulas (41), (43) and described approximation based IpDFT (49).

Damped sinusoidal signal -Bertocco-Yoshida algorithms
We start this section with derivation of DFT for damped sinusoidal signal (2). Let us rewrite signal (2) where ω k =(2π/N)k. Using the formula for the sum of the geometric series (58) becomes Going the same way for the second term in (57) In the following derivation of Bertocco-Yoshida IpDFT algorithms it is assumed that

Bertocco (BY-0) algorithm
Let us define the following ratio of complex DFT bins where V k is the DFT bin with the highest magnitude. From (63)  (65)

BY-1 algorithm
Let us define the ratio of the first order differences of the complex DFT bins in the form where V k is the DFT bin with the highest magnitude. Substituting (62) Damping and frequency are given by (65).

Yoshida (BY-2) algorithm
Let us define the ratio of the second order differences of the complex DFT bins in the form and λ is evaluated from (68). From (69) we get Damping and frequency are given by (65).
In the original derivation of Yoshida algorithm the ratio (69) is used and damping and frequency are given by (a) (b) Fig. 8. DFT of the damped sinusoidal signal. Solid circles denote DFT bins taken for Yoshida (BY-2) algorithm: a) ratio defined by (69), b) ratio defined by (71)

BY-3 algorithm
Let us define the ratio of the third order differences of the complex DFT bins in the form 2 2 21 1 11 2 Damping and frequency are given by (65).

Damped sinusoidal signal -RVCI windows
In the derivation of IpDFT algorithms for RVCI windows we treat damped signal with RVCI window as sinusoidal signal with damped window i.e.
where dn nn ww e − = is damped time window.

www.intechopen.com
Interpolation Algorithms of DFT for Parameters Estimation of Sinusoidal and Damped Sinusoidal Signals 21

Rectangular window (RVCI M=0)
Based on the spectrum of rectangular window (12) where jd ω ω =− . Let us define following ratios of the squares of frequency bins of the damped rectangular window spectrum where D=dN/(2π). In derivation of (79-80) sine functions were approximated by theirs arguments and it was assumed that (N-1)/N≈1.
By comparing D 2 in (79) and (80) we get desired frequency correction defined by (13) 12 12 1 2 1 22 The damping computed from (79) is and damping computed from (80) is For δ=0.5 (82) must not be used. In implementation if δ=0.5 then zero sample should be appended at the end of the signal to change δ. Equation (83) may always be used, as from definition (13) frequency correction is never equal -0.5.
Modulus of the frequency bin ω 0 may be computed from the proportion www.intechopen.com and signal's amplitude is Signal's phase may be computed as in case of sinusoidal signals (27 Sign '+' or '-' in (87) is taken the same way as in (13).

Hanning (Hann, RVCI M=1) window
The spectrum of damped Hanning window HH d n nn ww e − = is given by where We ω i s t h e s p e c t r u m o f d a m p e d r e c t a n g u l a r w i n d o w ( 7 7 ) a n d jd Inserting (77) where D=dN/(2π). From (90) and (91) we get desired frequency correction defined by (13) 12 12 1 2 3 24 2 The damping computed from (90) and (91) where D=dN/(2π) and M is the order of RVCI window. For M=0 (rectangular window) (97-98) becomes (79-80), and for M=1 (Hanning window) (97-98) is (90-91).

Some properties of IpDFT algorithms
In this section we present results of simulations that describe systematic errors and noise immunity of IpDFT methods. Because of space constrains, only the results of frequency or frequency and damping estimation are presented. Including results for amplitude and phase estimation would multiply the number of figures by three. Furthermore, in practice estimation of frequency and damping is of primary importance, and once having frequency and damping the amplitude and phase may be estimated by LS or FT. It is also true that amplitude and phase estimation errors, not shown in this section, behave similarly to frequency estimation errors.
First, systematic errors of IpDFT algorithms for sinusoidal and damped sinusoidal signals are presented, and then robustness against additive, zero-mean, Gaussian noise is shown. In all simulations number of samples N=512 was chosen.
Simulations were conducted in Matlab 64-bit floating point precision. Accuracy of this precision determined by the function eps (Matlab) is on the level 10 -15 -10 -16 , and estimation errors cannot be lower than this accuracy.

Systematic errors
In this section systematic errors of frequency estimation for sinusoidal signals and frequency and damping estimation for damped sinusoidal signals are presented. For each frequency ω 0 or damping d, test signals were generated with the phase from the interval <-π/2, π/2> changed with the step π/20 and the maximum absolute difference between estimated and true value was selected.
For obtaining general conclusions the frequency of the test signals was swept in the whole range from 0 to π rad. For easier interpretation the frequency of the test signal is also given in DFT index k. Fig. 9 shows two sinusoidal test signals with N=512 samples. The signal in Fig. 9a with frequency ω 0 =1.5·(2π/N) contains 1.5 periods, whereas signal in Fig. 9b with frequency ω 0 =249.5·(2π/N) rad contains 249.5 periods. Frequencies of those signals scaled in DFT index k are 1.5 and 249.5 and it means, that in the frequency spectrum those signals lie in the half way between DFT bins k=1 and k=2 and bins k=249 and k=250, respectively, and in both cases frequency correction δ (13) equals 0.5. The first signal is sampled approx 341 times per period and the second only approx 2.05 times per period.
In simulations the range of damping from d=0.0001 (Fig.10a) to d=0.01 (Fig.10b) is considered.
It is seen from Figs. 11-12 that 3p interpolation gives smaller systematic errors than 2p interpolation. Increasing the order of RVCI window results in significant reduction of systematic errors. High order RVCI M=6 window may give negligible small estimation errors as visible in Fig. 11 which is the effect of the fastest decay of the sidelobes. Still, due to wide main lobe of RVCI M=6 window, the signal has to contain sufficient number of periods. Systematic errors for Kaiser-Bessel and Dolph-Chebyshev windows are significantly higher then for RVCI M=6 window because sidelobes of those windows does not decay so fast. However, it is seen from Fig. 12 that for analysis of the short signal containing 2-4 periods it is advantageous to use, narrow main lobe, Dolph-Chebyshev 120 dB and Kaiser-Bessel =15.8 windows over RVCI windows. Local minima for integer values of k in Fig. 12 occur for coherent sampling and cosine windows. For the case of coherent sampling frequency correction is δ=0, DFT analysis is correct and there is no need to use IpDFT algorithms.
Figs. 13-14 present systematic errors of damped sinusoidal signal frequency and damping estimation for BY algorithms and IpDFT with RVCI windows. The frequency of the test signals in Fig. 13 was changed from ω 0 =1.5·(2π/N) to ω 0 =249.5·(2π/N) with the step 8·(2π/N), and damping was set to d=0.01 (compare Fig. 10b). In Fig. 14 the damping was swept from d=10 -4 to 10 -2 with equidistant steps in logarithmic scale, and the frequency was set to ω 0 =10.2·(2π/N) (i.e. δ=0.2). It is seen from Figs. 13-14 that by choosing high order RVCI windows significant reduction of systematic errors may be obtained, however, the price for this gain is the need for longer signal (in the sense of number of cycles) as explain previously for sinusoidal signals, and higher noise sensitivity as shown in the next section.

Noise
Noise performance is typically illustrated by comparison with Cramér-Rao Lower Bound (CRLB). Unbiased estimator that reaches CRLB is optimal Minimum Variance Unbiased (MVU) estimator. CRLB for sinusoidal signal (1) disturbed by zero-mean Gaussian noise with variance σ 2 is given by (Kay, 1993) www.intechopen.com and for damped sinusoidal signal (2) by (Yao & Pandit, 1995) where E in subscripts stands for estimated value, and η is signal to noise ratio defined as It is seen from (103-104) that the variance of frequency estimator is inverse proportional to the third power of signal length N. This strong dependence suggests high number of samples for frequency estimation, i.e. high sampling frequency or/and long observation time. Fig. 15 shows exemplary realization of the sinusoidal signal disturbed by the zero-mean Gaussian noise with S/N=10 dB (105) and its DFT spectrum. Mean value and standard deviation of estimation errors, shown in Figs. 16-18, were computed from 1000 realizations of the test signal. Signal's phase was generated as a random variable with uniform distribution on the interval <-π/2, π/2>. Fig. 16 presents results of frequency estimation of sinusoidal signal analyzed with different windows for three-point interpolation. It is seen from Fig. 16 that for S/N from approx 15 dB to 40 dB rectangular window (RVCI M=0) has the best noise immunity. For the lower disturbance (higher S/N) systematic errors become more significant than noise for rectangular window and better results are obtained with RVCI M=1 (Hanning) window, Kaiser-Bessel =4.86 window, and Dolph-Chebyshev 50 dB window. RVCI M=6 has the highest noise sensitivity due to the widest main lobe. Figs. 17-18 present results of frequency and damping estimation for damped sinusoidal signal. It is seen that for high S/N systematic error is dominant for BY-0, and RVCI M=0 window. BY-1 algorithm gives the best results in the wide range of S/N. It is seen from Fig. 18b that BY methods are better estimators of damping than IpDFT with RVCI windows.

Conclusion
This chapter describes DFT interpolation algorithms for parameters estimation of sinusoidal and damped sinusoidal signals. IpDFT algorithms have two main advantages: