# Patent

## Title
Constructing broad-band acoustic signals from lower-band acoustic signals

## Abstract

A method generates envelope spectra and harmonic spectra from an input broad-band training acoustic signal. Corresponding non-negative envelope bases are trained for the envelope spectra and non-negative harmonic bases are trained for the harmonic spectra using convolutive non-negative matrix factorization. Higher-band frequencies are generated for an input lower-band acoustic signal according to the non-negative envelope bases and the non-negative harmonic bases. Then, the input lower-band acoustic signal is combined with the higher-band frequencies to produce an output broad-band acoustic signal.

## Background
FIELD OF THE INVENTION
This invention relates generally to processing acoustic signals, and more particularly to constructing broad-band acoustic signals from lower-band acoustic signals.
BACKGROUND OF THE INVENTION
Broad-band acoustic signals, e.g., speech signals that contain frequencies from a range of approximately 0 kHz to 8 kHz are naturally better sounding and more intelligible than lower-band acoustic signals that have frequencies approximately less than 4 kHz, e.g., telephone quality acoustic. Therefore, it is desired to expand lower-band acoustic signals.
Various methods are known to solve this problem. Aliasing-based methods derive high-frequency components by aliasing low frequencies into high frequencies by various means, Yasukawa, H., “Signal Restoration of Broad Band Speech Using Nonlinear Processing,” Proc. European Signal Processing Conf. (EUSIPCO-96), pp. 987-990, 1996.
Codebook methods map a spectrum of the lower-band speech signal to a codeword in a codebook, and then derive higher frequencies from a corresponding high-frequency codeword, Chennoukh, S., Gerrits, A., Miet, G. and Sluijter, R., “Speech Enhancement via Frequency Bandwidth Extension using Line Spectral Frequencies,” Proc ICASSP-95, 2001.
Statistical methods utilize the statistical relationship of lower-band and higher-band frequency components to derive the latter from the former. One method models the lower-band and higher-band components of speech as mixtures of random processes. Mixture weights derived from the lower-band signals are used to generate the higher-band frequencies, Cheng, Y. M., O'Shaugnessey, D. O., and Mermelstein, P., “Statistical Recovery of Wideband Speech from Narrow-band Speech,” IEEE Trans., ASSP, Vol 2., pp 544-548, 1994.
Methods that use statistical cross-frame correlations can predict higher frequencies. However, those methods are often derived from complex time-series models, such as Gaussian mixture models (GMMs), hidden Markov models (HMMs) or multi-band HMMs, or by explicit interpolation, Hosoki, M., Nagai, T. and Kurematsu, A., “Speech Signal Bandwidth Extension and Noise Removal Using Subband HIGHER-BAND,” Proc. ICASSP, 2002.
Linear model methods derive higher-band frequency components as linear combinations of lower-band frequency components, Avendano, C., Hermansky, H., and Wand, E. A., “Beyond Nyquist: Towards the Recovery of Broad-bandwidth Speech from Narrow-bandwidth Speech,” Proc. Eurospeech-95, 1995.
SUMMARY OF THE INVENTION
A method estimates high frequency components, e.g., approximately a range of 4-8 kHz, of acoustic signals from lower-band, e.g., approximately a range of 0-4 kHz, acoustic signals using a convolutive non-negative matrix factorization (CNMF).
The method uses input training broad-band acoustic signals to train a set of lower-band and corresponding higher-band non-negative ‘bases’. The acoustic signals can be, for example, speech or music. The low-frequency components of these bases are used to determine high-frequency components and can be combined with an input lower-band acoustic signal to construct an output broad-band acoustic signal. The output broad-band acoustic signal is virtually indistinguishable from a true broad-band acoustic signal.



BRIEF DESCRIPTION OF THE DRAWING
FIG. 1 is a block diagram of a method for expanding an acoustic signal according to one embodiment of the invention.



DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
Convolutive Non-Negative Matrix Factorization
Matrix factorization decomposes a matrix V into two matrices W and H, such that:

V≈W·H,  (1)

where W is an M×R matrix, H is a R×N matrix, and R is less than M, while an error of reconstruction of the matrix V from the matrices Wand H is minimized. In such a decomposition, the columns of the matrix W can be interpreted as a set of bases, and the columns of the matrix H as the coordinates of the columns of V, in terms of the bases.

