# Patent

## Title
Noise variance estimator for speech enhancement

## Abstract

A speech enhancement method operative for devices having limited available memory is described. The method is appropriate for very noisy environments and is capable of estimating the relative strengths of speech and noise components during both the presence as well as the absence of speech.

## Background
TECHNICAL FIELD
The invention relates to audio signal processing. More particularly, it relates to speech enhancement and clarification in a noisy environment.
INCORPORATION BY REFERENCE
The following publications are hereby incorporated by reference, each in their entirety.

    [1] Y. Ephraim and D. Malah, “Speech enhancement using a minimum mean square error short time spectral amplitude estimator,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 32, pp. 1109-1121, Dec. 1984.
    [2] N. Virag, “Single channel speech enhancement based on masking properties of the human auditory system,” IEEE Tran. Speech and Audio Processing, vol. 7, pp. 126-137, Mar. 1999.
    [3] R. Martin, “Spectral subtraction based on minimum statistics,” in Proc. EUSIPCO, 1994, pp. 1182-1185.
    [4] P. J. Wolfe and S. J. Godsill, “Efficient alternatives to Ephraim and Malah suppression rule for audio signal enhancement,” EURASIP Journal on Applied Signal Processing, vol. 2003, Issue 10, Pages 1043-1051, 2003.
    [5] Y. Ephraim, H. Lev-Ari and W. J. J. Roberts, “A brief survey of Speech Enhancement,” The Electronic Handbook, CRC Press, Apr. 2005.


BACKGROUND ART
We live in a noisy world. Environmental noise is everywhere, arising from natural sources as well as human activities. During voice communication, environmental noises are transmitted simultaneously with the intended speech signal, adversely effecting the quality of a received signal. This problem is mitigated by speech enhancement techniques that remove such unwanted noise components, thereby producing a cleaner and more intelligible signal.
Most speech enhancement systems rely on various forms of an adaptive filtering operation. Such systems attenuate the time/frequency (T/F) regions of the noisy speech signal having low Signal-to-Noise-Ratios (SNR) while preserving those with high SNR. The essential components of speech are thus preserved while the noise component is greatly reduced. Usually, such a filtering operation is performed in the digital domain by a computational device such as a Digital Signal Processing (DSP) chip.
Subband domain processing is one of the preferred ways in which such adaptive filtering operation is implemented. Briefly, the unaltered speech signal in the time domain is transformed to various subbands by using a filterbank, such as the Discrete Fourier Transform (DFT). The signals within each subband are subsequently suppressed to a desirable amount according to known statistical properties of speech and noise. Finally, the noise suppressed signals in the subband domain are transformed to the time domain by using an inverse filterbank to produce an enhanced speech signal, the quality of which is highly dependent on the details of the suppression procedure.
An example of a prior art speech enhancer is shown in FIG. 1. The input is generated by digitizing an analog speech signal that contains both clean speech as well as noise. This unaltered audio signal y(n), where n=0,1, . . . , ∞ is the time index, is then sent to an analysis filterbank device or function (“Analysis Filterbank”) 2, producing multiple subbands signals, Yk(m), k=1,. . . , K, m=0,1, . . . , ∞, where k is the subband number, and m is the time index of each subband signal. The subband signals may have lower sampling rates compared with y(n) due to the down-sampling operation in Analysis Filterbank 2. The noise level of each subband is then estimated by using a noise variance estimator device or function (“Noise Variance Estimator”) 4 with the subband signal as input. The Noise Variance Estimator 4 of the present invention differs from those known in the prior art and is described below, in particular with respect to FIGS. 2a and 2b. Based on the estimated noise level, appropriate suppression gains gk are determined in a suppression rule device or function (“Suppression Rule”) 6, and applied to the subband signals as follows:

{tilde over (Y)}k(m)=gkYk(m), k=1, . . . , K.  (1)

Such application of the suppression gain to a subband signal is shown symbolically by a multiplier symbol 8. Finally, {tilde over (Y)}k(m) are sent to a synthesis filterbank device or function (“Synthesis Filterbank”) 10 to produce an enhanced speech signal {tilde over (y)}(n). For clarity in presentation, FIG. 1 shows the details of generating and applying a suppression gain to only one of multiple subband signals (k).

