Method and apparatus for estimating the fundamental frequency of a signal

[0002] The present invention relates, in general, to communication systems and, more particularly, to coding information signals in such communication systems.

[0003] Digital speech compression systems typically require estimation of the fundamental frequency of an input signal. The fundamental frequency 𝑓₀is usually estimated in terms of the pitch period τ₀(otherwise known as “lag”). The two are related by the expression $\begin{array}{cc}{\tau }_0=\frac{f_s}{f_0},& \left(1\right)\end{array}$

[0004] where the sampling frequency 𝑓_sis commonly 8000 Hz for telephone grade applications.

[0005] Since a speech signal is generally non-stationary, it is partitioned into finite length vectors called frames (e.g., 10 to 40 ms), each of which are presumed to be quasi-stationary. The parameters describing the speech signal are then updated at the associated frame length intervals. The original Code Excited Linear Prediction (CELP) algorithm further updates the pitch period (using what is called Long Term Prediction, or LTP) information on shorter subframe intervals, thus allowing smoother transitions from frame to frame. It was also noted that although τ₀could be estimated using open-loop methods, far better performance was achieved using the closed-loop approach. Closed-loop methods involve an exhaustive search of all possible values of τ₀(typically integer values from 20 to 147) on a subframe basis, and choosing the value that satisfies some minimum error criterion.

[0006] An enhancement to this method involves allowing τ₀to take on fractional values. An example of a practical implementation of this method can be found in the GSM half rate speech coder, and is shown in FIG. 1. Here, lags within the range of 21 to 22⅔ are allowed ⅓ sample resolution, lags within the range of 23 to 34⅚ are allowed ⅙ sample resolution, and so on. In order to keep the search complexity low, a combination of open-loop and closed loop methods is used. The open-loop method involves generating an integer lag candidate list using an autocorrelation peak picking algorithm. The closed-loop method then searches the allowable lags in the neighborhood of the integer lag candidates for the optimal fractional lag value. Furthermore, the lags for subframes 2, 3, and 4 are coded based on the difference from the previous subframe. This allows the lag information to be coded using fewer bits since there is a high intraframe correlation of the lag parameter. Even so, the GSM HR codec uses a total of 8+(3×4)=20 bits every 20 ms (1.0 kbps) to convey the pitch period information.

[0007] In an effort to reduce the bit rate of the pitch period information, an interpolation strategy was developed that allows the pitch information to be coded only once per frame (using only 7 bits→350 bps), rather than with the usual subframe resolution. This technique is known as relaxed CELP (or RCELP), and is the basis for the recently adopted enhanced variable rate codec (EVRC) standard for Code Division Multiple Access (CDMA) wireless telephone systems. The basic principle is as follows.

[0008] The pitch period is estimated for the analysis window centered at the end of the current frame. The lag (delay) contour is then generated, which consists of a linear interpolation of the past frame's lag to the current frame's lag. The linear prediction (LP) residual signal is then modified by means of sophisticated polyphase filtering and shifting techniques, which are designed such that the ⅛ sample interpolation boundaries are not crossed during perceptually critical instances in the waveform. The primary reason for this residual modification process is to account for errors introduced by the open-loop integer lag estimation process. For example, if the integer lag is estimated to be 32 samples, when in fact the true lag is 32.5 samples, the residual waveform can be in conflict with the estimated lag by as many as 2.5 samples in a single 160 sample frame. This can severely degrade the performance of the LTP. The RCELP algorithm accounts for this by shifting the residual waveform during perceptually insignificant instances in the residual waveform (i.e., low energy) to match the delay contour. In the event that there are no such opportunities for shifting, the shift count is accumulated and reserved for use during the next frame. By modifying the residual waveform to match the estimated delay contour, the effectiveness of the LTP is preserved, and the coding gain is maintained. In addition, the associated perceptual degradations due to the residual modification are claimed to be insignificant.

