Interpolated DFT for $\sin^{\alpha}(x)$ Windows

  • Published on
    25-Mar-2017

  • View
    214

  • Download
    2

Embed Size (px)

Transcript

  • 754 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 63, NO. 4, APRIL 2014

    Interpolated DFT for sin(x) WindowsKrzysztof Duda and Szymon Barczentewicz

    Abstract This paper describes interpolated discrete Fouriertransform (IpDFT) for parameter estimation of sinusoidal anddamped sinusoidal signals analyzed with a sin(x) window. For = 0, 2, 4, . . . sin(x) windows are RifeVincent class I (RVI)windows, for which IpDFT algorithms are known. We present anew IpDFTs for = 1, 3, 5, . . .. The bias-variance trade-off ofthe proposed IpDFT fits between results offered by RVI windows,e.g., for = 1, we get higher noise immunity than Hann (RVIorder 1, = 2) window and lower bias than rectangular (RVIorder 0, = 0) window.

    Index Terms Discrete Fourier transform (DFT), frequencyand damping estimation, frequency domain measurements, inter-polated DFT, signal processing, windowing.

    I. INTRODUCTION

    ESTIMATION of frequency, damping, amplitude, andphase of a single frequency and multifrequency signal isbased on the Fourier analysis or parametric modeling [1], [2].Tutorial comparisons of those methods are available, e.g., in[3] and [4].

    For nonperiodic signals, discrete Fourier transform (DFT)spectrum is affected by spectral leakage [5], [6], and estimatesof frequency, damping, amplitude, and phase are improvedby windowing and interpolated DFT (IpDFT) algorithms[7][17]. Typically, RifeVincent class I (RVI) windows (alsoknown as the maximum sidelobe decay windows) are used inIpDFT because it is easy to obtain IpDFT formulas for thosewindows. However, RVI windows are sin(x) windows for = 0, 2, 4, . . ..

    In this paper, we extend known IpDFTs for RVI windows tosin(x) = 0, 1, 2, 3, . . . windows. The contribution of thispaper is the new set of IpDFTs for sinusoidal and dampedsinusoidal signals analyzed with sin(x) = 1, 3, 5, . . .windows.

    II. SIGNAL MODEL AND sin(x) WINDOWS

    Let us consider a discrete time, damped sinusoidal signal inthe following form:

    x[n] = A cos(0n + )edn d 0 (1)where A > 0 is the signals amplitude, 0 < 0 < is itsangular frequency in radians, and 0 = rad corresponds to

    Manuscript received May 24, 2013; revised September 5, 2013; acceptedSeptember 6, 2013. Date of publication February 20, 2014; date of currentversion March 6, 2014. This work was supported by the Polish NationalScience Centre under Grant DEC-2012/05/B/ST7/01218. The Associate Editorcoordinating the review process was Dr. Dario Petri.

    The authors are with the Department of Measurement and Electronics,AGH University of Science and Technology, Krakw 30-059, Poland (e-mail:kduda@agh.edu.pl; barczent@agh.edu.pl).

    Color versions of one or more of the figures in this paper are availableonline at http://ieeexplore.ieee.org.

    Digital Object Identifier 10.1109/TIM.2013.2285795

    the half of the sampling rate in hertz, is the phase angle inradians, d 0 is the damping factor, n = 0, 1, 2, . . . , N 1is the index of the sample, and N is the number of samples.For d = 0 (1) simplifies to undamped sinusoid x[n] =A cos(0n+). IpDFT algorithms designed for damped signal(1) are still valid for undamped sinusoid, but the reversestatement is not true, because the damping changes the shapeof the signals spectrum [13].

    Before DFT analysis, the signal x[n] is multiplied by thetime window w[n]

    v[n] = w[n]x[n]. (2)DFT V [k] of signal v[n] is defined as

    V [k] =N1

    n=0v[n]e j 2N kn k = 0, 1, . . . , N 1. (3)

    Frequencies of DFT bins with indices k = 0, 1, 2, . . . , N 1are [k] = (2/N)k rad; thus we refer to the value 2 /N radas a DFT frequency step.

    The sin(x) windows are defined as [6]

    w[n] = sin(

    Nn)

    n = 0, 1, 2, . . . , N 1 (4)with normally being an integer.

    Using the well-known formulas for powers of trigonometricfunctions, see (5.72) and (5.70) in [18], definition (4) may berewritten separately for even and odd .

    For even = 2M (with M being nonnegative integer)window (4) is

    w=2M [n] =M

    m=0(1)ma2M [m]cos

    (2

    Nmn

    )(5)

    where

    a2M [0] = 122M

    (2MM

    )a2M [m] = 1

    22M1

    (2M

    M m)

    m = 1, 2, . . . , M. (6)Cosine windows defined by (5) are known as RVI windows

    of order M . For M = 0, = 0 we get a rectangular window,and for M = 1, = 2 we get Hann window.

    For odd = 2M 1 (with M being positive integer),window (4) is

    w=2M1[n]=M

    m=1(1)ma2M1[m]cos

    (2

    N(m0.5)n+

    2

    )

    (7)

    where

    a2M1[m] = 122M2

    (2M 1M m

    )m = 1, 2, . . . , M. (8)

    0018-9456 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

  • DUDA AND BARCZENTEWICZ: INTERPOLATED DFT 755

    TABLE I

    WINDOW COEFFICIENTS (6), (8)

    Window coefficients a[m] (6) and (8) for = 1, . . . , 8 aregiven in Table I.

    IpDFT algorithms described in the literature are designedfor the signal analyzed with RVI window (5).

    In the following, we derive new IpDFT algorithms forundamped and damped sinusoid analyzed with window definedby (7), i.e., for sin(x) window with odd .

    It is seen from (5) that for even , the time window isa weighted sum of rectangular window and M rectangularwindows modulated (i.e., shifted in frequency) by cosines withfrequencies being integer multiplies of a DFT frequency step.

    For odd , the time window (7) is a weighted sum of Mfrequency-shifted rectangular windows but the frequency shiftis always in the middle between the two successive DFT bins(i.e., the shift always contains the half of the DFT frequencystep) because there is (m0.5) in definition (7).

    The spectrum of (5) is

    W2M (ej) =

    M

    m=0

    [(1)m a2M [m]

    2W0(e

    j ( [m]))

    +(1)m a2M [m]2

    W0(ej (+ [m]))

    ](9)

    where [m] = (2/N)m rad, and W0(e j) is the spectrum ofthe rectangular window

    W0(ej) =

    N1

    n=0e j(N1)/2 sin(N/2)

    sin(/2). (10)

    The spectrum of (7) is

    W2M1(e j)

    = jM

    m=1

    [(1)m a2M1[m]

    2W0(e

    j ( [m]+/N))

    (1)m a2M1[m]2

    W0(ej (+ [m]/N))

    ]. (11)

    To derive IpDFT for the damped signal (1) we interpretthis signal as undamped sinusoid analyzed with the dampedwindow [13]

    v[n] = w [n]x [n] = w [n]A cos(0n + )edn= w [n]A cos(0n + ) (12)

    where w [n] = w [n]edn.

    The spectrum of damped window (5) is

    W2M (ej ) =

    M

    m=0

    [(1)m a2M [m]

    2W0(e

    j ( [m]))

    +(1)m a2M [m]2

    W0(ej (+ [m]))

    ](13)

    and the spectrum of damped window (7) is

    W2M1(e j )

    = jM

    m=1

    [(1)m a2M1[m]

    2W0(e

    j ( [m]+/N))

    (1)m a2M1[m]2

    W0(ej (+ [m]/N))

    ](14)

    where = jd and W0(e j ) is the spectrum of dampedrectangular window

    W0(ej ) = e j (N1)/2 sin(N/2)

    sin(/2). (15)

    In the Appendix, spectrum (14) is used for deriving newIpDFTs. For undamped signal = , and (13) and (14) reduceto (9) and (11).

    III. PROPOSED IPDFT ALGORITHMS

    In this section, we propose IpDFTs for sin(x) windowswith odd .

    The spectrum of the windowed sinusoidal signal v[n] =w[n]A cos(0n + ) is

    V (e j) = A2

    e jW (e j (0))+ A2

    e jW (e j (+0)) (16)

    where W (e j) is the spectrum of the used time windoww[n] (undamped or damped). If the signals frequency 0does not equal [k] = (2/N)k rad, then the IpDFTsare used for the accurate estimation of signals parame-ters [i.e., for finding the maximum of V (e j) (16) locatedbetween DFT bins, see Fig. 1]. Based on a few succes-sive DFT bins (3) of the windowed signal (2) frequencycorrection (displacement) is estimated that fulfils theequation

    0 = (kmax + )(2/N) 0.5 < 0.5 (17)

    where kmax is the index of the DFT bin with the highestmodulus (see Fig. 1).

    According to (16), the main lobe of amplitude spectrumof windowed sinusoidal signal depicted in Fig. 1 hasapproximately the same shape as amplitude spectrum of thewindow alone.

  • 756 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 63, NO. 4, APRIL 2014

    Fig. 1. Amplitude spectrum of windowed sinusoidal signal. Illustration ofDFT interpolation problem.

    Let us define the following relations:|V [kmax + 1]||V [kmax]| =

    |V ([kmax + 1])||V ([kmax])|

    = |V (0 2p2/N + 2/N)||V (0 2p2/N)| |W (2p2/N + 2/N)||W (2p2/N)| (18)

    |V [kmax]| + |V [kmax + 1]||V [kmax]| + |V [kmax 1]|= |V (0 3p2/N)| + |V (0 3p2/N + 2/N)||V (0 3p2/N)| + |V (0 3p2/N 2/N)| |W (3p2/N)| + |W (3p2/N + 2/N)||W (3p2/N)| + |W (3p2/N 2/N)| . (19)

    Solving (18) for 2p we obtain a two-point IpDFT (two DFTbins are used) for the signal analyzed with sin(x) window

    2p = sgn(|V [kmax + 1]| |V [kmax 1]|) (/2 + 1)|V [kmax 1]| (/2)|V [kmax]||V [kmax]| + |V [kmax 1]|

    = 0, 1, 2, 3, . . . (20)with |V [kmax + 1]| in (20) if |V [kmax + 1]| > |V [kmax 1]|(see Fig. 1) and |V [kmax 1]| otherwise, and sgn() denotingsign function that returns 1 for argument less than 0, 1 forargument greater than 0, and 0 for 0.

    Solving (19) for 3p we obtain a three-point IpDFT for thesignal analyzed with sin(x) window

    3p = sgn(|V [kmax + 1]| |V [kmax 1]|) |V [kmax + 1]| + |V [kmax 1]|2|V [kmax]| + ||V [kmax + 1]| |V [kmax 1]||

    = 0 (21a)3p = (/2 + 1) |V [kmax + 1]| |V [kmax 1]|

    2|V [kmax]| + |V [kmax 1]| + |V [kmax + 1]| = 1, 2, 3 . . . . (21b)

    Solutions (20) and (21) of (18) and (19) were only reportedfor RVI windows (5). In the Appendix, we prove that they arealso valid for windows defined by (7).

    For damped signal (1) analyzed with sin(x) window theIpDFT solution given in [13] is extended to the case of odd

    d = + 12

    R1 R2( + 2)R1 R2 R1 R2 (22)

    d1 = 2N

    ( + /2)2 R1( /2 1)2

    R1 1 = 0.5 (23)

    d2 = 2N

    ( /2)2 R2( + /2 + 1)2

    R2 1 = 0.5 (24)

    where = 0, 1, 2, 3, . . . and the ratios of DFT bins are definedby

    R1 = |V [kmax + 1]|2

    |V [kmax]|2

    R2 = |V [kmax 1]|2

    |V [kmax]|2 . (25)

    Equations (22)(25) were derived in [13] for RVI win-dows (5). In the Appendix, we prove that they are also validfor windows defined by (7).

    Systematic errors of damping estimation by (23) and (24)are practically the same. As suggested in [15], the mean valued = (d1 + d2)/2 may be used for damping estimation.

    IV. RESULTS

    In this section, we analyze the systematic errors and noiseimmunity of the proposed IpDFTs and compare it with theresults obtained for RVI windows. In all the figures, thewindow is specified by value. For even , the order ofRVI window is given, and for odd , it is stressed up that thenew IpDFT is used. For sinusoidal signal, a three-point IpDFT(21) was used. The lower subscript in E stands for estimatedfrequency.

    A. Theoretical Systematic Error

    For the case of a single undamped sinusoid expressions forerrors of estimation due to simultaneous spectral leakage andadditive noise are derived in [16] for even and multipointinterpolation. E.g. for Hann (w=2[n]) window and three-pointIpDFT, the maximum systematic error is [16]

    max2 =2(kmax + )||

    2kmax + (1

    2)(4 2)[(2kmax + )2 1][(2kmax + )2 4] . (26)

    Maximum systematic error of estimation for the three-point IpDFT is defined in [17]. For window (7), it is

    max2M1 = max |3p | = (/2 + 1)

    |W2M1(2kmax + + 1)| |W2M1(2kmax + 1)|2|W2M1()| + |W2M1( 1)| + |W2M1( + 1)|

    (27)

  • DUDA AND BARCZENTEWICZ: INTERPOLATED DFT 757

    Fig. 2. Theoretical maximum systematic errors of estimation for = 1, = 2 (Hann), and = 3.

    where is the true value, and windows spectrum may beapproximated from (A2) by

    W2M1(e j) sin( N/2 + /2)e j((N1)/2)

    M

    m=1

    (1)ma2M1[m](2m 1)(N2 m + 0.5

    ) (N2 + m 0.5

    ) (28)

    e.g., maximum systematic error for w=1[n] is given by

    max1 =

    32

    1(2kmax++1.5)(2kmax++0.5)

    1(2kmax+0.5)(2kmax+1.5)

    2(+0.5)(0.5)+

    1(0.5)(1.5)+

    1(+1.5)(+0.5)

    (29)

    Fig. 2 depicts the theoretical maximum systematic errorsof estimation for = 1, = 2 (Hann), and = 3.As expected, by increasing , we gain smaller systematicerrors. We also note that theoretical errors presented in Fig. 2are practically the same as the errors obtained during simula-tions described next.

    For even , sin(x) window is a linear combination ofcosines with frequencies being even multiplies of /N rad, i.e.,harmonics of 2 /N rad. For odd , sin(x) window is a linearcombination of sines with frequencies being odd multiplies of /N rad. It goes from the Fourier modulation theorem thatmultiplying coherently sampled signal by the window witheven shifts its spectrum by integer multiplies of 2 /N radto the neighboring DFT bins, and the signal is still coherentlysampled. For odd , the spectrum of the signal is shifted byodd multiplies of /N rad in the middle between neighboringDFT bins, and the signal is not coherently sampled any more;but if the frequency correction is = 0.5 the signal windowedwith odd window becomes coherently sampled. The abovediscussion explains negligible errors in Fig. 2 for even andodd .

    B. Simulations

    Figs. 36 depict systematic errors of frequency estimationfor undamped and damped sinusoid. Test signal has N = 512samples. The frequency of the test signal was changed from(2 /N) rad to 128 (2 /N) rad with the step 0.1 (2 /N) rad.

    Fig. 3. Systematic errors of frequency estimation for undamped sinusoidanalyzed with sin(x) = 0, . . . , 8 window. Three-point IpDFT (21).

    Fig. 4. Systematic errors of frequency estimation for undamped sinusoidanalyzed with sin(x) = 0, . . . , 8 window. Zoomed-in view of Fig. 3.

    For the given frequency value, the phase was changed from /2 rad to /2 rad with the step /20 and the highestabsolute value of the frequency estimation error was takenas a systematic error.

    The lower OX axis is labeled by frequency in radians,and the upper OX axis is scaled in frequencies of DFTbins k. The frequency of the kth DFT bin is k(2 /N) radand sinusoidal signal with frequency k contains k periods.For integer k, sinusoidal signal is coherently sampled andfrequency correction is zero.

    In Fig. 3, bold lines highlight envelops of the error showingthe worst case. The same errors are depicted in Fig. 4 withzoomed axes. It is seen from Fig. 4 that for even (i.e., forRVI windows) estimation error is practically zero for integer k,whereas for odd , the error is negligible for integer k plus 0.5.Results of simulations presented in Fig. 4 are practically thesame as theoretical errors shown in Fig. 2, i.e., max|E | max2M12/N .

    Figs 5 and 6 show systematic errors of frequency estimationfor damped sinusoid. Damped signal is not periodic, andeven for integer number of cycles the estimation error is notnegligible, although local minima are observed in Fig. 6.

    It is seen in Figs. 26 that systematic errors of the pro-posed IpDFTs are placed between errors obtained with RVIwindows for both undamped and damped sinusoid. Thus

  • 758 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 63, NO. 4, APRIL 2014

    Fig. 5. Systematic errors of frequency estimation for damped sinusoid (1),d = 0.001, analyzed with sin(x) = 0, . . . , 8 window. IpDFT (22).

    Fig. 6. Systematic errors of frequency estimation for damped sinusoid (1)analyzed with sin(x) = 0, . . . , 8 window. Zoomed-in view of Fig. 5.

    Fig. 7. Standard deviation of frequency estimation divided by CRLB forundamped sinusoid analyzed with sin(x) = 0, . . . , 8 window. Three-pointIpDFT (21).

    we gained possibility of additional adjustment of the IpDFTalgorithm.

    Figs. 7 and 8 show standard deviation of frequencyestimation for undamped and damped sinusoid respectively asa function of signal-to-noise ratio (S/N) divided by CramrRao Lower Bound...