The appropriate amount of suppression for each subband is strongly correlated to its noise level. This, in turn, is determined by the variance of the noise signal, defined as the mean square value of the noise signal with respect to a zero-mean Gaussian probability distribution. Clearly, an accurate noise variance estimation is crucial to the performance of the system.
Normally, the noise variance is not available, a priori, and must be estimated from the unaltered audio signal. It is well-known that the variance of a “clean” noise signal can be estimated by performing a time-averaging operation on the square value of noise amplitudes over a large time block. However, because the unaltered audio signal contains both clean speech and noise, such a method is not directly applicable.
Many noise variance estimation strategies have been previously proposed to solve this problem. The simplest solution is to estimate the noise variance at the initialization stage of the speech enhancement system, when the speech signal is not present (reference [1]). This method, however, works well only when the noise signal as well as the noise variance is relatively stationary.
For an accurate treatment of non-stationary noise, more sophisticated methods have been proposed. For example, Voice Activity Detection (VAD) estimators make use of a standalone detector to determine the presence of a speech signal. The noise variance is only updated during the time when it is not (reference [2]). This method has two shortcomings. First, it is very difficult to have reliable VAD results when the audio signal is noisy, which in turn affects the reliability of the noise variance estimation result. Secondly, this method precludes the possibility to update the noise variance estimation when the speech signal is present. The latter concern leads to inefficiency because the noise variance estimation can still be reliably updated during times wherein the speech level is weak.
Another widely quoted solution to this problem is the minimum statistics method (reference [3]). In principle, the method keeps a record of the signal level of historical samples for each subband, and estimates the noise variance based on the minimum recorded value. The rationale behind this approach is that the speech signal is generally an on/off process that naturally has pauses. In addition, the signal level is usually much higher when the speech signal is present. Therefore, the minimum signal level from the algorithm is probably from a speech pause section if the record is sufficiently long in time, yielding a reliable estimated noise level. Nevertheless, the minimum statistics method has a high memory demand and is not applicable to devices with limited available memory.
DISCLOSURE OF THE INVENTION
According to a first aspect of the invention, speech components of an audio signal composed of speech and noise components are enhanced. An audio signal is transformed from the time domain to a plurality of subbands in the frequency domain. The subbands of the audio signal are subsequently processed. The processing includes adaptively reducing the gain of ones of the subbands in response to a control. The control is derived at least in part from an estimate of variance in noise components of the audio signal. The estimate is, in turn, derived from an average of previous estimates of the amplitude of noise components in the audio signal. Estimates of the amplitude of noise components in the audio signal having an estimation bias greater than a predetermined maximum amount of estimation bias are excluded from or underweighted in the average of previous estimates of the amplitude of noise components in the audio signal. Finally, the processed audio signal is transformed from the frequency domain to the time domain to provide an audio signal in which speech components are enhanced. This aspect of the invention may further include an estimation of the amplitude of noise components in the audio signal as a function of an estimate of variance in noise components of the audio signal, an estimate of variance in speech components of the audio signal, and the amplitude of the audio signal.
According to a further aspect of the invention, an estimate of variance in noise components of an audio signal composed of speech and noise components is derived. The estimate of variance in noise components of an audio signal is derived from an average of previous estimates of the amplitude of noise components in the audio signal. The estimates of the amplitude of noise components in the audio signal having an estimation bias greater than a predetermined maximum amount of estimation bias are excluded from or underweighted in the average of previous estimates of the amplitude of noise components in the audio signal. This aspect of the invention may further include an estimation of the amplitude of noise components in the audio signal as a function of an estimate of variance in noise components of the audio signal, an estimate of variance in speech components of the audio signal, and the amplitude of the audio signal.
According to either of the above aspects of the invention, estimates of the amplitude of noise components in the audio signal having values greater than a threshold in the average of previous estimates of the amplitude of noise components in the audio signal may be excluded or underweighted.
The above mentioned threshold may be a function of ψ(1+{circumflex over (ξ)}(m)){circumflex over (λ)}d(m), where {circumflex over (ξ)} is the estimated a priori signal-to-noise ratio, {circumflex over (λ)}d is the estimated variance in noise components of the audio signal, and ψ is a constant determined by the predetermined maximum amount of estimation bias.
The above described aspects of the invention may be implemented as methods or apparatus adapted to perform such methods. A computer program, stored on a computer-readable medium may cause a computer to perform any of such methods.
It is an object of the present invention to provide speech enhancement capable of estimating the relative strengths of speech and noise components that is operative during both the presence as well as the absence of speech.
It is a further object of the present invention to provide speech enhancement capable of estimating the relative strengths of speech and noise components despite the presence of a significant noise component.
It is yet a further object of the present invention to provide speech enhancement that is operative for devices having limited available memory.
These and other features and advantages of the present invention will be set forth or will become more fully apparent in the description that follows and in the appended claims. The features and advantages may be realized and obtained by means of the instruments and combinations particularly pointed out in the appended claims. Furthermore, the features and advantages of the invention may be learned by the practice of the invention or will be obvious from the description, as set forth hereinafter.