Alternately, the columns of the matrix H represent weights with which the bases in the matrix W are combined to obtain a closest approximation to the columns of the matrix V.
Conventional factorization techniques, such as principal component analysis (PCA) and independent component analysis (ICA), allow the bases to be positive and negative, and the interaction between the terms, as specified by the components of the matrix H, can also be positive and negative.
In strictly non-negative data sets such as matrices that represent sequences of magnitude spectral vectors, neither negative components in the bases nor negative interaction are allowed because the magnitudes of spectral vectors cannot be negative.
One non-negative matrix factorization (NMF) constrains the elements of the matrices W and H to be strictly non-negative, Lee, D. D and H. S. Seung. “Learning the parts of objects with nonnegative matrix factorization,” Nature 401, pp. 788-791, 1999. They apply NMF to detect parts of faces in hand-aligned 2D images, and semantic features of summarized text. Another application applies NMF to detect individual notes in acoustic recordings of musical pieces, P. Smaragdis, “Discovering Auditory Objects Through Non-Negativity Constraints,” SAPA 2004, October 2004.
The NMF of Lee et al. treats all column bases in the matrix V as a combination of R bases, and assumes implicitly that it is sufficient to explain the structure within individual bases to explain the entire data set. This effectively assumes that the order in which the bases are arranged in the matrix V is irrelevant.
However, these assumptions are clearly invalid in data sets such as sequences of magnitude spectral bases, where structural patterns are evident across multiple bases, and an order in which the bases are arranged is indeed relevant.
Smaragdis describes a convolutive version of the NMF algorithm (CNMF), wherein the bases used to explain the matrix V are not merely singular bases, but actually short sequences of bases. This operation can be symbolically represented as:




  
    
      
        
          V
          ≈
          
            
              ∑
              
                t
                =
                0
              
              τ
            
            ⁢
            
              
            
            ⁢
            
              
                W
                t
                T
              
              ·
              
                H
                
                  t
                  →
                  T
                
              
            
          
        
        ,
      
    
    
      
        (
        2
        )
      
    
  




where each WtT is a non-negative M×R matrix, H is a non-negative R×N matrix, as above, the (t→) operator represents a right shift operator that shifts the columns of matrix H by t positions to the right. The T in the superscript of Equation 2 represents a transposition operator. The size of the matrix H is maintained by introducing zero valued columns at the leftmost position to account for columns that have been shifted out of the matrix.

We represent the jth vector in Wt as Wtj. Each set of vectors forms a sequence of spectral vectors wj, or a ‘spectral patch’ in an acoustic signal, e.g., a speech or music signal. These spectral patches form the bases that we use to ‘explain’ the data in the matrix V.
Equation 2 approximates the matrix V as a superposition of the convolution of these patches with the corresponding rows of the matrix H, i.e., the contribution of jth spectral patch to the approximation of the matrix V is obtained by convolving the patch with the jth row of the matrix H.
If τ=1, then this reduces to the conventional NMF. To estimate the appropriate matrices Wt, and the matrix H to estimate the matrix V, we can use the already existing framework of NMF.
We define a cost function as:




  
    
      
        
          D
          =
          
            
              
              
                
                  V
                  ⊗
                  
                    ln
                    ⁡
                    
                      (
                      
                        V
                        Λ
                      
                      )
                    
                  
                
                +
                Λ
                -
                V
              
              
            
            F
          
        
        ,
      
    
    
      
        (
        3
        )
      
    
  




where the norm on the right side is a Froebinus norm, {circle around (x)} represents a Hadmard component by component multiplication, Λ is the current reconstruction given by the right hand side of Equation 2, using the current estimates of H and the Wt matrices, and F is a lower cutoff frequency, e.g. 4000 Hz. The matrix division to the right is also per-component, and is the approximation to the matrix V given by the right hand side of Equation 2.

