US 6,351,729 B1Grant
Multiple-window method for obtaining improved spectrograms of signals
Issue Date:2002-02-26
•3 Claims
•2 Drawing Sheets
Abstract
There is disclosed a method for processing a time-varying signal to produce a high-resolution spectrogram that represents power as a function of both frequency and time. Data blocks of a time series, which represents of a sampled signal, are subjected to processing which results in a sequence of frequency-dependent functions referred to as eigencoefficients. Each eigencoefficient represents signal information projected onto a local frequency domain using a respective one of K Slepian sequences or Slepian functions. The spectrogram is derived from time- and frequency-dependent expansions formed from the eigencoefficients.
Metadata
Assignee
- Lucent Technologies Inc.
Inventor
- David James Thomson
Application Information
Application Number:US 09/352,417
Filing Date:1999-07-12
Priority Date:1999-07-12
Art Unit:7
Classifications
IPC:
G10L 1104H04B 1500
Field of Search:
704226704233704203704205704206704207704220381 9437386125324 7622324 7624367 62702 76
Patent Drawings (2 sheets)
Description
Field of the Invention
[0002] The invention relates to methods for the spectral analysis of time-sampled signals. More particularly, the invention relates to methods for producing spectrograms of human speech or other time-varying signals.
Art Background
[0003] It is useful, in many fields of technology, to determine the changing frequency content of time-dependent signals. For example, the spectral analysis of speech is useful both for automatic speech recognition and for speech coding. As a further example, the spectral analysis of marine sounds is useful for acoustically aided undersea navigation.
[0004] When an acoustic signal, or other signal of interest, is sampled at discrete intervals, a time series is produced. A time series is said to be stationary if its statistical properties are invariant under displacements of the series in time. Although few of the signals of interest are truly stationary, many change slowly enough that, for purposes of spectral analysis, they can be treated as locally stationary over a limited time interval.
[0005] The spectral analysis of stationary time series has been a subject of research for one hundred years. The earliest attempts to obtain a representation, or periodogram, of the power spectral density of the time series x(0), x(1), . . . , x(n), . . . , x(N−1) involved summing N terms of the form x(n)×einω and then taking the squared magnitude of the result. (The symbol ω represents frequency in radians per second. The symbol 𝑓, used below, represents frequency in cycles per second. Thus, ω=2π𝑓.) This operation was performed for each of N/2+1 discrete frequencies 𝑓. This was unsatisfactory for several reasons. One reason is that the result is not statistically consistent. That is, the variance of the resulting periodogram does not decrease as the sample size N is increased. A second reason is that the result can be severely biased by truncation effects, leading to inaccurate representation of processes having continuous spectra.
[0006] An improved spectrum estimate (it is an estimate because it is derived from a finite sample of the original signal) is obtained from the following method, which is conveniently described in two steps:
[0007] First, form the spectrum estimate {tilde over (S)}D(ω) using a data window D0, D1, . . . , Dn, . . . , DN−1to taper the sampled data sequence, according to:
[0008] The primary purpose of the data window is to control bias. That is, by tapering the sampled sequence, it is possible to mitigate the tendency of the frequency components where the power is highest to dominate the spectrum estimate.
[0009] Then, smooth the estimate {tilde over (S)}D(ω) by convolving it with a spectral window G(ω) to form the smoothed spectrum estimate {tilde over (S)}(ω) according to {tilde over (S)}(ω)={tilde over (S)}D(ω)*G(ω),
[0010] where * represents the convolution operation. The primary purpose of the spectral window is to make the spectrum estimate consistent. The spectral window is generally pulse-shaped in frequency space, and the width of this pulse is approximately the bandwidth of the spectrum estimate. Increasing the bandwidth decreases the variance of the resulting estimate, but it also reduces the frequency resolution of the estimate.
[0011] Although useful, the smoothed spectrum estimate {tilde over (S)}(ω) as described above has several drawbacks. The smoothing operation may obscure the presence of spectral lines. Moreover, the data window tends to give different weights to equally valid data points. The data window also tends to reduce statistical efficiency. That is, the amount of data needed to obtain a reliable estimate may exceed the theoretical ideal by a factor of two or more.
[0012] Recently, a new spectrum estimate having improved properties was proposed. This estimate is described, e.g., in D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proc. IEEE70 (September 1982) 1055-1096 (hereafter, “Thomson (1982)”). This estimate is computed using a sequence of window functions referred to as Slepian functions when expressed as functions of frequency, and as Slepian sequences when expressed as sequences in the time domain. Slepian functions are related to Slepian sequences through the Fourier transform. Because multiple window functions are used, such an estimate is referred to as a multitaper spectrum estimate, or occasionally as a multiple-window spectrum estimate.
[0013] The properties of Slepian functions and Slepian sequences are described in Thomson (1982), cited above, and in D. Slepian, “Prolate Spheroidal Wave Functions, Fourier Analysis, and Uncertainty—V: The Discrete Case,” Bell System Tech. J.57 (1978) 1371-1430, hereafter referred to as Slepian (1978). Briefly, the Slepian sequences depend parametrically on the size N of the data sample and on the chosen bandwidth W. (From practical considerations, the bandwidth is generally chosen to lie between 1/N and 20/N, and at least as a starting value it is typically about 5/N.) It should be noted that throughout this discussion, the well-known convention is used wherein all frequencies are normalized such that the Nyquist frequency equals 0.5.
[0014] Given values for these parameters, each Slepian sequence v(k)(N,W) is a k'th solution to a matrix eigenvalue equation Mv=λkv, where the element in the n'th row and m'th column of the matrix is given by:
[0015] n=1, 2, . . . , N, m=1, 2, . . . , N.
[0016] If the eigenvalues λkof this equation are arranged in descending order, approximately the first K of them are very close to (but less than) unity. K is the greatest integer less than or equal to 2NW. At least for moderate values of N, the solutions are readily computed using standard techniques. (For such purpose, it is advantageous to use an alternative representation of these sequences which uses a matrix in tridiagonal form. For further information, see Slepian (1978), which is hereby incorporated by reference.)
[0017] The Slepian functions Uk(N,W;𝑓) are computed from corresponding Slepian sequences through the formula
[0018] where ε is 1 when k is even, and i when k is odd.
[0019] Of any function which is the Fourier transform of an index limited sequence, the k=0 Slepian function has the greatest fractional energy concentration within the frequency range between −W and W. More generally, the k'th eigenvalue λkexpresses the fraction of energy retained within this frequency range by the corresponding Slepian function. As noted, this fraction is very close to unity for the first K Slepian functions.
[0020] The spectrum estimate of Thomson (1982) is computed from K eigencoefficients y0(𝑓), Y1(𝑓) , . . . , yK−1(𝑓), wherein the k'th such eigencoefficient is computed through the formula,
[0021] At a given frequency 𝑓=𝑓0, the spectrum estimate, denoted {overscore (S)}(𝑓), is band limited to a frequency range of ±W about 𝑓0. The spectrum estimate is computed from the eigencoefficients according to,
[0022] It will be appreciated that each term in this summation is individually a spectrum estimate of the usual kind, as represented, e.g., by Equation (1), in which a respective Slepian sequence is the data window. In fact, the k=0 term is the optimal spectrum estimate of that kind, but even so, it must be smoothed in order to make it statistically consistent. Smoothing, however, tends to increase the effective bandwidth to several times W, and it concomitantly increases the bias of the estimate. On the other hand, when the rest of the eigencoefficients are included (up to the k=K−1 term), consistency and good variance efficiency are achieved without decreasing the spectral resolution.
[0023] Multiple window spectrum estimates are discussed further in D. J. Thomson, “Time Series Analysis of Holocene Climate Data,” Phil. Trans. R. Soc. Lond. A330 (1990) 601-616 (hereafter, “Thomson (1990)”). That article introduces a slightly different definition of the Slepian function, which uses a more common definition of the Fourier transform than the one used, e.g., in Slepian (1978). The Slepian function Vk(𝑓) of Thomson (1990) may be computed by Fourier transforming the corresponding Slepian sequence according to:
[0024] This form of the Slepian function is related to Uk(N,W;𝑓) by the expression:
[0025] The same article also introduces an alternate form xk(𝑓) for the eigencoefficients, given by
[0026] The same article also describes a multiple-window spectrum estimate {overscore (S)}(𝑓) computed by summing the squared magnitudes of the eigencoefficients xk(𝑓), each weighted by an appropriately chosen weight coefficient wk:
[0027] Thomson (1990) also describes a procedure for subdividing the data sequence into overlapping blocks, the base time of each block advanced by some offset from the base time of the preceding block, and computing the multiple-window spectrum estimate on each block.
[0028] It should be noted that each of the preceding spectrum estimates implicitly assumes stationarity. That is, each assumes that {overscore (S)}(𝑓) does not involve time, except for the implicit time dependence that comes from defining the sample on the discretized time block spanning the interval (0, N−1). On the other hand, spectrograms dealing explicitly with nonstationary processes have been used for many years. An early paper describing such techniques is W. Koenig et al., “The Sound Spectrograph,” J. Acoust. Soc. Am.18:19 (1946). In essence, these techniques involve estimates of the kind expressed by Equation (1), above, with the further property that the sample is stepped along in time. Thus, such an estimate might be wiritten as
[0029] where b now represents the base time, that is, the time (measured from a fixed origin) at the beginning of a given sample block, and n represents relative (discrete) time within the block. Thomson (1990) updated this idea by replacing {tilde over (S)}D(𝑓) with {overscore (S)}(𝑓) as in Equation (8), above.
Summary of the Invention
[0030] Significantly, the bandwidth-limited signal in the frequency band (𝑓−W,𝑓+W) can be expanded in the time block [0, N−1] as
[0031] where xk(𝑓) is defined as in Equation (7), above. This observation is made, e.g., in D. J. Thomson, “Multi-Window Bispectrum Estimates,” Proc. Workshop on Higher-Order Spectral Analysis, Vail, Colo. (Jun. 28-30, 1989). However, it has not been appreciated, until now, that such an expansion may be useful for formulating an improved spectrum estimate.
[0032] I have found an improved spectrum estimate that is based on the expansion described by Equation (10), above. Because this spectrum estimate depends explicitly on both time and frequency, I refer to it as a spectrogram. The time resolution of this spectrogram is approximately ½W. Because in typical applications the product 2NW is equal to the number K of Slepian sequences, an alternately formulated estimate for this bandwidth is N/K. By contrast, the time resolution of conventional spectrograms is typically roughly equal to the block size, N. Thus, my improved spectrogram is a high-resolution spectrogram.
[0033] In a broad aspect, my invention involves a method for processing a time-varying signal to produce a spectrogram. The method includes sampling the signal at intervals, thereby to produce a time series x(n), wherein x represents sampled signal values and n represents discretized time. The method further includes obtaining plural blocks of data x0, x1, . . . , xN−1from the time series, wherein each block contains signal values x(n) taken at an integer number N of successive sampling intervals.
[0034] The method further includes calculating an integer number K of eigencoefficients xk(𝑓) on each said block, wherein each said eigencoefficient is dependent on frequency 𝑓 and has a respective index k, k=0, 1, . . . , K−1. The method further includes, for each said block, forming a time- and frequency-dependent expansion X(t,𝑓) from the eigencoefficients, wherein t represents time.
[0035] The method further includes taking a squared magnitude of the expansion, and outputting a spectrogram derived at least in part from the resulting squared magnitude. Significantly, each eigencoefficient represents signal information projected onto a local frequency domain using a respective one of K Slepian sequences or Slepian functions. Moreover, each expansion X(t,𝑓) is a sum of terms, each term containing the product of an eigencoefficient and a corresponding Slepian sequence.
Brief Description of the Drawings
[0036] FIG. 1 is a schematic diagram illustrating a procedure or apparatus for computing an eigencoefficient from a block of sampled data, using Slepian sequences, in accordance with Equation (7).
[0037] FIG. 2 is a schematic diagram illustrating a procedure or apparatus for computing a spectrogram in accordance with aspects of the present invention as represented by Equation (11).
[0038] FIG. 3 is a schematic representation of a process of obtaining spectral data from overlapping blocks of sampled data for the purpose of averaging, according to the invention in one embodiment.
Detailed Description
[0039] In one simple form, the improved spectrogram is an expression F(t,𝑓) for power as a function of time and frequency, related to X(t,𝑓) by
[0040] FIG. 1 shows a procedure, in accordance with Equation (7), for obtaining eigencoefficients xk(𝑓). Data block 10 is a sequence of N signal values, sampled at discrete times and digitized. The signal values are provided by any appropriate devices for sensing and conditioning of signals, such as microphones and associated electronic circuitry. Each of blocks 20.1-20.N represents a weighted complex sinusoid in frequency space. For each value of the index k, each of the weights in blocks 20.1-20.N is one scalar term from the k'th Slepian sequence. As shown, each sampled signal value is multiplied by a corresponding weighted sinusoid, and the results are summed. Through the frequency dependence of the complex sinusoids, each of the resulting eigencoefficients is a complex-valued function of frequency.
[0041] It should be noted that the raw eigencoefficients as given by Equation (7) tend to exhibit exterior bias. That is, the Slepian sequences are not strictly band-limited; instead, each has a certain energy fraction that lies outside of the bandwidth W. Uncorrected, this out-of-band energy fraction contributes bias, which can be particularly severe for the higher-order eigencoefficients, that is, for those whose index k is close to K. Accordingly, one way to suppress exterior bias is to limit k to values no greater than, e.g., K−2 or K−4.
[0042] Another way to suppress bias is to use the adaptive weighting procedure described in Thomson (1982). According to that process, a weight coefficient is obtained for each eigencoefficient xk(𝑓). Each of these weight coefficients is a function of frequency. In Equation (11), each eigencoefficient is modified by multiplying it by its respective weight coefficient. The adaptive weighting procedure, which is described at pages 1065-1066 of Thomson (1982), obtains optimized weight coefficients by minimizing an error function which measures bias in pertinent spectral estimates.
[0043] Yet another, and currently preferred, method for suppressing bias is a procedure that I refer to as coherent sidelobe subtraction. This procedure also obtains weight coefficients for the eigencoefficients.
[0044] Let X(𝑓) be the finite Fourier transform of the data. Then, very briefly, the coherent sidelobe subtraction procedure begins with the following estimate of dX(𝑓⊕ξ), where the special symbol ⊕ indicates that the absolute value of ξ must be less than W:
[0045] Here, each {circumflex over (x)}k(1)is an estimate of an eigencoefficient. Next, using weighted, overlapped estimates of dZ, a global estimate of dZ is formed, much in the manner of local regression smoothing. Then, using an exterior convolution, the coherent bias on the various {circumflex over (x)}k(1)is estimated and subtracted. Further details are provided in Appendix I attached hereto.
[0046] FIG. 2 shows the assembly of the raw or weighted eigencoefficients into the spectrogram F(t,𝑓). Each of eigencoefficients 30.1-30.K is multiplied by a corresponding Slepian sequence. This multiplication is carried out such that the k'th eigencoefficient is multiplied by the k'th Slepian sequence. Significantly, each eigencoefficient is a function of (continuous) frequency, and each Slepian sequence is a function of (discrete) time. Thus, each resulting product is a function of both frequency and time. The products are summed to form X(t,𝑓) in accordance with Equation (10). The figure shows the formation of F(t,𝑓) by multiplying X(t,𝑓) by its complex conjugate and normalizing by 1/K . The signal processing of FIGS. 1 and 2 is readily carried out by a digital computer or digital signal processor acting under the control of an appropriate hardware, software, or firmware program.
[0047] In many cases, it will be most useful to apply the high-resolution spectrogram to data that are sampled in overlapping blocks. Such blocks are conveniently described in terms of the base time b, the relative time t within a frame (which may be thought of as an offset from the base time of the frame), and the absolute time t0, which at a given position within a given frame is the sum of the corresponding base time and offset: t0=b+t. In these terms, an expression for eigencoefficients yk(b,𝑓) in which the base position is made explicit is given by:
[0048] A corresponding spectrogram F(b⊕t,𝑓), in which the symbol ⊕ indicates that the offset t may be included in the sum only if it lies in the interval [0, N−1], is given by:
[0049] It should be noted in this regard that because the expansion of Equation (10), above, extrapolates the signal to times lying beyond the interval [0, N−1], the above restriction on the sum in the time argument is merely advisable, but not strictly necessary.
[0050] At the edges of blocks, it is possible for the spectrogram to exhibit error related to the well-known Gibbs phenomenon. This is advantageously mitigated through an averaging procedure. For example, the spectrogram is readily averaged over two or more overlapping blocks. Where the blocks overlap, the constituent values that contribute to the average at each point in time are taken at positions in their respective blocks for which the corresponding base time and offset have a common sum; i.e., for computing an average at t0, the constituent values are taken at respective positions for which b+t=t0.
[0051] Those skilled in the art will appreciate that such an average over overlapping blocks is advantageously made a weighted average. Exemplary weighting procedures are described in the attached Appendix II.
[0052] Significantly, the spectrogram of Eq. (14) can be extended to include many overlapping data sections, so high-resolution spectrograms of long data sets can be formed by averaging.
[0053] FIG. 3 illustrates an averaging process for overlapping data blocks. Each of sheets 50.1-50.3 represents a spectrogram obtained from a respective data block. The first of these blocks has a base time of 0, the second a base time of b1>0, and the third a base time of b2>b1. Sections A-A′, B-B′, and C-C′ represent frequency spectra taken from sheets 50.1, 50.2, and 50.3, respectively, at values of the time, measured within the respective blocks, that all correspond to the same absolute time t0. These spectra are readily averaged, as discussed above, to provide an average spectrum for each given value of the absolute time.
[0054] Appendix I: Coherent Sidelobe Subtraction
[0055] Begin with Equation (12). Note that for any frequency 𝑓0there is a range of frequencies (𝑓0−W, 𝑓0+W) giving an estimate of d{circumflex over (X)}(1)(f0), specifically
[0056] nominally independent of the free parameter ξ. Here {circumflex over (x)}k(p)(𝑓) is the estimate of xk(𝑓) at the pthinteration.
[0057] We use a weighted sum of the free-parameter expansions to form an estimate of dX
[0058] where the weighting function Q may reflect nothing more than that the convergence of the orthogonal expansions is generally poorer near the ends of the domain than in the center or, in regions where the spectrum is changing rapidly, that some expansions are less reliable than others.
[0059] Next, estimate the exterior bias of xk(𝑓) using the convolution over the exterior domain
[0060] and subtract it from the raw eigencoefficients to form an improved estimate
[0061] {circumflex over (x)}k(p+1)(𝑓)=yk(𝑓)−{circumflex over (b)}k(p+1)(𝑓) . (18)
[0062] The integral in Equation (17) is taken between the limits −1/2 to 1/2, but excluding the range −W to W.
[0063] Appendix II: Weighting Procedures for Averages Over Overlapping Blocks
[0064] One possible approach is to use a scaled version of the Epanechnikov kernel, which is known to be optimum in certain pertinent problems. The Epanechnikov kernel is described, e.g., in J. Fan and I. Gijbels, Local Polynomial Modelling and its Applications, Chapman and Hall, London, 1996. Very briefly, the Epanechnikov kernel K0(t) is given by:
[0065] Thus, one appropriate weighted average {overscore (F)}E(t0,𝑓) is given by:
[0066] A second possibility is to weight by Fisher information as well. An estimate Î(b,𝑓) of Fisher information is given by:
[0067] Using this estimate, an adaptively weighted average {overscore (F)}A(t0,𝑓) can be taken according to:
[0068] Here, as well as in {overscore (F)}E(t0,𝑓), above, the summation represented by
[0069] can be replaced by a sum at the Nyquist rate
[0070] This would give, for example:
Claims
The invention claimed is:
1. A method for processing a time-varying signal to produce a spectrogram, comprising:
a) sampling the signal at intervals, thereby to produce a time series x(t), wherein x represents sampled signal values and t represents discretized time;
b) obtaining plural blocks of data x0,x1, . . . ,xN−1from the time series, wherein each block contains signal values x(t) taken at an integer number N of successive sampling intervals;
c) calculating an integer number K of eigencoefficients xk(𝑓) on each said block, wherein each said eigencoefficient is dependent on frequency 𝑓 and has a respective index k, k=0, 1, . . . , K−1;
d) for each said block, forming a time- and frequency-dependent expansion X(t,f) from the eigencoefficients;
e) taking a squared magnitude of the expansion; and
f) outputting a spectrogram derived at least in part from the result of step (e), wherein:
I) each eigencoefficient represents signal information projected onto a local frequency domain using a respective one of K Slepian sequences or Slepian functions; and
II) each expansion X(t,𝑓) is a sum of terms, each term containing the product of an eigencoefficient and a corresponding Slepian sequence.
2. The method of claim 1, wherein the signal information projected in each eigencoefficient is sampled at offsets 0, 1, . . . , N−1 from a base position b within the time series.
3. The method of claim 2, wherein:
each block overlaps at least one other block in an overlap region;
in each overlap region, the spectrogram is averaged over overlapping blocks; and
said averaging is carried out over respective combinations of base position and offset that have a common sum.
Patent Citations (10)
| Patent | Date | Inventor | Cited By |
|---|---|---|---|
| US3573446(A) | 1971-04-01 | Bergland | Examiner |
| US4802225(A) | 1989-01-01 | Patterson | Examiner |
| US4983906(A) | 1991-01-01 | Hiller | Examiner |
| US5179626(A) | 1993-01-01 | Thomson | Examiner |
| US5768392(A) | 1998-06-01 | Graupe | Examiner |
| US6044340(A) | 2000-03-01 | Van Hamme | Examiner |
| US6070137(A) | 2000-05-01 | Bloebaum et al. | Examiner |
| US6122610(A) | 2000-09-01 | Isabelle | Examiner |
| US6249762(B1) | 2001-06-01 | Kirsteins et al. | Examiner |
| US6263306(B1) | 2001-06-01 | Fee et al. | Examiner |
Non-Patent Literature (9)
- D. J. Thomson, “Spectrum Estimation and Harmonic Analysis,” Proc. IEEE 70, pp. 1055-1096 (1982).Applicant
- D. Slepian, “Prolate Spheroidal Wave Functions, Fourier Analysis, and Uncertainty-V: The Discrete Case,” Bell System Tech. J. 57, pp. 1371-1430 (1978).Applicant
- D. J. Thomson, “Time Series Analysis of Holocene Climate Data,” Phil Trans. R. Soc. Lond. A 330 pp. 601-616 (1990).Applicant
- W. Koenig et al., “The Sound Spectrograph,” J. Acoust. Soc. Am. 18:19 (1946).Applicant
- D. J. Thomson, “Multi-Window Bispectrum Estimates,” Proc. Workshop on Higher-Order Spectral Analysis, pp. 19-23, (1989).Applicant
- D. J. Thomson, “Multiple-Window Spectrum Estimates for Non-Stationary Data”, Proc. Workshop on Statistical Signal and Array, pp. 344-347 (1998).Applicant
- J. Fan et al. Local Polynominal Modelling and its Applications Chapman and Hall, London, 1996.Applicant
- D. J. Thomson “An Overview of Multiple-Window and Quadratic-Inverse Spectrum Estimation Methods”, IEEE Proc. on Acoustics, Speech and Signal Processing, pp. VI-185-VI-194, (1994).Applicant
- L.J. Lanzerotti, et al. “Statistical Properties Of Shock-Accelerated Ions In The Outer Heliosphere” The Astrophysical Journal, 380, pp. L93-L96, (1991).Applicant