DESCRIPTION OF THE DRAWINGS
FIG. 1 is a functional block diagram showing a prior art speech enhancer.
FIG. 2a is a functional block diagram of an exemplary noise variance estimator according to aspects of the present invention. Such noise variance estimators may be used to improve prior art speech enhancers, such as that of the FIG. 1 example, or may be used for other purposes.
FIG. 2b is a flow chart useful in understanding the operation of the noise variance estimator of FIG. 2a. 
FIG. 3 shows idealized plots of estimation of bias of noise amplitude as a function of the estimated a priori SNR for four values of real SNR.



BEST MODE FOR CARRYING OUT THE INVENTION
A glossary of acronyms and terms as used herein is given in Appendix A. A list of symbols along with their respective definitions is given in Appendix B. Appendix A and Appendix B are an integral part of and form portions of the present application.
A block diagram of an exemplary embodiment of a noise variance estimator according to aspects of the invention is shown in FIG. 2a. It may be integrated with a speech enhancer such as that of FIG. 1 in order to estimate the noise level for each subband. For example, the noise variance estimator according to aspects of the invention may be employed as the Noise Variance Estimator 4 of FIG. 1, thus providing an improved speech enhancer. The input to the noise variance estimator is the unaltered subband signal Y(m) and its output is an updated value of the noise variance estimation.
For purposes of explanation, the noise variance estimator may be characterized as having three main components: a noise amplitude estimator device or function (“Estimation of Noise Amplitude”) 12, a noise variance estimate device or function that operates in response to a noise amplitude estimate (“Estimation of Noise Variance”) 14, and a speech variance estimate device or function (“Estimate of Speech Variance”) 16. The noise variance estimator example of FIG. 2a also includes a delay 18, shown using z-domain notation (“Z−1”).
The operation of the noise variance estimator example of FIG. 2a may be best understood by reference also to the flow chart of FIG. 2b. It will be appreciated that various devices, functions and processes shown and described in various examples herein may be shown combined or separated in ways other than as shown in the figures herein. For example, when implemented by computer software instruction sequences, all of the functions of FIGS. 2a and 2b may be implemented by multithreaded software instruction sequences running in suitable digital signal processing hardware, in which case the various devices and functions in the examples shown in the figures may correspond to portions of the software instructions.
The amplitude of the noise component is estimated (Estimation of Noise Amplitude 12, FIG. 2a; Estimate N(m) 24, FIG. 2b). Because the audio input signal contains both speech and noise, such estimation can only be done by exploiting statistical differences that distinguish one component from the other. Moreover, the amplitude of the noise component can be estimated via appropriate modification of existing statistical models currently used for estimation of the speech component amplitude (references [4] and [5]).
Such speech and noise models typically assume that the speech and noise components are uncorrelated, zero-mean Gaussian distributions. The key model parameters, more specifically the speech component variance and the noise component variance, must be estimated from the unaltered input audio signal. As noted above, the statistical properties of the speech and noise components are distinctly different. In most cases, the variance of the noise component is relatively stable. By contrast, the speech component is an “on/off” process and its variance can change dramatically even within several milliseconds. Consequently, an estimation of the variance of the noise component involves a relatively long time window whereas the analogous operation for the speech component may involve only current and previous input samples. An example of the latter is the “decision-directed method” proposed in reference [1].
Once the statistical models and their distribution parameters for the speech and the noise components have been determined, it is feasible to estimate the amplitudes of both components from the audio signal. In the exemplary embodiment, the Minimum Mean Square Error (MMSE) power estimator, previously introduced in reference [4] for estimating the amplitude of the speech component, is adapted to estimate the amplitude of the noise component. The choice of an estimator model is not critical to the invention.
Briefly, the MMSE power estimator first determines the probability distribution of the speech and noise components respectively based on statistical models as well as the unaltered audio signal. The noise amplitude is then determined to be the value that minimizes the mean square of the estimation error.
Finally in preparation for succeeding calculations, the variance of the noise component is updated by inclusion of the current absolute value squared of the estimated noise amplitude in the overall noise variance. This additional value becomes part of a cumulative operation on a reasonably long buffer that contains the current and as well as previous noise component amplitudes. In order to further improve the accuracy of the noise variance estimation, a Biased Estimation Avoidance method may be incorporated.
Estimation of the Noise Amplitude  
Estimation of Noise Amplitude 12, FIG. 2a; Estimate N(m) 24, FIG. 2b 
As illustrated in FIGS. 1, 2a, and 2b (20), the input to the noise variance estimator (in this context, the “noise variance estimator” is block 4 of FIG. 1 and is the combination of elements 12, 14, 16 and 18 of FIG. 2a) is the subband:

Y(m)=X(m)+D(m)  (2)

where X(m) is the speech component, and D(m) is the noise component. Here m is the time-index, and the subband number index k is omitted because the same noise variance estimator is used for each subband. One may assume that the analysis filterbank generates complex quantities, such as a DFT does. Here, the subband component is also complex, and can be further represented as and

Y(m)=R(m)exp(jθ(m))  (3)

X(m)=A(m)exp(jα(m))  (4)

and

D(m)=N(m)exp(jφ(m))  (5)

where R(m), A(m) and N(m) are the amplitudes of the unaltered audio signal, speech and noise components, respectively, and θ(m), α(m) and φ(m) are their respective phases.

By assuming that the speech and the noise components are uncorrelated, zero-mean Gaussian distributions, the amplitude of X(m) may be estimated by using the MMSE power estimator derived in reference [4] as follows:

Â(m)=GSP(ξ(m),γ(m))·R(m)  (6)

where the gain function is given by





  
    
      
        
          
            
              G
              SP
            
            ⁡
            
              (
              
                
                  ξ
                  ⁡
                  
                    (
                    m
                    )
                  
                
                ,
                
                  γ
                  ⁡
                  
                    (
                    m
                    )
                  
                
              
              )
            
          
          =
          
            
              
                
                  ξ
                  ⁡
                  
                    (
                    m
                    )
                  
                
                
                  1
                  +
                  
                    ξ
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                
              
              ⁢
              
                (
                
                  
                    1
                    +
                    
                      υ
                      ⁡
                      
                        (
                        m
                        )
                      
                    
                  
                  
                    γ
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        7
        )
      
    
  
  
    
      
        
          υ
          ⁡
          
            (
            m
            )
          
        
        =
        
          
            
              ξ
              ⁡
              
                (
                m
                )
              
            
            
              1
              +
              
                ξ
                ⁡
                
                  (
                  m
                  )
                
              
            
          
          ⁢
          
            γ
            ⁡
            
              (
              m
              )
            
          
        
      
    
    
      
        (
        8
        )
      
    
  
  
    
      
        
          
            ξ
            ⁡
            
              (
              m
              )
            
          
          =
          
            
              
                λ
                x
              
              ⁡
              
                (
                m
                )
              
            
            
              
                λ
                d
              
              ⁡
              
                (
                m
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        9
        )
      
    
  
  
    
      
        
          γ
          ⁡
          
            (
            m
            )
          
        
        =
        
          
            
              R
              2
            
            ⁡
            
              (
              m
              )
            
          
          
            
              λ
              d
            
            ⁡
            
              (
              m
              )
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




Here λx(m) and λd(m) are the variances of the speech component and noise components respectively. ξ(m) and γ(m) are often interpreted as the a priori and a posteriori component-to-noise ratios, and that notation is employed herein. In other words, the “a priori” SNR is the ratio of the assumed (while unknown in practice) speech variance (hence the name “a priori) to the noise variance. The “a posteriori” SNR is the ratio of the square of the amplitude of the observed signal (hence the name “a posterori”) to the noise variance.
In the MMSE power estimator model, the respective variances of the speech and noise components can be interchanged to estimate the amplitude of the noise component:




  
    
      
        
          
            
              N
              ^
            
            ⁡
            
              (
              m
              )
            
          
          =
          
            
              
                G
                SP
              
              ⁡
              
                (
                
                  
                    
                      ξ
                      ′
                    
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                  ,
                  
                    
                      γ
                      ′
                    
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                
                )
              
            
            ·
            
              R
              ⁡
              
                (
                m
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        11
        )
      
    
  
  
    
      
        
          
            
              ξ
              ′
            
            ⁡
            
              (
              m
              )
            
          
          =
          
            
              
                λ
                d
              
              ⁡
              
                (
                m
                )
              
            
            
              
                λ
                x
              
              ⁡
              
                (
                m
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        12
        )
      
    
  
  
    
      
        
          
            γ
            ′
          
          ⁡
          
            (
            m
            )
          
        
        =
        
          
            
              R
              2
            
            ⁡
            
              (
              m
              )
            
          
          
            
              λ
              x
            
            ⁡
            
              (
              m
              )
            
          
        
      
    
    
      
        (
        13
        )
      
    
  




