• 沒有找到結果。

Estimation of single-tone signal frequency with special reference to a frequency-modulated continuous wave system

N/A
N/A
Protected

Academic year: 2021

Share "Estimation of single-tone signal frequency with special reference to a frequency-modulated continuous wave system"

Copied!
12
0
0

加載中.... (立即查看全文)

全文

(1)

This content has been downloaded from IOPscience. Please scroll down to see the full text.

Download details:

IP Address: 140.113.38.11

This content was downloaded on 28/04/2014 at 21:37

Please note that terms and conditions apply.

Estimation of single-tone signal frequency with special reference to a frequency-modulated

continuous wave system

View the table of contents for this issue, or go to the journal homepage for more 2012 Meas. Sci. Technol. 23 035002

(http://iopscience.iop.org/0957-0233/23/3/035002)

(2)

Chih-Fang Huang

, Hsiang-Pin Lu

and Wei-Hua Chieng

1Department of Information Communication, Yuan-Ze University, Taiwan, Republic of China 2Department of Mechanical Engineering, National Chiao Tung University, Taiwan, Republic of China

E-mail:jeffh@saturn.yzu.edu.tw

Received 1 September 2011, in final form 8 December 2011 Published 25 January 2012

Online atstacks.iop.org/MST/23/035002

Abstract

Estimating the frequency of a single-tone signal is a common problem in radar applications, such as frequency-modulated continuous wave radar. This study presents a frequency

estimation algorithm called the gradient search method using the derivative of discrete Fourier transform on the received signal samples. The analytical boundaries of the proposed method for different signal-to-noise ratios under the two conditions, with the rectangular window and with the Hann window, are derived. This study also compares the most appropriate algorithms available in the literature, including chirp-Z transform and other advanced methods.

Simulation and experimental results show that the proposed algorithm provides superior performance to previous methods.

Keywords:frequency estimation, FFT, DFT, Cramer–Rao bound (CRB), frequency-modulated continuous wave (FMCW)

(Some figures may appear in colour only in the online journal)

1. Introduction

Precisely estimating a sinusoidal signal frequency is an important task in signal processing. The frequency estimation problem is relevant to a wide range of areas, including radar, sonar and communications, and has consequently attracted considerable attention in the literature [1–4]. In liquid level measurements using a frequency-modulated continuous wave (FMCW) radar, the liquid level is converted from range data and is proportional to the frequency of received signals. Thus, the accuracy of frequency estimation is directly related to the accuracy of range data.

This study estimates the frequency of a single sinusoidal wave, which is given by

x(n) = a sin  2π fin fsn+ φ  + z(n), n= 0, 1, . . . , N − 1, (1)

where a, fin, fs and φ denote the signal amplitude, signal

frequency, sampling frequency and initial phase, respectively. The noise terms z(n) are assumed to be zero mean, additive white Gaussian noise (AWGN) with a variance ofσ2.

In a FMCW level measuring system, received signal errors are often produced by internal components, the environmental factor and quantization. The signal-to-noise ratio (SNR) is generally greater than 40 dB. This study analyzes the performance of each frequency estimator with an SNR of 0– 70 dB. Rife and Boorstyn [1] generated the Cramer–Rao lower bound and the maximum likelihood (ML) frequency estimator of this problem. The ML frequency estimator is a maximizer of the periodogram for calculation, for which the fast Fourier transform (FFT) and discrete Fourier transform (DFT) can be employed:

ˆfML= arg max

λ {Y (λ)} , (2a) 0957-0233/12/035002+11$33.00 1 © 2012 IOP Publishing Ltd Printed in the UK & the USA

(3)

where {Y (λ)} = N−1  k=0 x(k) e−j2πkλ    2 . (2b)

The ML Cramer–Rao bound (CRB) of frequency estimations can be treated as a fundamental (i.e. cannot be surpassed by any estimator) lower limit of the variance of any unbiased frequency estimator, operating on the signal modeled in equation (1). The CRB is given by

σ2 CR= 6 f2 s (2π )2ρN(N2− 1), (3) where SNRρ = a2/2σ2.

The length of the FFT needed to reach the CRB with the appropriate frequency resolution is too large for reasonable real-time implementation. The numerical maximization of (2) is never a computationally simple task and can lead to convergence and resolution problems [1]. Therefore, the frequency of a sinusoid is typically estimated using a two-step process involving a coarse estimator and a fine search. The coarse estimation stage is generally implemented using a maximum bin search to produce a coarse approximation of the periodogram [8,9]. This stage calculates the N-point FFT of the sampled signal and then locates the bin index with the largest magnitude.

Researchers have proposed a variety of fine frequency estimators [2–12]. Many multipoint interpolated DFT methods have been created. These methods are fast algorithms that interpolate the peak from DFT/FFT coefficients. In this approach, the main method of compensating the frequency is the ratio of DFT magnitudes for two or more frequency bins around the coarse estimation. Although Aboutanios and Mulgrew [8] developed a fast estimation method, their algorithm generates a high root-mean-square error (RMSE) when the value fin is close to 0 and fs/2. Moreover, this

algorithm cannot be combined with a windowing function, which generally increases the estimation accuracy [9]. The centrobaric energy method [9] is another interpolated DFT fine search algorithm that has been applied to the liquid level measuring system [13]. However, this algorithm is a one-shot scheme, that is, it cannot improve accuracy after several iterations.

Other fine search algorithms are based on the DFT magnitude comparison scheme. Following a few iterations, the estimated peak frequency is asymptotically close to the actual frequency. Zakharov and Tozer [5,6] developed simple algorithms consisting of iterative dichotomous searching procedures for precise signal frequency estimation.

The most prominent DFT-type frequency estimator is currently the chirp-Z transform (CZT) [11,12]. The CZT is based on an interpolation technique to resolve the spectrum between any two frequency components. This results in a higher frequency resolution than that of the original FFT method.

This study presents a novel search algorithm that resembles the dichotomous search algorithm developed by [5,6]. The algorithm performs the coarse-frequency estimation

using the position of the FFT peak. Subsequent iterations are carried out following the gradient search method (GSM). This study also analyzes and evaluates the RMSE of the GSM frequency estimation for different SNRs. The remainder of this paper is organized as follows. Section 2 presents the details of the proposed frequency estimator GSM including its underlying motivations. Section3analyzes the algorithms and assesses their asymptotic performance, and discusses convergence properties. Sections4and5present the simulation and experimental results, and section6 provides concluding remarks.

2. Gradient search method

The GSM consists of a coarse-frequency search followed by a fine estimation stage.

In the first step, we calculate the maximum bin from an initial coarse-frequency estimate. If the index of the maximum bin is denoted as m, the coarse peak frequency is fpeak= m f .

A search range [ flow, fhigh] can also be defined as flow = (m − 1) f and fhigh= (m + 1) f . Notably [8] showed that the search range can be refined in a narrower range of [ flow,

fhigh]= [(m − 0.5) f , (m + 0.5) f ].

In the second step, we derive the continuous DFT at fpeak.

Slope r indicates the curve direction.

In the third step, if r > 0, then the true peak lies to the right of fpeak. The flow can be moved to fpeak. If r < 0, then fhigh moves to fpeak. The new estimated peak, fpeak, is then ( flow+ fhigh)/2.

After the third step, we return to the second step to derive the continuous DFT at a new fpeak. Iterations may continue

until fpeakconverges to a defined accuracy. Figure1shows the

iteration steps of the GSM.

In the absence of noise, an observed sine wave, sampled at a known sampling frequency fs, is described as

x(n) = a sin(ωn + φ) ω = 2π fin

fs ∈ (0, π ), φ ∈ [0, 2π ), n = 0, 1, . . . , N − 1,

where a, fin, fs andφ are the amplitude, signal frequency,

sampling frequency and phase of the sine wave, respectively. The DFT is treated as a continuous Fourier transform to perform the differentiation on it. The N-point DFT of the samples, Xω, is defined as Xω≡ X(ω) = N−1  n=0 x(n)W (n) e−jωn= A ejθ, (4)

where X(ω) is the ω component of discrete Fourier transformation, A is the magnitude of X(ω), θ is the phase of X(ω) and W (n) is the windowing function.

The derivative of Xωis Xω = dXω dω = N−1  n=0 −jnxnW(n) e−jωn= Aωejθω+ jAωθω ejθω. (5) For a coarse peak frequency that is calculated using the FFT, the magnitude of DFT Xωis|Xω| = Aω. The gradient of|Xω|

(4)

Figure 1. Iteration process of the GSM.

of the gradient. To avoid resolving Aωdirectly, an alternative method was developed as follows. Dividing equation (5) by equation (4), we have Xω Xω = Aω Aω+ jθ  ω. (6) Thus, Aω= Re  Xω Xω  Aω,

where Re[·] denotes the operation of taking the real part. Because Aω> 0, the sign of Aωis identical to the sign of

Re

.

Let Xω= Rω+ jMω, Xω = Rω+ jMω. The gradient r is then given by

r= Re  Xω Xω  = Re[(Rω+ jMω)(Rω− jMω)] = RωRω+ MωMω. (7)

Equation (7) determines the direction of searching for the peak frequency. If r> 0, then the true peak lies on the right of fpeak,

whereas if r < 0, then the true peak lies on the left of fpeak. The steps of the GSM algorithm implementation are described as follows.

(a) Windowing samples are s(n) = x(n)w(n), n = 0, 1, . . . , N − 1, where w(n) is a Hann window; Ding [9] showed that applying a windowing function increases the accuracy of the interpolated DFT method.

(b) Given S= FFT(s) and Y (n) = |S(n)|, n = 0, 1, . . . , N − 1, find m = arg max {Y (n)}, coarse peak frequency

fm= m f , where  f = fs/N.

(c) The initial setting of peak frequency and search range are

fpeak= fm, flow = (m − 0.5) f , fhigh= (m + 0.5) f ;

For each iteration from 1 to Q iterations

Find derivative r at fpeak, where r= RfRf + MfMf;

if r> 0, then flow= fpeak;

else fhigh= fpeak; fpeak= ( flow+ fhigh)/2;

(5)

Figure 2. Iterative bounds of the GSM.

Figure 3. Estimation errors of GSM with different windowing functions.

3. Theoretical analysis

The results of this study show that the measurement error of the GSM comes from many sources, including the frequency leakage, iterative counts, AWGN sources, intrinsic error of GSM and quantization noise. The performance of the GSM is limited by these sources. Thus, this study derives the iterative bound, AWGN error bound and intrinsic error bound.

The iterative bound is dominated by the iterative counts Q. The frequency resolution of the coarse search is fs/N. After Q iterations, the final resolution is fs

2Q+2N. Thus, more Q iterations

can achieve a lower iterative bound, and the resulting error bound is closer to CRB. The iterative bound’s equation is

σQ=  fs 2Q+2 N 2 . (8)

In this paper, fs= 1000 Hz, N = 512, the resolution bandwidth

(RBW) is defined as RBW= fs

N = 1.953 Hz. The iterative

bounds of the simulation results (figure2) are consistent with the analytical formula given by (8). This study uses Q= 14 for comparison with the simulation results obtained by other algorithms.

(6)

Figure 4. Intrinsic bounds for GSM with and without windowing.

Figure 5. The simulated and analytical error bounds of GSM.

The question to ask at this stage is whether the variances of estimation can become any lower as the iteration Q increases. The answer is obviously impossible. Figure2shows that the RMSE is restricted by the iterative bound and the SNR, and AWGN influenced the measurement error. To estimate the accuracy of the GSM, it is necessary to derive the variance of

RfRf+ MfMf in (7).

Consider the sampled data with noise in the following form: ˜x(n) = x(n) + z(n) = a sin  2π fin fsn+ φ  + z(n), n= 0, . . . , N − 1, (9)

where ˜x(n) is the sampled signal and z(n) is WGN with the standard deviation σz. According to (9), the estimated

frequency ˆf is the best estimation while ˜Rˆf˜Rˆf+ ˜MˆfM˜ˆf = 0.

The tilde means that the item is sampled data in a noisy

environment. The DFTs of x(n) and z(n) are defined as

Xf = DFT(x(n)) = Rf+ jMf, Zf = DFT(z(n)) = Rf,Z+ jMf,Z.

The derivatives of Xf and Zf are defined as Xf = Rf+ jMf, Zf = Rf,Z+ jMf,Z. ˜Rˆf˜R ˆf+ ˜MˆfM˜ˆfcan be expressed as ˜Rˆf˜R ˆf+ ˜MˆfM˜ˆf = (Rˆf+ Rˆf,Z)(Rˆf+ Rˆf,Z) + (Mˆf+ Mˆf,Z)(M ˆf+ Mˆf,Z) = RˆfRˆf+ RˆfRˆf,Z+ Rˆf,ZRˆf+ Rˆf,ZRˆf,Z + MˆfMˆf+ MˆfMˆf,Z+ Mˆf,ZMˆf+ Mˆf,ZMˆf,Z. (10)

(7)

Figure 6. Plot of the RMSE of the frequency estimation for the GSM, and methods developed by Zakharov [5], Ding [9], Aboutanios [8] and CZT [11] under different SNRs.

Figure 7. Plot of the RMSE comparison of the GSM, Ding and Aboutanios methods over the full frequency range at SNR= 40 dB.

Because the WGN has a flat spectrum, Rf,ZRf,Z + Mf,ZMf,Z

is zero. Equation (10) then becomes

var( ˜Rˆf˜Rˆf+ ˜MˆfM˜ˆf) = R2ˆfvar(Rˆf,Z) + R2ˆfvar(Rˆf,Z) + M2

ˆfvar(Mˆf,Z) + M2ˆfvar(Mˆf,Z), (11)

where Rˆf,Z, Rˆf,Z, Mˆf,Z and Mˆf,Z are random variables, and they can be found by the definition of DFT and variance:

Zf = N−1  n=0 z(n)e−j2πfsfn var[Zf]= var(z(n)) · N−1  n=0 |e−j2πfsfn|2= Nσ2 z (12a) Zf = N−1  n=0 −j2π · n fs  z(n) e−j2πfsfn var[Zf]= var (z(n)) · N−1  n=0  −j2π · nfs · e−j2πfsfn 2 =2π2N(N − 1)(2N − 1) 3 f2 s σ2 z. (12b)

The variance of the real and imaginary parts is independent, and under the Gaussianity assumption, (12) can

(8)

(a) Double-precision is applied.

(b) Single-precision is applied.

Figure 8. Plot of the RMSE comparison of the GSM, CZT and Zakharov’s methods over the full frequency range at SNR= 70 dB. (a)

Double precision is applied. (b) Single precision is applied. be derived as var(Rˆf,Z), var(Rˆf,Z), var(Mˆf,Z) and var(Mˆf,Z) as follows: var[Rf]= var[Mf]= 2 z 2 , (13a) var[Rf]= var[Mf]= π2N(N − 1)(2N − 1) 3 f2 s σ2 z. (13b)

After substituting (13) into (11), the variance of ˜Rˆf˜Rˆf+ ˜MˆfM˜ˆf

becomes var( ˜Rˆf˜Rˆf+ ˜MˆfM˜ˆf) =  A2ˆfπ 2N(N − 1)(2N − 1) 3 f2 s +Xˆf2N 2  a2 2ρ. (14)

The estimation error ˆf− f0is related to ˜Rˆf˜Rˆf+ ˜MˆfM˜ˆf,

and the equation of estimation error can be derived as ˆf − f0= 1 2· ˜Rˆf˜R ˆf+ ˜MˆfM˜ˆf Aˆf · 1 Aϕ, (15)

where f0 is the true signal frequency, ˆf is the estimated

signal frequency, Af is the magnitude of X(ω), Af is the first

derivative of Af and Aϕis the second derivative of Af.

After some algebraic operations, the variance of estimation is var( ˆf − f0) ≈ 7π 4f2 s 96(4π − 16)2N3ρ = 0.6024 · fs2 N3ρ. (16)

(9)

(a)

(b)

Figure 9. Experimental settings. (a) Experimental apparatus. (b) Block diagram of experimental setting.

After derivation processing, the lower bound with the Hann window becomes

var( ˆf − f0)Hann≈ 1.2347 · f

2 s

N3ρ. (17)

Comparing (16) and (17) shows that the AWGN error bound with the Hann window seems larger than the one without the Hann window. However, the variance decreases as the SNR increases. Analysts normally apply a window to the data to alleviate the effects of leakage, but windowing also

distorts the original signals. The simulation results indicate that the windowing function improves the GSM accuracy at a higher SNR. The conditions of fs = 1000, N = 512, Q =

14 and an SNR of 0–60 dB are simulated. Figure3compares RMSE with different windowing functions, Hamming, Hann and Blackman Harris. In the case of low noise (SNR> 30 dB), processing with a windowing function can significantly reduce RMSE. Although these simulation results did not reach the CRB, the results obtained with a Hann window are better than

(10)

Figure 10. Experimental result of the RMSE comparison of GSM and Aboutanios’s algorithm.

Figure 11. Experimental results of RMSE comparison of GSM and algorithms of Zakharov and Ding.

those without windowing. Due to the obvious improvement in accuracy when using a Hann window at a high SNR, this study uses the Hann window for the following comparisons.

Another interesting characteristic for all estimators investigated in this study is how the estimation error would be if no noise exists. Simulation results show that even under noiseless conditions, there still are very low estimation errors for all estimators. These errors were labeled as the intrinsic bound in this paper. Intrinsic bounds are dependent on the windowing function and estimators, but unrelated to a noisy environment. Thus, it can be treated as the condition of a very high SNR. Figure 4 shows the simulation result of GSM’s intrinsic bound with and without Hann windowing.

The input frequency is fin ∈

0, fs

2

, and the lowest variance is at fin = f4s. The intrinsic bound with the

Hann window is lower than that without the Hann window. This is why the GSM with the Hann window shows better performance at a high SNR. Equations (18a) and (18b) are

the estimated intrinsic bounds for the rectangular and Hann windows, respectively, var( f − f0) = π 2f2 s N3  cos4D( f )(2πD( f ))2  

(rectangular window), (18a)

var( f − f0) = π 2f2 s N3 ×  cos2D( f )(2πD( f ))2  3 (Hann window), (18b) where D( f ) = f , f  N 4 =N 4 − f , f > N 4.

Combining the three error bounds, iterative bound, AWGN bound and intrinsic bound, figure5compares the simulation and analytical results for the GSM.

(11)

Figure 12. Comparison of simulation and experimental results for each algorithm.

4. Simulation results

Monte Carlo simulations were carried out to assess the performance of the proposed GSM frequency estimator. The performance of this method was compared with that of well-known methods, such as the iterative search method of Zakharov [5], the interpolated DFT of Ding [9], the iterative interpolation method of Aboutanios [8] and CZT [11,12]. All comparisons were accompanied by the benchmark ML CRB. Figure6presents the results.

The test parameters included a sampling rate of

fs = 1000 Hz and the FFT length of N = 512. For each

test point, the RMSE is the total error summation of simulated frequency f0 ∈

0, fs

2

. More than 2000 runs were simulated for each SNR. The iterative count Q is the same at Q= 14 for the Zakharov method and the proposed GSM.

Among these competitive methods, Ding’s and Aboutanios’s methods both adopt interpolation schemes that require few calculation resources. The disadvantage of these methods is that their RMSE is higher than that of the other estimators under low noise conditions. Figure7compares the RMSE of the GSM, Ding’s and Aboutanios’s algorithms at the SNR of 40 dB.

For an SNR exceeding 30 dB, figure 6 shows that the GSM, CZT and Zakharov’s methods are the three best estimators. It is interesting to investigating the performance of these estimators at a higher SNR. Figure 8(a) shows the simulation result at SNR = 70 dB, clearly showing that the GSM has the lowest RMSE. Note that the CZT and Zakharov methods require more precise computing capability

of calculating the amplitudes approaching the peak frequency. Experimental results (figure8(b)) show that using the data-type single precision (float) instead of double precision (double) degrades the performances of the CZT and Zakharov methods, though the GSM remains unaffected.

5. Experiments

This experiment compares the GSM with three other frequency estimators proposed by Zakharov [5], Ding [9] and Aboutanios [8]. This experiment was performed using an Agilent 33120A function generator as a signal source to simulate the beat frequency of a FMCW sensor. A low-cost micro control unit (MCU), Microchip dsPIC33F, served as the frequency measuring set. A personal computer (PC) was used to record the measuring data through serial communication with dsPIC33F. Figure9shows a block diagram of the experimental setting.

The Agilent 33120A facilitates changing frequency and magnitude. The testing signal frequencies ranged from 10 to 490 Hz in 10 Hz steps. For each frequency change, 50 measurements were carried out by all four algorithms. The measured frequencies were sent to a PC to analyze and record the RMSE. The dsPIC33 has a 12-bit internal analog-to-digital converter. The sampling rate was set to 1000 Hz. Because of being constrained by MCU memory, the maximum FFT length could only be 512 points.

The experimental results in figure 10 show that Aboutanios’s algorithm is highly frequency dependent. The

(12)

the GSM and that of other well-known frequency estimators were compared over the possible frequency range. Simulation results show that the proposed method is better than the interpolated DFT methods when the SNR exceeds 30 dB. Compared to the CZT and Zakharov dichotomous search method, the GSM requires less calculation precision capability. This study also shows that the GSM for the range measurement can be easily implemented in digital signal processors and low-cost MCUs, making it suitable for applications that require high accuracy and reliable level measurement.

Acknowledgments

The author would like to thank the National Science Council of the Republic of China, Taiwan, for financially supporting this research under contract no NSC 98-2221-E-009-047 and NSC 99-2410-H-155-035-MY2. U&U Engineering Inc. is commended for supporting experimental work, as well as Wen-Hsiang Lo and Yu-Chieh Lin for their assistance during the experiments.

frequency estimators with narrow acquisition range Proc.

Inst. Electr. Eng.—Commun. 148 1–7

[7] Aboutanios E 2004 A modified dichotomous search frequency estimator IEEE Signal Process. Lett. 11 186–8

[8] Aboutanios E and Mulgrew B 2005 Iterative frequency estimation by interpolation on Fourier coefficients IEEE

Trans. Signal Process. 53 1237–42

[9] Ding K and Jiang L 2001 Energy centrobaric correction method for discrete spectrum J. Vib. Eng. 14 354–8 [10] Huang Y and Huang Z 2008 A highly accurate iteratively

interpolated frequency estimator 4th Int. Conf. Wireless

Communications, Networking and Mobile Computing, 2008: WiCOM ‘08 pp 1–5

[11] Rabiner L R, Shafer R W and Rader C M 1969 The chirp-Z transform algorithm IEEE Trans. Audio Electroacoust.

Au-17 86–92

[12] Bluestein L I 1970 A linear filtering approach to the computation of the discrete Fourier transform Northeast

Electronics Research and Engineering Meeting Rec.

Au-10 451–5

[13] Liu J, Chen X and Zhang Z 2006 A novel algorithm in the FMCW microwave liquid level measuring system Meas.

數據

Figure 1. Iteration process of the GSM.
Figure 2. Iterative bounds of the GSM.
Figure 5. The simulated and analytical error bounds of GSM.
Figure 6. Plot of the RMSE of the frequency estimation for the GSM, and methods developed by Zakharov [ 5 ], Ding [ 9 ], Aboutanios [ 8 ] and CZT [ 11 ] under different SNRs.
+5

參考文獻

相關文件

Reading Task 6: Genre Structure and Language Features. • Now let’s look at how language features (e.g. sentence patterns) are connected to the structure

Wang, Solving pseudomonotone variational inequalities and pseudocon- vex optimization problems using the projection neural network, IEEE Transactions on Neural Networks 17

volume suppressed mass: (TeV) 2 /M P ∼ 10 −4 eV → mm range can be experimentally tested for any number of extra dimensions - Light U(1) gauge bosons: no derivative couplings. =&gt;

Define instead the imaginary.. potential, magnetic field, lattice…) Dirac-BdG Hamiltonian:. with small, and matrix

incapable to extract any quantities from QCD, nor to tackle the most interesting physics, namely, the spontaneously chiral symmetry breaking and the color confinement.. 

• Formation of massive primordial stars as origin of objects in the early universe. • Supernova explosions might be visible to the most

(Another example of close harmony is the four-bar unaccompanied vocal introduction to “Paperback Writer”, a somewhat later Beatles song.) Overall, Lennon’s and McCartney’s

Microphone and 600 ohm line conduits shall be mechanically and electrically connected to receptacle boxes and electrically grounded to the audio system ground point.. Lines in