[0009] But, while this last claim may be true for medium bit rate coders such as the EVRC full rate mode (i.e., 8.5 kbps), it is less apparent for the EVRC half rate mode, which operates at 4.0 kbps. This is because of the relative ability of the fixed codebooks to model the associated inverse error signal. That is, if coding distortions are introduced by inefficiencies in the LTP, and those distortions can be effectively modeled by the fixed codebook, then the net effect is that the distortion will be canceled. So, while the EVRC full rate mode allocates 120 of 170 its per frame for fixed codebook gain and shape, the half rate mode can afford only 42 of 80 bits per frame for the same. This results in a disproportionate performance degradation due, in part, to the fixed codebook's inability to model the coding distortion introduced by the LTP.

[0010] Therefore, there is a need for an improved method of open-loop pitch estimation that provides subsample resolution.

[0011] FIG. 1 generally depicts fractional lag values for a GSM half-rate speech coder.

[0012] FIG. 2 generally depicts a speech compression system employing open-loop lag estimation as is known in the prior art.

[0013] FIG. 3 generally depicts a open-loop lag estimation system in accordance with the invention.

[0014] FIG. 4 generally shows the structure of an ith adaptive filter element within the filter bank of FIG. 3.

[0015] FIG. 5 generally depicts the process of variable length sequencing, variable offset, and subsequent windowing in accordance with the invention.

[0016] FIG. 6 generally depicts an example of an m=7 bit trained scalar quantization table in accordance with the invention.

[0017] FIG. 7 generally depicts a comparison of voiced speech lag estimation between a prior art method and lag estimation in accordance with the invention.

[0018] FIG. 8 generally depicts a comparison of average absolute accumulated shift between a prior art method and lag estimation in accordance with the invention.

[0019] Stated generally, a method and apparatus for improved pitch period estimation in a compression system is disclosed. The system uses original estimates of integer lag and open-loop prediction gain as input to an adaptive filter parameter initialization block which supplies inputs to a plurality of adaptive filter elements. Adaptive filter elements provide information regarding the harmonics of the residual signal to an adaptive filter parameter analysis block. Adaptive filter parameter analysis block estimates the fundamental frequency of the residual signal based on the analysis of the harmonics and outputs a pitch period for eventual use in a delay contour computation.

[0020] Stated more specifically, a method for estimating a fundamental frequency of a signal includes the steps of analyzing harmonics of the signal and estimating the fundamental frequency of the signal based on the analysis of the harmonics. In the preferred embodiment, the step of analyzing harmonics of the signal further comprises phase locking to the harmonics of the signal and the step of estimating further comprises quantizing the fundamental frequency of the signal and converting the quantized fundamental frequency of the signal to the lag domain for use as a pitch period τ. In this embodiment, the signal further comprises a residual of a speech coded signal.

[0021] FIG. 2 generally depicts a speech compression system 200 employing open-loop lag estimation as is known in the prior art. As shown in FIG. 2, the input speech signal s(n) is processed by an linear prediction (LP) analysis filter 202 which flattens the short-term spectral envelope of input speech signal s(n). The output of the LP analysis filter is designated as the LP residual ε(n). The LP residual signal ε(n) is then used by the open-loop lag estimator 204 as a basis for estimating the delay contour τ_c(n) and the open-loop pitch prediction gain β_ol. The RCELP residual modification process 206 uses this information to map the LP residual to the delay contour, as described above. The modified residual signal is then passed through a weighted synthesis filter 207 before being processed by the long term predictor 208 and eventually by the fixed codebook 210, which characterizes the synthesizer excitation sequence. At the decoder side, the fixed codebook index/gain is input to an excitation generator 212 which outputs an excitation sequence. The excitation sequence is passed through long and short term synthesis filters 214 and 216 to produce the reconstructed speech output.

[0022] In the preferred embodiment of the present invention, a more accurate estimate of the RCELP delay contour is produced which results in a more accurate mapping of the delay contour to the LP residual signal ε(n). FIG. 3 generally depicts a open-loop lag estimation system 300 in accordance with the invention. As shown in FIG. 3, the LP residual signal ε(n) is used as the input to the autocorrelation analysis block 302 which uses the prior art method shown in section 4.2.3 of IS-127 as “pre-optimized” values for the integer lag τ₀and open-loop prediction gain β_ol. These values are then used to calculate the initial parameters for the adaptive harmonic filter bank. The filter bank is used to estimate the residual signal's fundamental frequency 𝑓₀by “phase locking” to the LP residual harmonic frequencies. The fundamental frequency is then quantized and converted to the lag domain for use as the optimal pitch period τ in accordance with the invention.