Estimation of the Speech Variance  
Estimation of Speech Variance 16, FIG. 2a; Estimate {circumflex over (λ)}x(m) 22, FIG. 2b 
The estimation of the speech component variance {circumflex over (λ)}x(m) may be calculated by using the decision-directed method proposed in reference [1]:

{circumflex over (λ)}x(m),μÂ2(m−1)+(1−μ)max(R2(m)−{circumflex over (λ)}d(m),0)  (14)

Here

0<<μ<1  (15)

is a pre-selected constant, and Â(m) is the estimation of the speech component amplitude. The estimation of the noise component variance {circumflex over (λ)}d(m) calculation is described below.

Estimation of the Noise Amplitude (Continued from Above) 
The estimation of the amplitude of the noise component is finally given by




  
    
      
        
          
            
              N
              ^
            
            ⁡
            
              (
              m
              )
            
          
          =
          
            
              
                G
                SP
              
              ⁡
              
                (
                
                  
                    
                      
                        ξ
                        ^
                      
                      ′
                    
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                  ,
                  
                    
                      
                        γ
                        ^
                      
                      ′
                    
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                
                )
              
            
            ·
            
              R
              ⁡
              
                (
                m
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        16
        )
      
    
  
  
    
      
        
          
            
              
                ξ
                ^
              
              ′
            
            ⁡
            
              (
              m
              )
            
          
          =
          
            
              
                
                  λ
                  ^
                
                d
              
              ⁡
              
                (
                m
                )
              
            
            
              
                
                  λ
                  ^
                
                x
              
              ⁡
              
                (
                m
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        17
        )
      
    
  
  
    
      
        
          
            
              γ
              ^
            
            ′
          
          ⁡
          
            (
            m
            )
          
        
        =
        
          
            
              R
              2
            
            ⁡
            
              (
              m
              )
            
          
          
            
              
                λ
                ^
              
              x
            
            ⁡
            
              (
              m
              )
            
          
        
      
    
    
      
        (
        18
        )
      
    
  




Although a complex filterbank is employed in this example, it is straightforward to modify the equations for a filterbank having only real values.
The method described above is given only as an example. More sophisticated or simpler models can be employed depending on the application. Multiple microphone inputs may be used as well to obtain a better estimation of the noise amplitudes.
Estimation of the Noise Variance 
Estimation of Noise Variance 14, FIG. 2a; Estimate λd(m) 26, FIG. 2b 
The noise component in the subband input at a given time index in is, in part, determined by its variance λd(m). For a zero-mean Gaussian, this is defined as the mean value of the square of the amplitude of the noise component:

λd(m)=E{N2(m)}  (19)

Here the expectation E{N2(m)} is taken with respect to the probability distribution of the noise component at time index m.

By assuming the noise component is stationary and ergodic, λd(m) can be obtained by performing a time-averaging operation on prior estimated noise amplitudes. More specifically, the noise variance λd(m+1) of time index m+1 can be estimated by performing a weighted average of the square of the previously estimated noise amplitudes:




  
    
      
        
          
            
              λ
              ^
            
            d
          
          ⁡
          
            (
            
              m
              +
              1
            
            )
          
        
        =
        
          
            
              ∑
              
                i
                =
                0
              
              ∞
            
            ⁢
            
              
            
            ⁢
            
              
                w
                ⁡
                
                  (
                  i
                  )
                
              
              ⁢
              
                
                  
                    N
                    ^
                  
                  2
                
                ⁡
                
                  (
                  
                    m
                    -
                    i
                  
                  )
                
              
            
          
          
            
              ∑
              
                i
                =
                0
              
              ∞
            
            ⁢
            
              
            
            ⁢
            
              w
              ⁡
              
                (
                i
                )
              
            
          
        
      
    
    
      
        (
        20
        )
      
    
  




where w(i), i=0, . . . , ∞ is a weighting function. In practice w(i) can be chosen as a window of length L: w (i)=1, i=0, . . . , L−1. In the Rectangle Window Method (RWM), the estimated noise variance is given by:





  
    
      
        
          
            
              λ
              ^
            
            d
          
          ⁡
          
            (
            
              m
              +
              1
            
            )
          
        
        =
        
          
            1
            L
          
          ⁢
          
            
              ∑
              
                i
                =
                0
              
              
                L
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  N
                  ^
                
                2
              
              ⁡
              
                (
                
                  m
                  -
                  i
                
                )
              
            
          
        
      
    
    
      
        (
        21
        )
      
    
  