The cost function of Equation 3 is a modified Kullback-Leibler cost function. Here, the approximation is given by the convolutive NMF decomposition of Equation 2, instead of the linear decomposition of Equation 1.
Equation 2 can also be viewed as a set of NMF operations that are summed to produce the final result. From this perspective, the chief distinction between Equations 1 and 2 is that the latter decomposes the matrix V into a combination of τ+1 matrices, while the former uses only two matrices.
This interpretation permits us to obtain an iterative procedure for the estimation of the matrices Wt and H matrices by modifying the NMF update equations of Lee et al. The modified iterative update equations are given by:




  
    
      
        H
        =
        
          H
          ⊗
          
            
              
                ∑
                t
              
              ⁢
              
                
                  W
                  t
                  T
                
                ·
                
                  [
                  
                    
                      V
                      
                        ←
                        t
                      
                    
                    Λ
                  
                  ]
                
              
            
            
              
                ∑
                t
              
              ⁢
              
                
                  W
                  t
                  T
                
                ·
                1
              
            
          
        
      
    
    
      
        (
        4
        )
      
    
  
  
    
      
        
          W
          t
        
        =
        
          
            W
            t
          
          ⊗
          
            
              
                [
                
                  V
                  Λ
                
                ]
              
              ·
              
                H
                
                  t
                  →
                  T
                
              
            
            
              1
              ·
              
                H
                
                  t
                  →
                  T
                
              
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




where {circle around (x)} represents a component-by-component Hadamard multiplication, and the division operations are also component-into-component. The (←t) operator represents a left shift operator, the inverse of to the right shift operator in Equation 2. The overall procedure for estimating the Wt and H matrices, thus, is as follows:

Initialize all matrices, e.g., use a random initialization, thereafter iteratively update all terms using Equations 4 and 5.
The spectral patches Wj, comprising the jth columns of all the matrices Wtj trained by the CNMF, represent salient spectrographic structures in the acoustic signal.
When applied to speech signals as described below, the trained bases represent relevant phonemic or sub-phonetic structures.
Constructing High Frequency Structures of a Band Limited Acoustic Signal
As shown in FIG. 1, a method 100 for constructing higher-band frequencies for a narrow-band signal includes the following components:
A signal processing component 110 generates, from an input broad-band training acoustic signal 101, representations for low-resolution spectra and high-resolution spectra, hereinafter ‘envelope spectra’ 111, and the ‘harmonic spectra’ 112, respectively.
A training component 120 trains corresponding non-negative envelope bases 121 for the envelope spectra, and non-negative harmonic bases 122 for the harmonic spectra using the convolutive non-negative matrix factorization.
A construction component 130 constructs higher-band frequencies 131 for an input lower-band acoustic signal 132, which are then combined 140 to produce an output broad-band acoustic signal 141.
Signal Processing
A sampling rate for all of the acoustic signals is sufficient to acquire both lower-band and higher-band frequencies. Signals sampled at lower frequencies are upsampled to this rate. We use a sampling rate of 16 kHz, and all window sizes and other parameters described below are given with reference to this sampling rate.
We determine a short-time Fourier transform of the acoustic signals using a Hanning window of 512 samples (32 ms) for each frame, with an overlap of 256 samples between adjacent frames, timed-synchronously with the corresponding input broad-band training acoustic signal.
A matrix S represent a sequence of complex Fourier spectra for the acoustic signal, a matrix Φ represent the phase, and a matrix V represents the component-wise magnitude of the matrix S. Thus, the matrix V represents the magnitude spectrogram of the signal.
In the matrices V and Φ, each column represents respectively the magnitude spectra and phase of a single 32 ms frame of the acoustic signal. If there are M unique samples in the Fourier spectrum for each frame, and there are N frames in the signal, then the matrices V and Φ are M×N matrices.
We determine the envelope spectra 111 and the harmonic spectra 112 of the training acoustic signal 101 by cepstral weighting or ‘liftering’ the matrix V. The matrix Ve represents the sequence of envelope spectra derived from the matrix V, and the matrix Vh represents the sequence of corresponding harmonic spectra. The matrices Ve and Vh are both M×N matrices derived from the matrix V according to:

Vh=exp(IDCT(DCT((log(V)){circle around (x)}Zh)))  (6)

Ve=exp(IDCT(DCT((log(V)){circle around (x)}Ze)))  (7)

