Copyri ghl © IFAC Idel1l ificali oll and SySlclll Paralllelcr Eslimalioll. Budapesl. Hungary 199 1
A REAL-TIME PROCEDURE FOR HIGH-RESOLUTION MUL TIPLE-TONE DETECTION M. Bcrtocco, C. Offclli and D. Petri Diparcimcnto di Elcllronica e In/o rmalica. Universild di Padova. Via Cradenigo 61A. 135131 Padova. ltaly
Abstract The paper describes the use of first-order Discrete Prolate Spheroidal Sequences (DPSS) in the implementation of high resolution multiple tones detectors in real-time spectral analysis. The incoming samples are weighted by a DPSS and the Discrete Time Fourier Transform (DTFT) is evaluated. While in the case of a single line the windowing operation yields a spectrum exhibiting a sharp null at the frequency of the spectral line, no null is observed in the spectrum of multiple lines closely spaced in frequency. This different beha vior provided by the windowing allows an easy real-time discrimination between the two situations. Due to the small sidelobe levels of the adopted windows, multiple tones sufficiently spaced in frequency can be thought not interacting so that the application of the proposed procedure gives good results also in these cases. Kevwords: Spectral analyzers; Signal detection; Estimation; processing; Fourier analysis; Real-time computer systems.
Digital
signal
of resolving two equal level spectral lines closely spaced in frequency. This severe limitation is highlighted in Fig. 1, in which the spectrum of a two-tone signal is represented. The contribution of each tone is roughly summarized with a dotted line, while the resulting transform envelope is indicated wi th a solid line. If the spectral lines are sufficiently spaced in frequency, two peaks are present (Fig. la); on the contrary, when they are closely spaced a unique peak is observed (Fig. 1b). There is a minimum in the frequency distance of two adjacent tones that allows their theoretical detection; nevertheless, for currently available automatic methods based on a Fourier transform the actual distance is rather high.
INTRODUCTION In real-time spectral analysis of a sampled deterministic process composed of pure tones, the frequency-domain approach via fast Fourier transform algorithms is of great interest both for its high computational efficiency and for its robustness with respect to model mismatches. In this approach, an iterative scheme is commonly adopted in which it is given a finite data record of the input signal sampled at liT frequency, a suitable processing of the acquired data is performed in order to extract some useful information, and then a new record consecutive to the previous one is considered. The finite length of the input records is often determined by the processing system which physically performs the elaboration; it is also limited by the maximum delay allowed before a result is given.
Two Loca 1 Hax Ima
In spectral analysis, the finite length of the data record introduces several drawbacks; a first disadvantage is due to the implicit windowing of the analyzed signal with a rectangular window, which produces the "spectral leakage effect". One of the consequences of the leakage is that strong spectral lines often mask weak signal tones. To overcome this problem windows other than the rectangular one are commonly used to weight the time sequence before computing the tr a nsform (Harris,1978; Nuttall,1981). Another disadvantage consists in the poor ability 1001
(a) Non Resolvable Tones
Fig. 1. spectral tones.
Cb) Resolvable Tones
resolution
of
nearby
the selection of the "best" window will be introduced, and computational effort issues will be addressed. Finally some selected simulation examples will be discussed.
In the paper, a simple Fourier-based spectral estimate S () is obtained from the weighted samples of the i nput signal x(), i.e.
METHOD OUTLINE
Sx (i\.) (1)
1 tIT v[n] x(nT) e- j2rr i\.n/N 12
DPSSs Definition and Their Fundamental Properties for the Detection Problem. , i\.
E
[O,N] The DPSSs have been defined by Slepian (1978) for a given window length N. An analytical closed form for these sequences is not available. In the literature they are not defined in a straightforward form, but rather as the N solutions of the following eigenvector problem
n=O
where i\.=fNT is the normalized frequency expressed in bin, NT is the record length, and v[] is a weighting window function. The aim of this work is to detect the possibility that somewhere in the frequency axis two or more tone closely spaced coexist, against the case of a single spectral line. In order to achieve this goal, the window v [ ] is selected in such a way that remarkably changes in the spectral estimate Sx() are observed with respect to the two different situations.
N- I
sin[2rr(n-m)W] -------vh[n] rr(n-m) n =O
L
(2)
n,h = 0, .. ,N-1 that satisfies the conditions, for h=0, ... ,N-1:
normalization
N -I
L v~[n]
Among the possible windowing functions, in the paper it is proposed the use of the first-order Discrete Prolate Spheroidal Sequence (DPSS) as windows suitable for the detection of multiple lines in a known frequency interval. The utilized windows have an odd-symmetry in the time-domain, and hence the data are multiplied by positive and negative values. This implies that, when the analyzed sequence consists of a single tone, the corresponding spectrum exhibits a deep null at the input signal frequency. On the contrary, if two or more complex sinusoids closely spaced are present, the resulting transform drastically changes, and no deep null is found, allowing a clear distinction of the latter case from the former.
=
1
(3)
n =0
N -I
"\' v [n]
Lh n =O
N -I
>00 ,
L(
N;l - n)vh[n] >0 0
(4)
n ::: 0
and they are ordered on the basis of the eigenvalue magnitudes: (5)
From the numerous simulations carried out it has been found that the best results for the detection are achieved if h = 1 is chosen; besides the W parameter which affects the window shape gives optimal results when it is in the range from 2 to 5.
The adopted windows are characterized by very small sidelobe levels. This means that they also reduce the long range spectral leakage and allow to perform simultaneous tests on the presence of multiple lines in different frequency ranges. Furthermore, the evaluated spectral estimate gives information on the rough location of tones in terms of the S () peaks, regardless the x number of spectral lines which compose the input signal.
In order to give an idea of the window behavior, in Fig. 2 it is shown the magnitude transform, IV I () I, of VI [], for N = 64, and W = 4; the plot is normalized with respect to its maximum. The whole shape of IV I () I depends on w, in particular
A similar approach is discussed by Griffiths (1986), where the use of a linear FM chirp window has been proposed. A disadvantage of such a kind of weighting functions is due to the great amount of spectral leakage. As a consequence, when signals with multiple lines are analyzed, quite similar spectra are found regardless the tones are well spaced in frequency or not. Moreover, no information is given on the rough location of the frequency lines, thus requiring a further use of a classical taper. no rma lized frequen cy (bin)
In the following sections at first the DPSS definition will be recalled, and some fundamental properties of the DPSSs with respect to the detection problem will be highlighted. Then some figures of merit for
Fig. 2. Module of the Transform of the DPSS VI [n]; N=64 W=4.
1002
the mainlobe width is approximately equal to 2W; while the sidelobe levels tend to increase inversely with W. When the leakage is not a critical issue, low values of W can be chosen, otherwise good results were observed with W taken approximately equal or greater than 4. For the selected value W in Fig. 1, it is worth noticing the great amount of sidelobe attenuation, near to 80 dB, which guarantees a great immunity with respect to the interference of tones located at a distance greater than the mainlobe width.
Q)
'"0
....Z c:
00
ro
E
Three properties of the DPSSs, particularly useful for the detection algorithm development, are reported. The first is the symmetry property, in general expressed as : h,n=O, .. ,N-l
(6)
so that for odd h Vh () is zero for A = O. Furthermore, a DPSS transform of order h admits h zeros in (-W,W); this means that Vi (A) has no other null inside the window main lobe , except for A=O. The third important property is that the DPSSs are the orthogonal index limited sequences with the greatest fractional energy in (-W, W) , and consequently with the minimum energy outside the mainlobe (Slepian, 1978). These last two properties constitute the main reason for the choice of the DPSSs with h = 1 and W e [2,5] for the detection algorithm design. In the following, our discussion will be on the prolate sequence with h = 1, later on simply named vC], dropping the index h = 1.
no rmalized fr equency (bin)
Fig. 3. Spectrum of W=2,3,4,5 origin
vel , and
IV() I, for near A
N=64, the
an integer multiple k of the tone period To' the corresponding discrete spectrum S () exhibits a deep null for A = k. If the x T = kT condition is not satisfied, because w 0 of frequency quantization, the deep null of S () is replaced by a local minimum, and x the higher is the frequency resolution, the lower is the minimum level. In order to give a quantitative meaning to this effect it is advantageous to fix both a suitable threshold A, below which a local minimum is accepted as a nu 11, and a frequency resolution II at which spectral samples are evaluated. It is supposed that the single-tone signal: (7) x(nT) = eXp{j2ITA on/N)
In this paper, the Wand N parameters which determine the v h [ ] sequences are regarded as design quantities for the detection algorithm, and hence it will be assumed that the v[] functions are tabulated for some significant values of the parameters and stored somewhere for a subsequent use. In any case, reliable numerical algorithms for the v[] evaluation have been proposed for every value of Nand W of practical interest; see (Thompson, 1982) or also (Greenhall, 1990) for further details.
is given, where AO= T/T o ' According to the relationship (1), the spectral estimate becomes, S (A) x
IV(A-A ) 12 o
(8)
is actually evaluated at the normalized frequencies Ak= kll, with k an integer number between 0 and the smallest integer greater than 1/ ll. The worst case for the local minimum is observed when the tone is located in a middle point between two successive frequency samples (AO = mll + 1l/2, me71). Supposing that with a given II frequency resolution the minimum has a relative amplitude, with respect to the window peak equal to A, as a consequence for any other value of AO the resul ting relative minimum is deeper than A. For this reason A can be chosen as the threshold at which a relative minimum is detected as a window nUll.
Windows Selection .
As the main window characteristic for the detection is the presence of a null at the frequency origin, it is better to take a deeper investigation of the window behavior near A=O. The remarkably sharpness of this deep null is better highlighted in Fig. 3, in which a logarithmic scale has been adopted also for the frequency axis. Four curves have been plotted accordingly to four different values of W: 2, 3, 4, and 5. Each transform is normalized in order to have a unitary maximum. For numerical analysis it is required that the normalized A frequency A is quanti zed and a finite number of values is considered; we in this paper we call frequency resolution (ll) the distance between two successive frequency samples. When the input signal consists of a single tone and the observation interval Tw= NT is
In order to give a quantitative meaning to the proposed method, some quantities will now be defined; they will determine the window detection performances as well as the required frequency resolution. The quoted threshold A in this paper is named null confidence level. Another index of merit is the null bandwidth BA' defined as
1003
twice the distance from the origin to a point of IV() I A dB above the window maximum. The relationship between the two parameters is expressed by the relationship
Table 1. Null bandwidth BA o f
IV()I ex pressed in
bin, for diffe r ent values of th e null c onfidence value A, and W.
A
20 log la [
(9)
From inspection of Fig. 3, after a null conf idence level A has been chosen, the corresponding null bandwidth BA is determined. The required frequency resolution at which S () must be evaluated x is equal or less than BA' A standard FFT algorithm based on N points gives a normalized frequency resolution of 1 bin, corresponding to l/NT Hz. Greater resolution values can be achieved for instance if the windowed N-points data record vel x() is padded with M-N zeros to form an M-points data record. Choosing M = LN, in such a way that M is a power of two integer number, the corresponding obtained normalized frequency resolution becomes /',. = l/L bin equal to l/NLT HZ; the L quantity in this context is called zero padding factor. The theoretical freedom in the choice of L allows to obtain whichever value for the resolution /',.. Because the computational effort of a FFT algorithm is proportional with M log M, it is interesting to deal with values of L as small as possible. In any case, for high values of L an FFT algorithm can still be exploited but in order to limit the number M = LN of frequency samples involved, it is better to observe at one time only restricted intervals of the frequency axis.
~(dB)
10
ZO
Z.O 3.0 4.0 5.0
3.Z4E-1 3.80E-1 4. 44E-1 4. 44E-1
9. 34E-Z 1. 09E-1 1. Z7E- 1 1. 48E-1
30 3. 14E-Z 3. 66E-Z 4. Z8E- Z 4. Z8E-Z
40 9. 0ZE-3 1.Z3E-Z 1. ZsE-Z 1. 43E- Z
Offelli, Petri, 1990). Once A is given by the relationship L = 1//',. '" l/B A, the zero padding factor L can be derived from inspection of table 1. A compromise solution between better detectability degree (great values of A) and computational effort (great values of L) has to be carried out. The choice of A and L can eventually lead to reconsider the initial value of W, and then to adjust the couple A, L again.
SIMULATION RESULTS
Chosen a suitable value for the null confidence level A, and found the corresponding null bandwidth, a lower bound for the frequency resolution /',. is expressed by the obvious relation /',. ~ BA' Finally, once fixed /',., the zero-padding factor L is the greatest integer value which best approximates the ratio 1//',.. While a rough estimate of BA' and hence /',., can be obtained from inspection of Fig. 3, quantitative results are summarized in table 1. This table shows that BA increases as A and W increase; in fact for N sufficiently large the DPSSs can be regarded as a sampled version of the Continuous-time Prolate Wave Functions (Slepian, 1978). Hence, IV() I changes only for a normalization factor, which has no consequence as long as relative characteristics, e.g. bandwidth, are measured. It is observed that, from a practical point of view, the choice of the design parameters for the detection may be as follows. It is assumed a given value for the data record length N, usually imposed by real-time constraints, such as the maximum admitted delay for which an estimate is obtained. An initial value for W can be chosen on the basis of the desired window mainlobe width and window sidelobe level (Thompson, 1982; Harris, 1978;
In this section some examples which highlight the performances of the proposed method are presented. As the main interest is in the effect of the proposed weighting, no other processing is adopted, even if it is clear that further manipulations may improve the overall results. In each simUlation a data record length of N = 64 samples is chosen, a value often found in the literature; the bandwidth W is taken equal to 4, thus obtaining good results both with respects to the resolution and the long range spectral leakage effect. It has been fixed a sufficiently high null confidence level, i.e. at least 10, or 20 dB. As it can be seen from inspection of table 1, for W = 4 it is required a null bandwidth BAE< 0.12 bin (or 0.127/64 E< 3 1. 9 .10- Hz) , corresponding to a zero padding factor L = 8 '" l/B A; this requires the evaluation of an FFT based on M = 512 points. In Fig. 4 it is shown the windowed transform of a unitary sinusoid at the normalized frequency Aa= 10+2/L bin (i.e. 0.15674/T Hz) i since Aa is located in the middle of two consecutive FFT samples, this can be regarded as a worst case, in the sense that the minimum of the evaluated spectrum will exhibit the maximum possible magnitude. The evaluated transform is equal to the translated version at A of the o window V(), as expected; the null is more than 20 dB below the mainlobe maximum, according to table 1. In Fig. 4 it is also well shown the little influence of the window at frequencies far from Aa.
1004
Finally, in Fig. 6 the results due to the same test signal of Fig. 5, plus a third unit level tone at A 2 = 26 bin are plotted. The effect of the additional line is manifested in the lobe centered at 26 bin; it presents a sharp null, thus indicating a single line. This example also shows how the adopted windows make negligible the interference between the two lower frequency tones and the third one. Furthermore a rough estimate of the single tone location and the closed tones can be deduced.
norma l ized frequen cy (bin)
Fig. 4.
Spectrum signal.
of
a
single
tone
input
In Fig. 5 the effect of windowing is represented for a test signal consisting of two spectral lines with unitary amplitude and frequencies Aa= 10 and A1 = 10.5 bin respectively. The initial relative phase between the two tones is 0.2 9rr rad. The comparison of Fig. 4 and Fig. 5, clearly shows the difference of the evaluated spectra: the sharp null of Fig. 4 is not present in Fig. 5 and it is substituted by a flat portion, so that a multiple tone signal in the 10 bin neighborhood is detected with no doubt. As a consequence, a very simple post-processing technique applied to the S () estimate can x distinguish in an effective and reliable manner the latter case from the former. Besides, the whole detection process requires a computational effort of the same order of a single FFT, and thus it is very attractive for implementation in real-time systems.
en
'U Q)
'U
"
...., .~
0:
00
ro
E
mtlVmIIlVIIVVIIVlmll'~~ 56
60
mul tiple
Fig. 6. Spectrum of the input signal of Fig. 5 plus an additional tone well-spaced in frequency.
CONCLUDING REMARKS It has been proposed the use of the DPSSs of order one for the detection of multiple tones in an observed signal. A null in the origin of the transform of these windows allows to yield considerable changes in the evaluated spectrum when the input signal changes from one to two or more closely spaced sinusoids. The low sidelobe level makes effective the use of such windows in the analysis of a signal with a more complex harmonic structure because spectral lines sufficiently spaced in frequency have li ttle interference. One possible application would be to obtain a preliminary information on the harmonic contents of the input signal, in order to select those frequencies that have to be analyzed with computationally expensive methods.
6-1
norma lized frequency (bin)
Fig. 5. Spectrum of a input signal.
n ormal ized frequ e ncy (b i n)
sinusoid
The proposed method based on a simple and efficient FFT scheme, it is advantageous in the applications in which the processing complexity and the short data record length constitute a critical issue, such as in real-time harmonic analysis.
1005
REFERENCES Griffiths,L.J. (1986). High-Resolution Spectral Estimation: Rethinking the Fourier Transform. Proc. ICASSP, 169-172. Greenhall, C.A. (1990). Orthogonal Sets of Data Windows Constructed from Trigonometric Polynomials. IEEE Trans. ASSP, 38 no.5, 870-872. Harris, F.J.(1978}. On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform. Proc. IEEE, 66 no.l, 51-83. Nuttall, A.H. (1981). Some Windows With Very Good Sidelobe Behavior. IEEE Trans. ASSP 29, 84-91. Offelli, C., Petri, D. (1990). Interpolation Techniques for Real-Time Mul tifrequency Waveform Analysis. IEEE Trans. Instrum. Heas., IH-39 no.l, 3106-3111. Slepian, D. (1978). Prolate Spheroidal Wave Functions, Fourier Analysis and uncertainty V: the Discrete Case. Bell Syst. Tech. Journ. 57 no.5, 1371-1430. Thompson, D.J. (1982). Spectrum Estimation and Harmonic Analysis. Proc. IEEE, 70 no. 9, 1055-1096.
1006