It is also possible to use an exponential window:

w(i)=βi+1  (22)

where

0<β<1  (23)

In the Moving Average Method (MAM), the estimated noise variance is the moving average of the square of the noise amplitudes:

{circumflex over (λ)}d(m+1)=(1−β)/{circumflex over (λ)}d(m)+β{circumflex over (N)}k2(m)  (24)

where the initial value {circumflex over (λ)}d(0) can be set to a reasonably chosen pre-determined value.

Bias Estimation Avoidance  
Occasionally, the model is unable to provide an accurate representation of the speech and noise components. In these situations, the noise variance estimation can become inaccurate, thereby producing a very biased result. The Bias Estimation Avoidance (BEA) method has been developed to mitigate this problem.
In essence, the BEA assigns a diminished weight to noise amplitude estimates {circumflex over (N)}(m) such that:

bias(m)=E{N2(m)−{circumflex over (N)}2(m)}/E{N2(m)}  (25)

where the bias, bias(m), is larger than a pre-determined maximum Bmax, i.e.:

|bias(m)|>Bmax  (26)

The accuracy of the noise amplitude estimation {circumflex over (N)}(m) is subject to the accuracy of the model, particularly the variances of the speech and the noise components as described in previous sections. Because the noise component is relatively stationary, its variance evolves slowly with time. For this reason, the analysis assumes:

{circumflex over (λ)}d(m)=λd(m)  (27)

By contrast, the speech component is transient by nature and prone to large errors. Assuming the real a priori SNR is

ξ*(m)=λx(m)/λd(m)  (28)

while the estimated a priori SNR is

{tilde over (ξ)}(m)={circumflex over (λ)}x(m)/λd(m)  (29)

the estimation bias of {circumflex over (N)}2(m) is actually given by





  
    
      
        
          
            bias
            ⁡
            
              (
              m
              )
            
          
          =
          
            
              
                
                  ξ
                  ~
                
                ⁡
                
                  (
                  m
                  )
                
              
              -
              
                
                  ξ
                  *
                
                ⁡
                
                  (
                  m
                  )
                
              
            
            
              
                (
                
                  1
                  +
                  
                    
                      ξ
                      ~
                    
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                
                )
              
              2
            
          
        
        ⁢
        
          

        
        ⁢
        
          Clearly
          ,
          if
        
      
    
    
      
        (
        30
        )
      
    
  
  
    
      
        
          
            ξ
            ~
          
          ⁡
          
            (
            m
            )
          
        
        =
        
          
            ξ
            *
          
          ⁡
          
            (
            m
            )
          
        
      
    
    
      
        (
        31
        )
      
    
  




one has an unbiased estimator and

E{{circumflex over (N)}2(m)}=E{N2(m)}=λd(m)  (32)

As seen in FIG. 3, the estimation bias is asymmetric with respect to the dotted line in the figure, the zero bias line. The lower portion of the plot indicates widely varying values of the estimation bias for varying values of ξ* whereas the upper portion shows little dependency on either {tilde over (ξ)} or ξ*.
For the SNR range of interest, under-estimation of noise amplitude, i.e.:

E{{circumflex over (N)}2(m)}<E{N2(m)}  (33)

will result in a positive bias, corresponding to the upper portion of the plot. As can be seen, the effect is relatively small and therefore not problematic.

The lower portion of the plot, however, corresponds to cases wherein the variance of the speech component is underestimated, resulting in a large negative estimation bias as given by Eqn. (30), i.e.:

λx(m)>{circumflex over (λ)}x(m)  (34)

and

λd(m)>{circumflex over (λ)}x(m)  (35)

or, alternatively

ξ*(m)>{tilde over (ξ)}(m)  (36)

and

{tilde over (ξ)}(m)<1  (37)

as well as a strong dependency on different values of ξ*. These are situations in which the estimate of the noise amplitude is too large. Consequently, such amplitudes are given diminished weight or avoided altogether.

In practice, experience has taught that such suspect amplitudes R(m) satisfy:

R2(m)>ψ(1+{circumflex over (ξ)}(m))λd(m)  (38)

where ψ is a predefined positive constant. This rule provides a lower bound for the bias:





  
    
      
        
          
            bias
            ⁡
            
              (
              m
              )
            
          
          >
          
            1
            -
            
              
                1
                2
              
              ⁢
              ψ
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        39
        )
      
    
  
  
    
      
         
        ⁢
        
          ψ
          =
          
            2
            ⁢
            
              (
              
                
                  B
                  max
                
                +
                1
              
              )
            
          
        
      
    
    
      
        (
        40
        )
      
    
  