[0023] To explain open-loop lag estimation in accordance with the invention, it is given [6] that the complex conjugate roots of the recursive digital filter are of the form: $\begin{array}{cc}H\left(z\right)=\frac{\gamma }{\left(1-{\mathrm{re}}^{{\mathrm{j\omega }}_0}{z}^{-1}\right)\left(1-{\mathrm{re}}^{{\mathrm{j\omega }}_0}{z}^{-1}\right)}=\frac{\gamma }{1-2r\cos \left({\omega }_0\right)z^{-1}+r^2{z}^{-2}}& \left(2\right)\end{array}$

[0024] which results in a bandpass frequency response with center frequency ω₀. As the value of r approaches 1, the filter response becomes more and more narrow, and the filter structure can be classified as a digital resonator. A value of r=1 results in oscillation, and values of r>1 result in instability. The value of γ that yields unity gain at the center frequency ω₀is given as: $\begin{array}{cc}\gamma =\left(1-r\right)\sqrt{1+r^2-2r\cos 2{\omega }_0}.& \left(3\right)\end{array}$

[0025] Furthermore, by modifying H(z) to include a unit delay in the numerator: $\begin{array}{cc}H\left(z\right)=\frac{\gamma z^{-1}}{1-2r\cos \left({\omega }_0\right)z^{-1}+r^2{z}^{-2}}& \left(4\right)\end{array}$

[0026] the phase response is such that there is a phase shift of −π/2(−90°) at the center frequency ω₀, and an asymptotical approach to 0 and −π on the low and high frequency sides of ω₀, respectively. This is an important point because for an input signal x(n), the output signal y(n) would be orthogonal to the input when evaluated at ω₀. Extending this idea to an adaptive filter context, let us consider the following time varying transfer function: $\begin{array}{cc}H\left(q,n\right)=\frac{\gamma q^{-1}}{1-a\left(n\right)q^{-1}-{\mathrm{bq}}^{-2}},& \left(5\right)\end{array}$

[0027] where q⁻is a unit delay operator and n represents the sampled time index. For a sinusoidal input, if a(n) were allowed to vary until the output of the filter were orthogonal to the input, the frequency of the input ω_incould be determined at time n by the relation:

[0028] ω_in=cos⁻¹(a(n)/2r).  (6)

[0029] This is the basic premise of the invention. The input signal ε(n) is used as an input to a filterbank of adaptive filter elements 306-308. Each of the adaptive filter elements 306-308 pre-filters the residual at a fixed harmonic frequency. FIG. 4 shows the corresponding structure of an ith adaptive filter element, for example adaptive filter element 306, contained within the filter bank of FIG. 3. The pre-filtered output is windowed 404 to prevent spurious signal conditions, and then used as input to an adaptive resonator 406 of the form given in Eq. (5). Once the data for the current frame has been processed, the adaptive filter coefficients for each of the harmonics is analyzed, and an estimate of the fundamental frequency is generated.

[0030] The LP residual signal ε(n) is first filtered by the zero-state harmonic pre-filter 402, given as:

[0031] p_i(n)=γ_1iε(n+j)+c_ip_i(n−1)+d_ip_i(n−2), 0≦n<L, 1≦i≦N,  (7)

[0032] where N is the number of harmonics to be analyzed, and the filter coefficients are given as: $\begin{array}{cc}\begin{array}{c}{c}_i=2r_1\cos {\omega }_i\\ d_i=-r_1^2\end{array}\},1\le i\le N.& \left(8\right)\end{array}$