The matrix Ze has the lower K frequency components of each row are set to one, and the rest of the frequency components are set to zero. The matrix Zh has the higher frequency components set to one and the rest of the frequency components set to zero, i.e.,

Zh=1−Ze.

The discrete cosine transform (DCT) and the inverse DCT operations in Equations 6 and 7 are applied separately to each row of the respective matrix arguments.
With an appropriate selection of the lower frequency K components, e.g., K=M/3, the matrices Ve and Vh model the structure of the envelope spectra and harmonic spectra of the training signal 101.
Lower frequencies of the envelope spectra of the lower-band portion of the training acoustic signal, and upper frequencies of the envelope spectra of the training acoustic signal can be combined to compose a synthetic envelope spectral matrix. Similarly, lower frequencies of the harmonic spectra of the lower-band training signal, and upper frequencies of the harmonic spectra of the input broad-band training signal can be combined to compose a synthetic harmonic spectral matrix.
Training Spectral Bases
The first stage of the training step 120 trains the matrices Ve, Vh, and Φ from the training signal 101. The training signal can be speaker dependent or speaker independent, because characteristics of any speaker or group of speakers can be acquired by relatively short signals, e.g., five minutes or less.
The matrices are obtained in a two-step process. In the first step, the training signal is filtered to a frequency band expected in the lower-band acoustic signal 132, and then down-sampled to an expected sampling rate of the lower-band signal 132, and finally upsampled to the sampling rate of the higher-band signal 131. This signal is a close approximation to the signals that is obtained by up-sampling the lower-band signal.
Harmonic, envelope and phase spectral matrices Vhn, Ven, and Φn are obtained from the upsampled lower-band training signal.
Envelope, harmonic and phase spectral matrices Vew, Vhw and Φw are derived from the wide-band training signal 101. The matrices Vh, Ve and Φ are formed from frequency components less than a predetermined cutoff frequency F, from the spectral matrices for the lower-band, and the higher frequency components of the matrices derived from the broad-band signal as:

Ve=ZwVew+ZnVen 

Vh=ZwVhw+ZhVen 

Φ=ZwΦw+ZnΦn  (8)

The matrix Zw is a square matrix with the first diagonal elements set to one and the remaining elements set to zero. The matrix Zn is also a square matrix with the last diagonal elements set to one and the remaining elements set to zero. The parameter L is a frequency index that corresponds to the cutoff frequency F.
The spectral patch bases Wte for t=1, . . . , τe are derived for the envelope spectra Ve using the iterative update process specified by Equations 4 and 5. The matrix H is discarded.
The set of lower-band spectral envelope bases, Wte,l derived from the envelope spectra Ve, are obtained by truncating all the matrices at the Lth row, such that each of the resulting matrices is of size L×R:

Wte,l=ZLWte  (9)

The matrix ZL is a L×M matrix, where the L leading diagonal elements are one, and the remaining elements are zero.

The set of lower-band spectral harmonic bases, Wth,l are obtained similarly. The set of matrices, Wte, Wtl,t, Wth form the spectral patch bases to be used for construction.
The phase matrix Φ is separated into a L×N low-frequency phase matrix Φl and a M−(L×N) high-frequency matrix Φu.
A linear regression between the matrices is obtained:

AΦ=Φu·pseudoinverse(Φh)  (10)

Constructing Broad-Band Acoustic Signals
The input lower-band acoustic signal 132 is upsampled to the sampling rate of the broad-band training signal 101, and the phase, envelope and harmonic spectral matrices Φ, Vh, and Ve, are derived from upsampled signal. The lower frequency components of the matrices are separated out as Ve=ZLVe and Vh=ZLVh.
CNMF approximations are obtained for the matrices Vel and Vhl, based on the Wte,l and Wth,l bases obtained from the training signal. This approximates Vel and Vhl as:




  
    
      
        
          V
          h
          l
        
        ≈
        
          
            ∑
            
              t
              =
              0
            
            
              τ
              h
            
          
          ⁢
          
            
          
          ⁢
          
            
              
                
                  (
                  
                    W
                    t
                    
                      h
                      ,
                      l
                    
                  
                  )
                
                T
              
              ·
              
                
                  (
                  
                    H
                    h
                  
                  )
                
                
                  t
                  →
                  T
                
              
            
            ⁢
            
              
            
            ⁢
            and
            ⁢
            
              
            
            ⁢
            
              V
              e
              l
            
          
        
        ≈
        
          
            ∑
            
              t
              =
              0
            
            
              τ
              e
            
          
          ⁢
          
            
          
          ⁢
          
            
              
                (
                
                  W
                  t
                  
                    e
                    ,
                    l
                  
                
                )
              
              T
            
            ·
            
              
                (
                
                  H
                  e
                
                )
              
              
                t
                →
                T
              
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