In summary, a positive bias is negligible. A negative bias is tenable if estimated noise amplitudes {circumflex over (N)}(m) defined in Eqn. (16) and consistent with Eqn. (38) are given diminished weight. In practical application, since the value of λd(m) is unknown, the rule of Eqn. (38) can be approximate by:




  
    
      
        
          
            
              R
              2
            
            ⁡
            
              (
              m
              )
            
          
          >
          
            
              ψ
              (
              
                1
                +
                
                  
                    ξ
                    ^
                  
                  ⁡
                  
                    (
                    m
                    )
                  
                
              
              )
            
            ⁢
            
              
                
                  λ
                  ^
                
                d
              
              ⁡
              
                (
                m
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        41
        )
      
    
  
  
    
      
        
          
            ξ
            ^
          
          ⁡
          
            (
            m
            )
          
        
        =
        
          
            
              
                λ
                ^
              
              x
            
            ⁡
            
              (
              m
              )
            
          
          
            
              
                λ
                ^
              
              d
            
            ⁡
            
              (
              m
              )
            
          
        
      
    
    
      
        (
        42
        )
      
    
  




Two such examples of the BEA method are the Rectangle Window Method (RWM) with BEA and the Moving Average Method (MAM) with BEA. In the former implementation, weight given to samples that are consistent with Eqn. (38) is zero:




  
    
      
        
          
            
              λ
              ^
            
            d
          
          ⁡
          
            (
            
              m
              +
              1
            
            )
          
        
        =
        
          
            1
            L
          
          ⁢
          
            
              ∑
              
                i
                ∈
                
                  Φ
                  m
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  N
                  ^
                
                2
              
              ⁡
              
                (
                i
                )
              
            
          
        
      
    
    
      
        (
        43
        )
      
    
  




where Φm is a set that contains L nearest {circumflex over (N)}2(i) to time index m that satisfy

R2(i)≦ψ(i+{circumflex over (ξ)}(i)){circumflex over (λ)}d(i)  (44)