[0033] In this expression, the pole radius is given as the constant r₁=0.989 and the harmonic frequency ω_iis given as: $\begin{array}{cc}{\omega }_i=\frac{2\pi }{{\tau }_0}\left(h+i-1\right),1\le i\le N,& \left(9\right)\end{array}$

[0034] where τ₀is the pre-optimized lag from autocorrelation analysis block 302, h is the harmonic number corresponding to the smallest frequency greater than the minimum frequency 𝑓_min: $\begin{array}{cc}h=\lceil \frac{{\tau }_0{f}_{\min }}{f_s}\rceil ,& \left(10\right)\end{array}$

[0035] and 𝑓f_sis the sampling frequency. In this embodiment, 𝑓_min=800 Hz,𝑓_s=8000 Hz, and N=4 harmonics. So, for example, if the lag τ₀were chosen to be 40, the harmonic filter bank would be configured to analyze the 4th through 7th harmonics of the fundamental frequency 𝑓₀=𝑓_s/τ₀=200 Hz, or more explicitly 800, 1000, 1200, and 1400 Hz.

[0036] Referring back to Eqs. (3) and (7), the filter gain variable γ_1ican then be calculated as:

[0037] $\begin{array}{cc}{\gamma }_{1i}=\left(1-r_1\right)\sqrt{1+r_1^2-2r_1\cos 2{\omega }_i},1\le i\le N.& \left(11\right)\end{array}$

[0038] In some applications, however, it may be more desirable (from a complexity standpoint) to hold this value constant across the range of allowable frequencies, however, care must be taken to assure that the adaptive filter gains are not biased as a result.