The Hh and He matrices are obtained through iterations of Equation 4.
Then, broad-band spectrograms are constructed by applying the estimated matrices Hh and He to the complete bases Wte and Wth obtained by the training:




  
    
      
        
          
            V
            _
          
          h
        
        =
        
          
            
              ∑
              
                t
                =
                0
              
              
                τ
                h
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  
                    (
                    
                      W
                      t
                      h
                    
                    )
                  
                  T
                
                ·
                
                  
                    (
                    
                      H
                      h
                    
                    )
                  
                  
                    t
                    ->
                    T
                  
                
              
              ⁢
              
                
              
              ⁢
              and
              ⁢
              
                
              
              ⁢
              
                
                  V
                  _
                
                e
              
            
          
          =
          
            
              ∑
              
                t
                =
                0
              
              
                τ
                e
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  (
                  
                    W
                    t
                    e
                  
                  )
                
                T
              
              ·
              
                
                  (
                  
                    H
                    e
                  
                  )
                
                
                  t
                  ->
                  T
                
              
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




The higher-band frequencies 131 and input lower-band frequencies 132 are obtained according to:

{circumflex over (V)}h=Zw Vh+ZhVh and {circumflex over (V)}e=Zw Ve+ZeVe.  (13)

The complete magnitude spectrum for the output broad-band signal 141 is obtained as a combination (C):

{circumflex over (V)}={circumflex over (V)}h{circle around (x)}{circumflex over (V)}e.

A phase for output the broad-band signal is:

{circumflex over (Φ)}=(Zh+ZUAΦZL)  (14)

where ZU is a M×L matrix, with (M−L) leading diagonal elements set to one, and the remaining elements set to zero.

Then, the complete output broad-band signal 141 is obtained by determining an inverse short-time Fourier transform of {circumflex over (V)}ejΦ.
Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications may be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.

## Claims

1. A method for constructing a broad-band acoustic signal from a lower-band acoustic signal, comprising:
generating envelope spectra and harmonic spectra from an input broad-band training acoustic signal;
generating corresponding non-negative envelope bases for the envelope spectra and non-negative harmonic bases for the harmonic spectra using convolutive non-negative matrix factorization;
generating higher-band frequencies for an input lower-band acoustic signal according to the non-negative envelope bases and the non-negative harmonic bases; and
combining the input lower-band acoustic signal with the generated higher-band frequencies to produce an output broad-band acoustic signal.

2. The method of claim 1, in which the input broad-band training acoustic signal and the input lower-band acoustic signal are speaker dependent.

3. The method of claim 1, in which the input broad-band training acoustic signal and the input lower-band acoustic signal are speaker independent.

4. The method of claim 1, in which the input broad-band training acoustic band signal and the output broad-band acoustic signal include frequencies in a range of approximately 0 khZ to 8 kHz, and the input lower-band acoustic signal includes frequencies in a range of approximately 0 kHz to 4 kHz, and the higher-band acoustic signal includes frequencies approximately in a range of 4 kHz to 8 kHz.

5. The method of claim 1, in which a sampling rate for the input broad-band training acoustic signal is sufficient to acquire both the lower-band and higher-band frequencies.

6. The method of claim 5, in which the input broad-band training signal is low-pass filtered to a frequency expected in the lower-band acoustic signal, and further comprising:
downsampling the low-pass filtered signal to a lower sampling rate; and
upsampling the downsampled signal back to the sampling rate of the input broadband training acoustic signal, to generate a lower-band training acoustic signal.