In the latter implementation, such samples may be included with a diminished weight:




  
    
      
        
          
            
              
                λ
                ^
              
              d
            
            ⁡
            
              (
              
                m
                +
                1
              
              )
            
          
          =
          
            
              
                (
                
                  1
                  -
                  β
                
                )
              
              ⁢
              
                
                  
                    λ
                    ^
                  
                  d
                
                ⁡
                
                  (
                  m
                  )
                
              
            
            +
            
              β
              ⁢
              
                
                  
                    N
                    ^
                  
                  k
                  2
                
                ⁡
                
                  (
                  m
                  )
                
              
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        45
        )
      
    
  
  
    
      
        β
        =
        
          {
          
            
              
                
                  
                    β
                    0
                  
                
                
                  
                    
                      
                        R
                        2
                      
                      ⁡
                      
                        (
                        m
                        )
                      
                    
                    ≤
                    
                      
                        ψ
                        (
                        
                          1
                          +
                          
                            
                              ξ
                              ^
                            
                            ⁡
                            
                              (
                              m
                              )
                            
                          
                        
                        )
                      
                      ⁢
                      
                        
                          
                            λ
                            ^
                          
                          d
                        
                        ⁡
                        
                          (
                          m
                          )
                        
                      
                    
                  
                
              
              
                
                  
                    β
                    1
                  
                
                
                  
                    else
                    .
                  
                
              
            
            ⁢
            
              

            
            ⁢
            and
          
        
      
    
    
      
        (
        46
        )
      
    
  
  
    
      
        
          β
          1
        
        <
        
          β
          0
        
      
    
    
      
        (
        47
        )
      
    
  




Completing the description of the FIG. 2b flowchart, the time index m is then advanced by one (“m←m+1” 56) and the process of FIG. 2b is repeated.
Implementation
The invention may be implemented in hardware or software, or a combination of both (e.g., programmable logic arrays). Unless otherwise specified, the processes included as part of the invention are not inherently related to any particular computer or other apparatus. In particular, various general-purpose machines may be used with programs written in accordance with the teachings herein, or it may be more convenient to construct more specialized apparatus (e.g., integrated circuits) to perform the required method steps. Thus, the invention may be implemented in one or more computer programs executing on one or more programmable computer systems each comprising at least one processor, at least one data storage system (including volatile and non-volatile memory and/or storage elements), at least one input device or port, and at least one output device or port. Program code is applied to input data to perform the functions described herein and generate output information. The output information is applied to one or more output devices, in known fashion.
Each such program may be implemented in any desired computer language (including machine, assembly, or high level procedural, logical, or object oriented programming languages) to communicate with a computer system. In any case, the language may be a compiled or interpreted language.
Each such computer program is preferably stored on or downloaded to a storage media or device (e.g., solid state memory or media, or magnetic or optical media) readable by a general or special purpose programmable computer, for configuring and operating the computer when the storage media or device is read by the computer system to perform the procedures described herein. The inventive system may also be considered to be implemented as a computer-readable storage medium, configured with a computer program, where the storage medium so configured causes a computer system to operate in a specific and predefined manner to perform the functions described herein.
A number of embodiments of the invention have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the invention. For example, some of the steps described herein may be order independent, and thus can be performed in an order different from that described.
APPENDIX A
Glossary of Acronyms and Terms


    BEA Biased Estimation Avoidance
    DFT Discrete Fourier Transform
    DSP Digital Signal Processing
    MAM Moving Average Method
    RWM Rectangle Window Method
    SNR Signal to Noise ratio
    T/F time/frequency
    VAD Voice Activity Detection


APPENDIX B
List of Symbols


    y(n), n=0, 1, . . . , ∞ digitized time signal
    {tilde over (y)}(n) enhanced speech signal
    Yk(m), k=1, . . . , K, m=0, 1, . . . , ∞ subband signal k
    {tilde over (Y)}k(m) enhanced subband signal k
    X(m) speech component of subband k
    D(m) noise component of subband k
    gk suppression gain for subband k
    R(m) noisy speech amplitude
    θ(m) noisy speech phase
    A(m) speech component amplitude
    Â(m) estimated speech component amplitude
    α(m) speech component phase
    N(m) noise component amplitude
    {circumflex over (N)}(m) estimated noise component amplitude
    φ(m) noise component phase
    GSP gain function
    λx(m) speech component variance
    {circumflex over (λ)}x(m) estimated speech component variance
    λd(m) noise component variance
    {circumflex over (λ)}d(m) estimated noise component variance
    ξ(m) a priori speech component-to-noise ratio
    γ(m) a posteriori speech component-to-noise ratio
    ξ′(m) a priori noise component-to-speech ratio
    γ′(m) a posteriori noise component-to-speech ratio
    α pre-selected constant
    β pre-selected for bias estimation

## Claims

1. A method for enhancing speech components of an audio signal composed of speech and noise components, comprising
transforming the audio signal from the time domain to a plurality of subbands in the frequency domain,
wherein each of said plurality of subbands is presumed to have a speech component and a noise component, said noise component having an amplitude and a variance at time index m, wherein said amplitude of the noise component is estimated by exploiting statistical differences that distinguish between the speech component and the noise component,

processing each of said plurality of subbands, said processing including applying a gain factor, wherein said gain factor is derived at least in part from an estimation of said variance in noise components, wherein the estimation comprises
at each time index m, updating said estimation of variance in noise components of the subband signal from an average of past estimates of the amplitude of noise components in the subband signal, and
wherein said past estimates of the amplitude of noise components in the subband signal having values greater than a threshold are excluded from or underweighted in said weighted average, and

transforming the processed subband signal from the frequency domain to the time domain to provide an audio signal in which speech components are enhanced.

2. A method according to claim 1 wherein the average of past estimates of the amplitude of noise components is a weighted average of the square of the past estimate of the amplitude of a noise component and the past estimated variance in noise components.

3. A method according to claim 2 wherein the weighting function of the weighted average is a preselected constant.

4. A method according to claim 1 wherein each estimate of the amplitude of noise components in the subband signal is a function of an estimate of variance in noise components of the subband signal, an estimate of variance in speech components of the subband signal, and the amplitude of the subband signal.

5. A method according to claim 1 wherein said threshold is a function of ψ(1+{circumflex over (ξ)}(m)){circumflex over (λ)}d(m), where {circumflex over (ξ)} is the estimated a priori signal-to-noise ratio, {circumflex over (λ)}d is the estimated variance in noise components of the subband signal, and ψ is a constant determined by a predetermined maximum amount of an estimation bias.

6. A method according to claim 5 wherein each estimate of the amplitude of noise components in the subband signal is a function of an estimate of variance in noise components of the subband signal, an estimate of variance in speech components of the subband signal, and the amplitude of the subband signal.

7. Apparatus adapted to perform the methods of any one of claims 1 through 6.

8. A non-transitory computer-readable storage medium encoded with a computer program for causing a computer to perform the method of any one of claims 1 through 6.