[0039] Also in Eq. (7), two symbols require definition. First, the sequence length L is defined such that at least three full pitch periods must be contained within the LP residual signal ε(n), up to a given maximum. This guarantees a meaningful input to the adaptive filters 306-308. The sequence length L can then be defined as: $\begin{array}{cc}L=\{\begin{array}{cc}3L_f/2,& {\tau }_0\le L_f/2\\ 3{\tau }_0,& L_f/2<{\tau }_0<2L_f/3,\\ 2L_f,& \mathrm{otherwise}\end{array}& \left(12\right)\end{array}$

[0040] where L_𝑓=160 is the frame length of this embodiment. The associated variable offset j into the input sequence ε(n) can then be expressed as:

[0041] $\begin{array}{cc}j=\{\begin{array}{cc}{L}_f/2,& {\tau }_0\le L_f/2\\ 2L_f-3{\tau }_0,& L_f/2<{\tau }_0<2L_f/3.\\ 0,& \mathrm{otherwise}\end{array}& \left(13\right)\end{array}$

[0042] The process of variable length sequencing, variable offset, and subsequent windowing can be more readily observed in FIG. 5. In this embodiment, the LPC residual sequence ε(n) consists of 320 samples of information, the first 80 of which represent the last half of the last frame of information. These are used to “prime” the state of the LP analysis filter 202 and are used here to extend the pitch analysis frame for low frequencies. The next 160 samples represent the current frame of information, and the following 80 samples represent the “look-ahead” to the next frame. Since RCELP attempts to interpolate the lag from frame-to-frame, it is desirable to estimate the lag corresponding to the point at the end of the current frame. To do this adequately, a “look-ahead” to the next frame thereby centering the analysis frame towards the end of the current frame must occur. In the event that the initial lag is greater than half a frame length, however, the analysis frame is elongated towards past samples. This is a conscious effort to not increase algorithmic delay by not increasing the length of the look-ahead. Using this method, the sequence length L is equal to 240 for lags ≦80. For lags >80 and ≦120, the analysis window is stretched back in time successively until the beginning of the look-back of the previous frame is reached.

[0043] Also implicit in FIG. 5 is the windowing 404 of the pre-filtered output sequence p_i(n), which can be expressed as:

[0044] x_i(n)=w(n)p_i(n), 0≦n<L,1≦i≦N,  (14)

[0045] where ω(n) is the window: $\begin{array}{cc}w\left(n\right)=w\left(L-1-n\right)=\{\begin{array}{cc}{\sin }^2\left(\pi \left(n+0.5\right)/L_f\right),& 0\le n<L_f/2\\ 1,& L_f/2\le n<L/2\end{array}.& \left(15\right)\end{array}$

[0046] The window ω(n) can be described as a smoothed trapezoid window. Other window types may be used with varying degrees of performance, however, keeping the window “tails” constant is a computational advantage since only L_𝑓/2 samples need be stored or calculated. Other window types that have dependence on L in the trigonometric function e.g., sin²(π(n+0.5)/L)) would require recalculation and/or different storage memory for each value of L.

[0047] The windowed pre-filter output x(n) is then used as input to a zero-state adaptive harmonic resonator 406 of the form:

[0048] y_i(n)=γ_2ix(n)+a_i(n)y_i(n−1)+b_iy_i(n−2), 0≦n<L,1≦i≦N,  (16)

[0049] where the fundamental difference between this and the pre-filter in Eq. (7) is that the filter coefficient corresponding to the filter resonant frequency a_i(n) in Eq. (16) is allowed to vary with time. The initial coefficient values are calculated as in Eq.(8) as: $\begin{array}{cc}\begin{array}{c}{a}_i\left(0\right)=2r_2\cos {\omega }_i\\ b_i=-r_2^2\end{array}\},1\le i\le N.& \left(17\right)\end{array}$

[0050] where the pole radius is r₂=0.989 and the resonant frequency ω_iis given in Eqs. (9) and (10). The filter gain term is also calculated as in Eq. (11) by: $\begin{array}{cc}{\gamma }_{2i}=\left(1-r_2\right)\sqrt{1+r_2^2-2r_2\cos 2{\omega }_i},1\le i\le N.& \left(18\right)\end{array}$

[0051] As the windowed pre-filter sequence is passed through the second filter 406, which in the preferred embodiment is an adaptive harmonic resonator 406, the resonant frequency coefficient for each harmonic is modified according to:

[0052] a_i(n+1)=a_i(n)+α_ix_i(n)y_i(n−1), 0≦n>L,1≦i≦N,  (19)

[0053] where a_iis the normalized adaptation gain given by: $\begin{array}{cc}{\alpha }_i=\frac{{\alpha }_1{\beta }_{\mathrm{ol}}}{k+\sum _{n=0}^{L-1}{x}_i^2\left(n\right)},1\le i\le N.& \left(20\right)\end{array}$

[0054] In this embodiment, the gain constant α₁=0.49, the reciprocal limiting constant k=150, and the open-loop prediction gain (obtained from the autocorrelation analysis block in FIG. 3) has the range 0≦β_ol≦1. These values assume the input speech to be normalized within the 16 bit range ±32767. Also in this embodiment, the gain term depends on the entire L samples to be passed through the pre-filter/window prior to the execution of the adaptive filter.

[0055] Once all L samples have been passed through the adaptive harmonic filter 406, the coefficients a_i(L) can be analyzed for deviation from the initial center frequencies. The inverse of Eqs. (8) and (9) can be applied to each of the harmonic elements of the filter bank to form the following estimation of the unquantized fundamental frequency 𝑓₀: $\begin{array}{cc}{f}_0=\frac{f_s}{2\pi }\sum _{i=1}^N\lambda \left(i\right)\frac{{\cos }^{-1}\left(a_i\left(L\right)/2r_2\right)}{h+i-1},& \left(21\right)\end{array}$

[0056] where λ(i) is a weighting function that appropriately weights the importance of each of the harmonic elements, such that the sum of all elements of λ(i) equal unity. In this embodiment, λ(i) is equal to a linear average, or λ(i)=1 /N, 1≦i≦N. λ(i) is highly dependent on the input data; other functions may or may not yield better performance.

[0057] The quantized value of the fundamental frequency 𝑓* is then found by choosing the value of the index k that minimizes the following:

[0058] min {(𝑓₀−𝑓_table(k))²}, 0≦k>2^m,  (22)

[0059] where f_table(k) are the allowable values of the quantized fundamental frequency and m is the number of bits allocated to the lag parameter. An example of an m=7 bit trained scalar quantization table is given in FIG. 6.

[0060] The quantized fundamental frequency 𝑓* can then be used to generate the corresponding delay contour τ_c(n) which is output from the lag contour computation block 312 of FIG. 3, for which the following example is given as: $\begin{array}{cc}{\tau }_c\left(n\right)=\left\{\begin{array}{cc}\tau \left(m-1\right)+\left(n\left(\tau \left(m-1\right)-\tau \left(m\right)\right)\right)/L_f,& \left|\tau \left(m-1\right)-\tau \left(m\right)\right|<16\\ \tau \left(m\right),& \mathrm{otherwise}\end{array}\right\},& \left(23\right)\\ 0\le n<L_f& \end{array}$

[0061] where τ=𝑓_s/𝑓*, m represents the current frame and m−1 is the previous frame. In this equation, the interpolated lag is considered valid only if the trajectory does not change too abruptly.

[0062] In the preferred embodiment, a database of over 80,000 frames of speech and music signals was used to generate the data shown in FIG. 6 based on the dataset's probability density function. While this data is statistically optimal over the various fundamental frequency values from the training set, it is interesting to note that the distribution more closely reflects the properties of the human auditory system. That is, psychoacoustics principles reveal that the critical bands of hearing are uniform in frequency below 500 Hz; in the prior art open-loop lag estimator shown in FIG. 2, the distribution of quantization range is uniform over pitch period, which is inversely proportional to frequency. Thus, for a given fundamental frequency range of say 66 to 400 Hz, the psychoacoustically optimal distribution for a 7 bit quantizer would consist of a uniform frequency distribution spaced at 2.6 Hz intervals. The table shown in FIG. 6 yields relatively constant spacing through about 250 Hz (τ=32), but then sharply increases to about 20 Hz at the end frequency of about 400 Hz (τ≈20). This decrease in resolution is due to the diminished probability of encountering very high frequency talkers. So, the present invention facilitates the combination of both statistically and psychoacoustically joint optimal datasets by allowing arbitrary quantization levels for the fundamental frequencies. This is not readily achievable in the prior art.

[0063] The support for objective improvement can be observed in FIG. 7, where the lag trajectory for a short passage of strongly voiced speech is shown. While the prior art shows distinct “staircase” effects during transitional (frames 52 to 58 ) and steady state (frames 37 to 45 ) passages, the lag estimation in accordance with the invention effectively smoothes out the rough edges associated with the integer lag boundaries.

[0064] As described above, improvements to the RCELP algorithm can be evaluated objectively by measuring the accumulated shift that results from the inability of the LP residual signal ε(n) to be appropriately mapped to the estimated delay contour. Since one purpose of the preferred embodiment in accordance with the invention is to more accurately estimate the RCELP delay contour, the efficiency of the RCELP algorithm is improved since lag estimation in accordance with the invention requires less error to be tolerated by lowering the accumulated average shift factor. The improvement can be observed in FIG. 8, which was generated from a 80,000+ frame database.

[0065] Additionally, the subjective performance improvement is highly audible. Testing shows consistent preference to lag estimation in accordance with the invention during blind A/B tests, and it is estimated that the inventive method and apparatus provides 0.1 to 0.2 Mean Opinion Score (MOS) points improvement in Absolute Category Rating (ACR) tests when used with the EVRC half-rate maximum mode of operation (4.0 kbps).

[0066] While the invention has been particularly shown and described with reference to a particular embodiment, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention. For example, one skilled in the art will recognize that lag estimation in accordance with the invention can additionally benefit other, more general algorithms/vocoders which require accurate open-loop estimation of the fundamental frequency of an input signal. Such algorithms/vocoders include, but are not limited to, harmonic vocoders, sinusoidal transform coders (STC), and homomorphic vocoders. In addition to cellular communication systems, other applications which may benefit include digital hearing aids, audio speech coders, voice mail systems, etc. The corresponding structures, materials, acts and equivalents of all means or step plus function elements in the claims below are intended to include any structure, material, or acts for performing the functions in combination with other claimed elements as specifically claimed.