7. The method of claim 5, further comprising:
determining a short-time Fourier transform of the input broad-band training acoustic signal using a Hanning window of 512 samples for each frame, with an overlap of 256 samples between adjacent frames, and in which, for the input broad-band training acoustic signal, a matrix S represents a sequence of complex Fourier spectra, a matrix Φw represents a phase, and a matrix Vw represents a component-wise magnitude of the matrix S such that the matrix Vw represents a magnitude spectrogram of the input broad-band training acoustic signal.

8. The method of claim 7, in which the input broad-band training acoustic signal includes M unique samples in the Fourier spectrum for each frame, and there are N frames in the an input broad-band training acoustic signal, and the matrices Vw and Φw are M×N matrices.

9. The method of claim 8, further comprising:
determining the envelope spectra and the harmonic spectra of the input broad-band training acoustic signal by cepstral weighting of the matrix Vw.

10. The method of claim 6, further comprising:
determining a short-time Fourier transform of the lower-band training acoustic signal using a Hanning window of 512 samples for each frame, with an overlap of 256 samples between adjacent frames, timed-synchronously with the corresponding input broad-band training acoustic signal.

11. The method of claim 10, in which the input lower-band training acoustic signal includes M unique samples in a Fourier spectrum for each frame, and there are N frames in the lower-band training acoustic signal, resulting in an M×N spectral matrix, from which a matrix Φn representing a phase, and a matrix Vn representing a component-wise magnitude are derived.

12. The method of claim 11, further comprising:
determining the envelope spectra and the harmonic spectra of the lower-band training acoustic signal by cepstral weighting of the matrix Vn.

13. The method of claims 9 or 12, further comprising:
combining lower frequencies of the envelope spectra of the lower-band training acoustic signal, and upper frequencies of the envelope spectra of the input broad-band training acoustic signal to compose a synthetic envelope spectral matrix.

14. The method of claim 13, further comprising:
learning non-negative envelope bases for the synthetic envelope spectral matrix.

15. The method of claims 9 or 12, further comprising:
combining lower frequencies of the harmonic spectra of the lower-band training signal, and upper frequencies of the harmonic spectra of the input broad-band training signal to compose a synthetic harmonic spectral matrix.

16. The method of claim 15, further comprising:
learning non-negative harmonic bases for the synthetic harmonic spectral matrix.

17. The method of claims 8 or 11, in which a linear transformation AΦ is determined between lower frequencies of the matrix Φw and upper frequencies of the matrix Φw.

18. The method of claim 1, further comprising:
upsampling the input lower-band acoustic signal to a sampling frequency of the input broad-band training acoustic signal.

19. The method of claim 18, further comprising
determining a short-time Fourier transform of the input lower-band acoustic signal using a Hanning window of 512 samples for each frame, with an overlap of 256 samples between adjacent frames to generate a Fourier spectral matrix; and
deriving an envelope spectrum and a harmonic spectrum from the Fourier spectral matrix by cepstral weighting.

20. The methods of claim 14, further comprising:
deriving optimal weights of the non-negative envelope bases from the envelope spectrum of the input lower-band acoustic signal.

21. The method of claim 20, further comprising:
combining the upper frequencies of the envelope bases with the optimal weights to derive a reconstructed upper-frequency envelope spectrum.

22. The method of claim 16, further comprising:
deriving optimal weights of the non-negative harmonic bases from the harmonic spectrum of the input lower-band acoustic signal.

23. The method of claim 22, further comprising:
combining the upper frequencies of the harmonic bases with the optimal weights to derive a reconstructed upper-frequency harmonic spectrum.

24. The method of claim 21, further comprising:
multiplying the reconstructed upper-frequency envelope and harmonic spectra to derive a reconstructed upper-frequency magnitude spectrum.

25. The methods of claims 17, further comprising:
multiplying a phase of the lower frequencies of the lower-band signal by the linear transformation AΦ to derive a reconstructed phase of the upper-frequency magnitude spectrum.

26. The methods of 24, further comprising:
combining the reconstructed phase and magnitude of the upper-frequency magnitude spectrum;
determining an inverse Fourier transform to derive the upper frequency signal; and
combining the upper frequency signal with the input lower-band signal to produce an output broad-band acoustic signal.