# Patent

## Title
Vehicle engine sound extraction and reproduction

## Abstract

An RPM signal representing a current RPM of a vehicle engine is received. An instant magnitude and an instant phase are determined for each of N harmonics of a target engine sound based on the RPM signal. The instant magnitude and the instant phase are each smoothed to generate a smoothed harmonic magnitude and a smoothed harmonic phase for each of the N harmonics. A magnitude perturbation signal and a phase perturbation signal are generated for each of the N harmonics. An individual harmonic signal is also generated for each of the N harmonics based, at least in part, on the smoothed harmonic magnitude, the magnitude perturbation signal, the smoothed harmonic phase, and the phase perturbation signal. The individual harmonic signals are summed to generate an engine harmonic enhancement signal, and the engine harmonic enhancement signal is transduced to acoustic energy.

## Background
BACKGROUND
This disclosure relates to reproducing a character of an engine sound in a vehicle via engine harmonic enhancement. This method can provide a natural sound for applications in which a high level of enhancement is required.
Some known engine harmonic enhancement methods used harmonic generation based on magnitude only. In some cases, an initial phase difference between the harmonics was introduced. However, such methods may be insufficient to generate a natural sounding, high level of enhancement, or to be able to reproduce a particular engine sound character.
SUMMARY
This disclosure is based, in part, on the realization that the respective phase of each of a plurality of harmonics (e.g., the 0.5 harmonic order through 40th harmonic order) over a RPM range (e.g., 600 RPM to 7500 RPM) can be taken into account during harmonic generation in an engine harmonic enhancement (EHE) process. At the same time small variations in phase and magnitude, which confer naturalness to a generated engine sound, can be preserved by introducing a distortion signal in the generation process. This distortion signal has the same standard deviation as the variations identified in the analysis of the original engine sound.
In one aspect, a method includes performing a harmonic decomposition on a target engine sound, thereby to extract each of N harmonics of the target engine sound over an RPM range; and out of the extracted harmonics, extracting phase and shape information for each of the N harmonics over the RPM range for reproducing the target engine sound. The method also includes configuring an engine harmonic enhancement (EHE) system to utilize the extracted phase and shape information to generate an engine harmonic enhancement signal, he(t), to be added on top of a baseline engine sound.
Implementations may include one of the following features, or any combination thereof.
In some implementations, performing the harmonic decomposition includes obtaining a set of frequency weighting coefficients ak(t) and bk(t) for each of the N harmonics over the RPM range. It is considered that as time is passing the RPM is increasing, such as at time t=0 the engine is at the lowest RPM and the last moment in time corresponds to the highest RPM.
In certain implementations, performing the harmonic decomposition includes utilizing an adaptive algorithm to obtain the harmonic frequency weighting coefficients ak(t) and bk(t).
In some implementations, the adaptive algorithm is a least mean square (LMS) algorithm.
In certain implementations, extracting the phase and shape information for each of the N harmonics includes smoothing out in time the harmonic frequency weighting coefficients ak(t) and bk(t), thereby obtaining smoothed harmonic frequency weighting coefficients ãk(t) and {tilde over (b)}k(t), for each of the N harmonics over the RPM range.
In some implementations, the smoothing is performed by low pass filtering the harmonic frequency weighting coefficients ak(t) and bk(t).
In certain implementations, the low pass filtering is performed with a finite impulse response (FIR) filter.
In some implementations, the finite impulse response (FIR) filter is based on a window (e.g., Hamming, Hanning, Blackman, Blackman-Harris, Kaiser, Chebyshev, etc.), where a length of the window determines the degree by which the harmonic frequency coefficients will be smoothed.
In certain implementations, the low pass filtering is performed with an infinite impulse response (IIR) filter.
In some implementations, extraction of the phase and shape information includes: computing a shape, {tilde over (c)}k(t), and a phase, {tilde over (p)}k(t), of each of the N harmonics from the smoothed harmonic frequency weighting coefficients ãk(t) and {tilde over (b)}k(t); assuming that an RPM-dependent shape, {tilde over (c)}k(RPM), and an RPM-dependent phase, {tilde over (p)}k(RPM), of each of the N harmonics over the RPM range is represented by the shape, {tilde over (c)}k(t), and the phase, {tilde over (p)}k(t), computed in time, such that {tilde over (c)}k(t)={tilde over (c)}k(RPM) and {tilde over (p)}k(t)={tilde over (p)}k(RPM); and defining each of the N harmonics as a pair of vectors [Ck] and [Pk] that contain samples of the RPM-dependent shape, {tilde over (c)}k(RPM), and of the RPM-dependent phase, {tilde over (p)}k(RPM), respectively, over a subset of RPM values.
In certain implementations, the shape of each of the N harmonics is computed as:




  
    
      
        
          c
          ~
        
        k
      
      ⁡
      
        (
        t
        )
      
    
    =
    
      
        
          
            
              a
              ~
            
            k
            2
          
          ⁡
          
            (
            t
            )
          
        
        +
        
          
            
              b
              ~
            
            k
            2
          
          ⁡
          
            (
            t
            )
          
        
      
    
  
  ;




and the phase of each of the N harmonics is computed as:





  
    
      
        p
        ~
      
      k
    
    ⁡
    
      (
      t
      )
    
  
  =
  
    unwrap
    ⁢
    
      
        {
        
          atan
          ⁡
          
            (
            
              
                
                  
                    
                      b
                      ~
                    
                    k
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
                
                  
                    
                      a
                      ~
                    
                    k
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
              
              -
              
                
                  π
                  ·
                  max
                
                ⁢
                
                  {
                  
                    
                      sign
                      ⁡
                      
                        (
                        
                          
                            a
                            k
                          
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                        )
                      
                    
                    ,
                    0
                  
                  }
                
              
            
            )
          
        
        }
      
      .
    
  




In some implementations, the method further includes operating the EHE system to transduce the engine harmonic enhancement signal, he(t), to acoustic energy.
In certain implementations, configuring the engine harmonic enhancement system to utilize the extracted phase and shape information to generate the engine harmonic enhancement signal, he(t), includes configuring the engine harmonic enhancement system to: receive an RPM signal representing a current RPM of a vehicle engine; based on the RPM signal, determine an instant magnitude, ck(t), for each of the N harmonics by interpolating the instant magnitude, ck(t), from a magnitude table [Ck] wherein the magnitude table [Ck] includes samples of an RPM-dependent shape, {tilde over (c)}k(RPM), for each of the N harmonics over a subset of RPM values; based on the RPM signal, determine an instant phase, pk(t), for each of the N harmonics by interpolating the instant phase, pk(t), from a phase table [Pk], wherein the phase table [Pk] includes samples of an RPM-dependent phase {tilde over (p)}k(RPM), for each of the N harmonics over the subset of RPM values; smooth the instant magnitude, ck(t), to generate a smoothed harmonic magnitude, ĉk(t), for each of the N harmonics; smooth the instant phase, pk(t), to generate a smoothed harmonic phase, {circumflex over (p)}k(t), for each of the N harmonics; generate a magnitude perturbation signal, dck(t), for each of the N harmonics; generate a phase perturbation signal, dpk(t), for each of the N harmonics; and generate an individual harmonic signal, hk(t), for each of the N harmonics according to: hk(t)=(ĉk(t)+dck(t))·sin(2πflkt+{circumflex over (p)}k(t)+dpk(t)); and sum the individual harmonic signals, hk(t), and thereby generate the engine harmonic enhancement signal, he(t).
In some implementations, N is an integer from 40 to 80.
In certain implementations, the RPM range is 600 RPM to 7500 RPM.
In another aspect, a non-transitory computer-readable medium bears instructions to cause a processor to perform a harmonic decomposition on a target engine sound, thereby to extract each of N harmonics of the target engine sound over an RPM range; and out of the extracted harmonics, extract phase and shape information for each of the N harmonics over the RPM range for reproducing the target engine sound.
Implementations may include one of the above and/or below features, or any combination thereof.
In yet another aspect, a method includes: in a vehicle engine harmonic enhancement system, receiving an RPM signal representing a current RPM of a vehicle engine; based on the RPM signal, determining an instant magnitude, ck(t), for each of N harmonics of a target engine sound by interpolating the instant magnitude, ck(t), from a vector [Ck], wherein the vector [Ck] includes samples of an RPM-dependent shape, {tilde over (c)}k(RPM), for each of the N harmonics over a subset of RPM values; based on the RPM signal, determining an instant phase, pk(t), for each of the N harmonics by interpolating the instant phase, pk(t), from a vector [Pk], wherein the vector [Pk] includes samples of an RPM-dependent phase {tilde over (p)}k(RPM), for each of the N harmonics over the subset of RPM values; smoothing the instant magnitude, ck(t), to generate a smoothed harmonic magnitude, ĉk(t), for each of the N harmonics; smoothing the instant phase, pk(t), to generate a smoothed harmonic phase, {circumflex over (p)}k(t), for each of the N harmonics; generating a magnitude perturbation signal, dck(t), for each of the N harmonics; generating a phase perturbation signal, dpk(t), for each of the N harmonics; and generating an individual harmonic signal, hk(t), for each of the N harmonics according to: hk(t)=(ĉk(t)+dck(t))·sin(2πflkt+{circumflex over (p)}k(t)+dpk(t)); summing the individual harmonic signals, hk(t), and thereby generating an engine harmonic enhancement signal, he(t) to be added on top of a baseline engine sound; and transducing the engine harmonic enhancement signal, he(t), to acoustic energy.
Implementations may include one of the above and/or below features, or any combination thereof.
In some implementations, generating the magnitude perturbation signal includes generating the magnitude perturbation signal, dck(t), for each of the N harmonics according to:




  
    
      d
      ⁢
      
        
      
      ⁢
      
        
          c
          k
        
        ⁡
        
          (
          t
          )
        
      
    
    =
    
      
        
          ncf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              c
              k
            
          
        
        
          σ
          
            ncf
            k
          
        
      
    
  
  ,




where σΔck, is the standard deviation of fluctuations in the harmonic magnitude of each of the N harmonics of the target engine sound; ncfk(t) is filtered random noise used for generating the magnitude perturbation signal; and σncfk is the standard deviation for the filtered random noise used for generating the magnitude perturbation signal.

In certain implementations, generating the phase perturbation signal includes generating the phase perturbation signal, dpk(t), for each of the N harmonics according to:




  
    
      
        dp
        k
      
      ⁡
      
        (
        t
        )
      
    
    =
    
      
        
          npf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              k
            
          
        
        
          σ
          
            npf
            k
          
        
      
    
  
  ,




where σΔpk, is the standard deviation of fluctuations in the harmonic phase of each of the N harmonics of the target engine sound; npfk(t) is filtered random noise used for generating the phase perturbation signal; and σnpfk is the standard deviation for the filtered random noise used for generating the phase perturbation signal.

In another aspect, an engine harmonic enhancement system, includes: a digital signal processor; and a memory coupled to the digital signal processor and containing instructions that, when executed, cause the digital signal processor to: receive an RPM signal representing a current RPM of a vehicle engine; based on the RPM signal, determine an instant magnitude, ck(t), for each of N harmonics of a target engine sound by interpolating the instant magnitude, ck(t), from a vector [Ck], wherein the vector [Ck] includes samples of an RPM-dependent shape, {tilde over (c)}k(RPM), for each of the N harmonics over a subset of RPM values; based on the RPM signal, determine an instant phase, pk(t), for each of the N harmonics by interpolating the instant phase, pk(t), from a vector [Pk], wherein the vector [Pk] includes samples of an RPM-dependent phase {tilde over (p)}k(RPM), for Each of the N harmonics over the subset of RPM values; smooth the instant magnitude, ck(t), to generate a smoothed harmonic magnitude, ĉk(t), for each of the N harmonics; smooth the instant phase, pk(t), to generate a smoothed harmonic phase, {circumflex over (p)}k(t), for each of the N harmonics; generate a magnitude perturbation signal, dck(t), for each of the N harmonics; generate a phase perturbation signal, dpk(t), for each of the N harmonics; and generate an individual harmonic signal, hk(t), for each of the N harmonics according to: hk(t)=(ĉk(t)+dck(t))·sin(2πflkt+{circumflex over (p)}k(t)+dpk(t)); sum the individual harmonic signals, hk(t), and thereby generating an engine harmonic enhancement signal, he(t) to be added on top of a baseline engine sound; and transduce the engine harmonic enhancement signal, he(t), to acoustic energy.
Implementations may include one of the above features, or any combination thereof.



BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a flow diagram of an example process for extracting parameters that define the character of a target engine sound and which can be used for reproducing the character of a target engine sound.
FIG. 2 is a block diagram of an adaptive linear combiner.
FIG. 3 is a flow diagram of an example process for extracting phase and shape information for preserving the character of the target engine sound.
FIG. 4 is a block diagram of a vehicle including a vehicle engine sound enhancement system.
FIG. 5 is a flow diagram of an example process, which may be employed by the vehicle engine sound enhancement system of FIG. 4, for reproducing the desired, target engine sound in a subject automobile.



DETAILED DESCRIPTION
Though the elements of several figures may be shown and described as discrete elements in a block diagram and may be referred to as “circuitry”, unless otherwise indicated, the elements may be implemented as one of, or a combination of, analog circuitry, digital circuitry, or one or more microprocessors executing software instructions. The software instructions may include digital signal processing (DSP) instructions. Operations may be performed by analog circuitry or by a microprocessor executing software that performs the mathematical or logical equivalent to the analog operation. Unless otherwise indicated, signal lines may be implemented as discrete analog or digital signal lines, as a single discrete digital signal line with appropriate signal processing to process separate streams of audio signals, or as elements of a wireless communication system. Some of the processes may be described in block diagrams. The activities that are performed in each block may be performed by one element or by a plurality of elements, and may be separated in time. The elements that perform the activities of a block may be physically separated. One element may perform the activities of more than one block. Unless otherwise indicated, audio signals or video signals or both may be encoded and transmitted in either digital or analog form; conventional digital-to-analog or analog-to-digital converters may be omitted from the figures.
The degree by which the character of a target sound can be reproduced in a vehicle depends on the vehicle's engine baseline sound. FIG. 1 illustrates an exemplary process that can be particularly useful for applications where large enhancements, relative to the vehicle baseline, are desired. The psychoacoustic parameters that define the character of a target engine sound can be extracted directly from the target engine sound.
Referring to FIG. 1, an exemplary process 100 of reproducing the character of a particular engine sound (a target engine sound). The process 100 includes performing a harmonic decomposition on the target engine sound (120) to extract the harmonics of the target engine sound, and, out the extracted harmonics, extracting sufficient phase and shape information (e.g., for each of the extracted harmonics) to reproduce the character of target engine sound (130). Then, an enhancement signal can be generated (140) which can be added on top of a baseline engine sound to reproduce the target engine sound.
The degree by which the target sound is reproduced can depend on the baseline sound. If the baseline sound harmonics are louder than all the target sound harmonics, there will be no change in the obtained sound. However, the level of enhancement can be increased to a degree that will dominate the baseline, in which case the obtained sound can get very close to the target engine sound.
Engine Sound Harmonic Decomposition
The engine sound of a vehicle, as perceived in the cabin of the vehicle, can be assumed to be the sum of engine speed harmonics, broadband engine noise, wind noise, tire noise, and noise generated by other sources. To simplify, the engine sound can be represented by the sum of engine speed harmonics plus broadband noise. Considering only the harmonics of the engine sound, each harmonic amplitude and phase can be assumed to vary in time according to the following formula:

hk(t)=ck(t)·sin(ωk(t)·t+φk(t))

where




  
    
      ω
      k
    
    =
    
      
        2
        ⁢
        π
        ⁢
        
          
        
        ⁢
        
          f
          k
        
      
      =
      
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            l
            k
          
          ⁢
          
            f
            1
          
        
        =
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            l
            k
          
          ⁢
          
            RPM
            60
          
        
      
    
  
  ;




f1 is the fundamental frequency; lk is the harmonic order that usually is a multiple of 0.5, but, in general, can be a positive non-zero real number; RPM is the engine speed; and fk is the frequency, in Hz, of harmonic k.

The target engine sound can be written as the sum of all contributing harmonics:

x(t)=Σk=1Nck(t)·sin(ωk(t)·t+φk(t))

FIG. 2 illustrates an implementation of an adaptive linear combiner 200 that can be utilized to extract the harmonics. The adaptive linear combiner 200 can be used to identify the engine harmonics and their evolution in time. The adaptive linear combiner 200 can allow for the continuous identification under any engine state and engine mode of operation.
The adaptive linear combiner 200 accepts an RPM signal 210 and a microphone signal 212 (i.e., the target engine sound signal, s(t), e.g., as measured by one or more microphones in a passenger compartment of the vehicle which produces the target engine sound) and processes the inputs 210, 212 to generate an engine noise-free signal 224 (i.e., the engine harmonics hest(t)) and the broadband engine noise, bn(t) 214.
The adaptive linear combiner 200 includes N harmonic frequency generators 216, N sets of harmonic frequency weighting coefficients 218, N summers 219, a harmonics summer 220, and an adaptation control module 222.
Each harmonic can be expressed as follows:




  
    
      
        
          
            h
            k
          
          ⁡
          
            (
            t
            )
          
        
        =
        
        ⁢
        
          
            
              c
              k
            
            ⁡
            
              (
              t
              )
            
          
          ·
          
            sin
            ⁡
            
              (
              
                
                  
                    
                      ω
                      k
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                +
                
                  
                    φ
                    k
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
              
              )
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              
                c
                k
              
              ⁡
              
                (
                t
                )
              
            
            ·
            
              sin
              ⁡
              
                (
                
                  
                    
                      ω
                      k
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
            ·
            
              cos
              ⁡
              
                (
                
                  
                    φ
                    k
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
                )
              
            
          
          +
          
            
              
                c
                k
              
              ⁡
              
                (
                t
                )
              
            
            ·
            
              cos
              ⁡
              
                (
                
                  
                    
                      ω
                      k
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
            ·
            
              sin
              ⁡
              
                (
                
                  
                    φ
                    k
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
                )
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              
                a
                k
              
              ⁡
              
                (
                t
                )
              
            
            ·
            
              sin
              ⁡
              
                (
                
                  
                    
                      ω
                      k
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
          
          +
          
            
              
                b
                k
              
              ⁡
              
                (
                t
                )
              
            
            ·
            
              cos
              ⁡
              
                (
                
                  
                    
                      ω
                      k
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
          
        
      
    
  




Where: ak(t)=ck(t)·cos(φk(t)) and bk(t)=ck(t)·sin(φk(t))
Thus, to generate the unweighted sine and cosine components of each of N harmonics, each of the N harmonic frequency generators 216 includes both a sine component generator 216a and a cosine component generator 216b. When the harmonic frequency generators 216 receive the RPM signal 210, the RPM signal 210 is converted to a frequency and used by the sine and cosine component generators 216a, 216b to generate sine and cosine components for the frequency. For example, based on the RPM signal 210, the frequencies of both the sine component 216a and the cosine component 216b of the first harmonic generator 216 are set to the fundamental angular frequency (i.e., ω1=2πf1) of the target engine sound, the frequencies of sine and cosine components of the second harmonic frequency generator are set to twice the fundamental angular frequency (i.e., 2·ω1=ω2=4πf1)) of the target engine sound, the frequencies of the sine and cosine components of the third harmonic frequency generator are set to three times the fundamental angular frequency (i.e., 3·ω1=ω3=6πf1) of the target engine sound, and so on. This scheme includes full order, half orders, and can also include fractional orders. In some examples the sine and cosine component generators 216a, 216b include look-up tables for generating the sine and cosine components.
The N sine components 216a and the N cosine components 216b generated by each of the harmonic frequency generators 216 are then passed to the N sets of harmonic frequency weighting coefficients 218. Each set of harmonic frequency weighting coefficients 218 includes a first weighting coefficient ak(t) 218a for weighting the sine component 216a of the generated harmonic frequency and a second weighting coefficient bk(t) 218b for weighting the cosine component 216b of the generated harmonic frequency.
After the generated harmonic frequency components 216a, 216b for each of the harmonic frequency generators 216 are weighted by their corresponding harmonic frequency weighting coefficients 218 the pairs of weighted harmonic frequency components are summed by respective summers 219 to produce N harmonics h1(t)−hN(t), which are then summed at harmonics summer 220 to produce an estimate of the engine noise 224 (i.e., the estimated harmonic engine sound, hest(t)). The estimate of the engine noise 224 is combined with the microphone signal 212, resulting in the broadband noise, bn(t).
In operation, the adaptive linear combiner 200 adaptively determines the estimate of the engine noise 224. That is, the adaptive linear combiner 200 adapts the N sets of harmonic frequency weighting coefficients 218 to converge to the engine noise-free signal 214 (i.e., the broadband noise, bn(t)). Thus, the engine noise-free signal 214 can also be viewed as an error signal for an adaptive filter. This error signal is passed to an adaptation control module 222 which uses the error signal to determine a coefficient update for the N sets of harmonic frequency weighting coefficients 218.
In some examples, the adaptation control module 222 determines the coefficient update using a least mean squares (LMS) approach. The following exemplifies the use of the LMS algorithm.
Using the following notation: X(t) is the input vector at moment t; W(t) is the coefficients vector; hest(t) is the estimated harmonic engine sound is the output of the linear combiner; s(t) the microphone signal is the desired signal; and the error in estimation is the broadband noise bn(t).




  
    
      X
      ⁡
      
        (
        t
        )
      
    
    =
    
      [
      
        
          
            
              sin
              ⁡
              
                (
                
                  
                    
                      ω
                      1
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
          
        
        
          
            
              cos
              ⁡
              
                (
                
                  
                    
                      ω
                      1
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
          
        
        
          
            ⋮
          
        
        
          
            
              sin
              ⁡
              
                (
                
                  
                    
                      ω
                      N
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
          
        
        
          
            
              cos
              ⁡
              
                (
                
                  
                    
                      ω
                      N
                    
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                  ·
                  t
                
                )
              
            
          
        
      
      ]
    
  
  ,
  
    

  
  ⁢
  
    
      W
      ⁡
      
        (
        t
        )
      
    
    =
    
      [
      
        
          
            
              
                a
                1
              
              ⁡
              
                (
                t
                )
              
            
          
        
        
          
            
              
                b
                1
              
              ⁡
              
                (
                t
                )
              
            
          
        
        
          
            ⋮
          
        
        
          
            
              
                a
                N
              
              ⁡
              
                (
                t
                )
              
            
          
        
        
          
            
              
                b
                N
              
              ⁡
              
                (
                t
                )
              
            
          
        
      
      ]
    
  




The update equations for the adaptive linear combiner are:

hest(t)=X(t)T·W(t)

bn(t)=s(t)−y(t)

W(t+1)=W(t)+μ·bn(t)·X(t)

The step size of the adaptive algorithm, μ, can be set to a value that extracts the harmonics only, separating these from the broadband noise. The step size μ can be selected based on subjective listening tests, by listening to the broadband noise, bn(t)=s(t)−hest(t). If one cannot hear any harmonic content, the identification of the harmonics is considered successful. A visual identification of a spectrogram can validate the subjective findings.
Alternatively, a good identification of the harmonics can be recognized by comparing loudness spectrograms of the extracted harmonics and of the residual broadband noise for different step sizes. For the residual broadband noise, a well spread loudness, without peaks in the frequency domain, over the duration of the sound recording is identified. The lowest step size that shows this property is the one that can be considered for analysis. When analyzing the extracted harmonics, the loudness spectrogram that shows all the harmonic details, for the lowest step size is identified. This can be recognized by the fact that higher step sizes do not add more definition to the loudness spectrogram. This is the step size that can be considered in the analysis. Finally, comparing the step sizes identified via the previous tests, the lowest that satisfies both tests is selected. This is the step size that will be used in the identification process. The rate of change for lower versus higher order harmonics is different; therefore choosing a different step size per harmonic is helpful.
Other complex order analysis methods can also be used to identify the harmonics. Methods based on fast Fourier transform (FFT) constitute a good example. In other examples, the adaptation linear combiner 200 determines the coefficient update using a recursive least squares (RLS) approach.
Extraction of Phase and Shape Information
Start with a recording of the engine sound that is desired to be reproduced. Preferably, the recording is taken under wide open throttle (WOT), starting from the lowest RPM possible up to red line, preferably in the 3rd gear. Running a harmonic decomposition, as described in the previous section, one can obtain the harmonic frequency weighting coefficients ak(t) and bk(t) for each of N harmonics over an RPM range.
After the harmonic frequency weighting coefficients ak(t) and bk(t) are obtained via the harmonic decomposition process described in the previous section, the next step is to extract the minimum phase and shape information (130, FIG. 1) needed to preserve the character of the target engine sound. Referring to FIG. 3, the first step in extraction of the phase and shape information (130) is to smooth out in time the harmonic frequency weighting coefficients ak(t) and bk(t), obtaining ãk(t) and {tilde over (b)}k(t) (310), for each of the N harmonics over the RPM range.
The smoothing is performed by low pass filtering ak(t) and bk(t). The LP filter can be a finite impulse response (FIR) filter based on a Hamming window, where the length of the window determines the degree by which the coefficients will be smoothed. The shape (magnitude) of each of the N harmonics can then be computed (312) as:




  
    
      
        c
        ~
      
      k
    
    ⁡
    
      (
      t
      )
    
  
  =
  
    
      
        
          
            a
            ~
          
          k
          2
        
        ⁡
        
          (
          t
          )
        
      
      +
      
        
          
            b
            ~
          
          k
          2
        
        ⁡
        
          (
          t
          )
        
      
    
  




and the phase is computed (312) as:




  
    
      
        p
        ~
      
      k
    
    ⁡
    
      (
      t
      )
    
  
  =
  
    unwrap
    ⁢
    
      {
      
        atan
        ⁡
        
          (
          
            
              
                
                  
                    b
                    ~
                  
                  k
                
                ⁡
                
                  (
                  t
                  )
                
              
              
                
                  
                    a
                    ~
                  
                  k
                
                ⁡
                
                  (
                  t
                  )
                
              
            
            -
            
              
                π
                ·
                max
              
              ⁢
              
                {
                
                  
                    sign
                    ⁡
                    
                      (
                      
                        
                          
                            a
                            ~
                          
                          k
                        
                        ⁡
                        
                          (
                          t
                          )
                        
                      
                      )
                    
                  
                  ,
                  0
                
                }
              
            
          
          )
        
      
      }
    
  




For the FIR filter, a window length between of 5 ms to 20 ms is suitable for this smoothing. Furthermore, while an FIR LP filter with a Hamming window has been described, any type of window can be used, e.g., Hanning, Blackman, Blackman-Harris, Kaiser, Chebyshev, etc. Furthermore, any type of LP filter can be used including, for example, an infinite impulse response (IIR) filter. For the IIR filter there will be no windowing used.
The engine sound can be approximated by the sum of harmonics expressed as:

hk(t)={tilde over (c)}k(t)·sin(2πflkt+{tilde over (p)}k(t))

The computed phase {tilde over (p)}k(t) can also be looked at as a deviation in frequency from the actual harmonic frequency. This frequency deviation is:




  
    Δ
    ⁢
    
      
    
    ⁢
    
      
        f
        k
      
      ⁡
      
        (
        t
        )
      
    
  
  =
  
    
      
        
          p
          ~
        
        k
      
      ⁡
      
        (
        t
        )
      
    
    
      2
      ⁢
      π
      ⁢
      
        
      
      ⁢
      
        l
        k
      
      ⁢
      t
    
  




Then, it will be assumed that the shape and phase of each harmonic over the RPM range is represented by the shape and phase obtained in the measurement in time (314):

{tilde over (c)}k(t)={tilde over (c)}k(RPM)

{tilde over (p)}k(t)={tilde over (p)}k(RPM)

where, RPM, is the engine speed over the sweep. The two signals are defined in a very high resolution, with samples taken at the sampling rate.
This can be too much information to store; however, the amount of information can be reduced, by defining each harmonic as two vectors that contain samples of {tilde over (c)}k(RPM) and {tilde over (p)}k(RPM). First, the samples for which the RPM values increase monotonically are selected. Then, those signals are decimated to obtain vectors [Ck] and [Pk] (also known as the magnitude and phase tables), which define the magnitude (shape) and phase of each of the N harmonics over a subset of RPM values (316). The points can be spaced equidistantly over the RPM range, or selected, in such a way, to maximize the retained information.
An alternative approach is to further smooth {tilde over (c)}k(t) and {tilde over (p)}k(t), obtaining k(t) and k(t), which can then be used to compute the magnitude and phase tables, [Ck] and [Pk].
Generation of Enhancement Signal
The phase and shape information about the target engine sound that was derived via the decomposition (120, FIG. 1) and extraction (130, FIG. 1) steps discussed in the previous sections can then be utilized to reproduce the desired, target engine sound in a subject vehicle. As used herein, the term “subject vehicle” refers to a vehicle in which an engine harmonic enhancement is provided to reproduce the target engine sound. The subject automobile may, for example, be a particular automobile model that is different from another automobile model which produced the original, recorded target engine sound.
FIG. 4 is a block diagram of a subject vehicle including a vehicle sound enhancement system 400. An engine harmonic enhancement (EHE) processor 412 (e.g., a digital signal processor) is coupled to a summer 414. Also coupled to the summer 414 by an entertainment audio equalizer and spatial processor 416 is an entertainment audio signal source 418. The summer 414 is coupled through a multi-channel amplifier 420 to a number of loudspeakers 422-1-422-4 positioned about the vehicle cabin, and in some implementations, may be positioned, for example, as loudspeaker 424 to radiate acoustic energy to the exterior of the vehicle. The couplings between the EHE processor 412, the entertainment audio equalizer and spatial processor 416, the summer 414, and the amplifier 420 may be multi-channel, as indicated by the multiple lines. Signal lines may be implemented as discrete analog or digital signal lines, as a single discrete digital signal line with appropriate signal processing to process separate streams of audio signals, or as elements of a wireless communication system.
In operation, the entertainment audio source 418 and entertainment audio equalizer and spatial processor 416 operate conventionally, to provide equalized and spatially processed audio entertainment to occupants of the vehicle cabin. In some implementations, the entertainment audio signal source 418 can include announcement audio signals, for navigation, warning signals, and the like. The EHE processor 412 receives an RPM signal 426 representative of a number of revolutions per minute of an engine 427 of the subject vehicle as input. The EHE processor 412 uses the input signal 426 to form an engine harmonic enhancement signal, he(t), which is provided to speakers 422-1-422-4 and 424 for presentation to an automobile operator. The engine harmonic enhancement signal, he(t), and the processed entertainment audio signals are summed at summer 414, amplified by amplifier 420, and transduced to acoustic energy by the loudspeakers 422-1-422-4, and 424
FIG. 5 illustrates an example process 500 that can be employed (e.g., via the EHE processor 412 of FIG. 4) to reproduce the desired, target engine sound in a subject automobile. The RPM signal 426 is received (510) by the EHE processor 412 of FIG. 4. An angular frequency




  (
  
    
      ω
      k
    
    =
    
      
        2
        ⁢
        π
        ⁢
        
          
        
        ⁢
        
          f
          k
        
      
      =
      
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            l
            k
          
          ⁢
          
            f
            1
          
        
        =
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            l
            k
          
          ⁢
          
            RPM
            60
          
        
      
    
  
  )




is calculated for each of N harmonics based on based on the RPM signal 426 (512).

Based on the current RPM, the instant magnitude, ck(t), can be determined for each harmonic by interpolating the value from the vector [Ck] (514), and the instant phase, pk(t), can be determined for each harmonic by interpolating the value from the vector [Pk] (516); the vectors [Ck] and [Pk] having been previously determined (e.g., by the vehicle and/or sound system manufacturer) via the extraction process.
Because the RPM is generally updated only once every 10 to 30 ms, regardless of the number of points used in the magnitude and phase tables, the interpolated signal will show abrupt changes when the RPM changes. For this reason, the interpolation of the instant magnitude, ck(t), (514) is followed by a first order magnitude smoothing step (518) to eliminate abrupt changes to the instant magnitude, ck(t), when the RPM value is updated. Similarly, the interpolation of the instant phase, pk(t), (516) is followed by a first order phase smoothing step (520) to eliminate abrupt changes to the instant phase, pk(t), when the RPM value is updated. It is beneficial to smooth both the phase and the magnitude of each harmonic, but the phase impacts more the quality of the generated signal. The smoothed generated harmonic magnitude will be denoted by ĉk(t), and the smoothed generated harmonic phase by {circumflex over (p)}k(t).
Because of the smoothing done in the table computation step (the extraction process (130)), the naturally occurring variations in the magnitude and phase are filtered out. Because the actual fluctuations in the magnitude of each harmonic are not available in the harmonic generation stage, fluctuations can be created that have similar statistical properties as those measured in the extraction stage.
To estimate the fluctuations that where lost in the approximation of each harmonic magnitude, {tilde over (c)}k(t) is compared to k(t) to compute Δck according to:

Δck(t)={tilde over (c)}k(t)−k(t)

The standard deviation of Δck, σΔck, can be determined as part of the extraction process or at some other point prior to the generation process. The standard deviation, σΔck, will be used as a measure of the magnitude fluctuation for each of the N harmonics, in the generation process.
Starting with the generation of random noise, nck(t), with a normal distribution, that is low pass filtered with a cutoff frequency around 25 Hz to obtain ncfk(t). The noise signals need to be different for each harmonic. The standard deviation for the filtered noise can be estimated based on the generation method, and will be denoted here as σnfck.
The fluctuations in the magnitude are obtained by generating the following magnitude perturbation signal for each of the N harmonics (522):




  
    
      d
      ⁢
      
        
      
      ⁢
      
        
          c
          k
        
        ⁡
        
          (
          t
          )
        
      
    
    =
    
      
        
          ncf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              c
              k
            
          
        
        
          σ
          
            ncf
            k
          
        
      
    
  
  ,




where ncfk(t) is the filtered noise used in the magnitude generation, and where dck(t) is the distortion (perturbation) signal that creates the fluctuations in the generated harmonic magnitude.
Similarly a fluctuation signal can be generated for the phase (524). To estimate the fluctuations that where lost in the approximation of each harmonic phase, {tilde over (p)}k(t) is compared to k(t) to compute Δpk(t) according to:

Δpk(t)={tilde over (p)}k(t)−k(t),

The standard deviation of Δpk, σΔpk, can be determined as part of the extraction process or at some other point prior to the generation process. The standard deviation, σΔpk, will be used as a measure of the magnitude fluctuation for each of the N harmonics, in the generation process.
Starting with the generation of random noise, npk(t), with a normal distribution, that is low pass filtered with a cutoff frequency around 25 Hz to obtain npfk(t). The noise signals need to be different for each harmonic. The standard deviation for the filtered noise can be estimated based on the generation method, and will be denoted here as σnpfk.
The fluctuations in the phase are obtained by generating the following phase perturbation signal for each of the N harmonics (524):




  
    
      
        dp
        k
      
      ⁡
      
        (
        t
        )
      
    
    =
    
      
        
          npf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              k
            
          
        
        
          σ
          
            npf
            k
          
        
      
    
  
  ,




where npfk(t) is the filtered noise used in the phase generation, and where dpk(t) is the distortion (perturbation) signal that creates the fluctuations in the generated harmonic phase.
An individual harmonic signal, hk(t), is generated (526) for each of the N harmonics according to:

hk(t)=(ĉk(t)+dck(t))·sin(2πflkt+{circumflex over (p)}k(t)+dpk(t))

Notably, the phase information for each of the N harmonics is included in the generation of the individual harmonic signals (526). The individual harmonic signals, hk(t), are added together to generate an engine harmonic enhancement signal, he(t)=Σk=1Nhk(t), that will be played via the sound system of the subject vehicle to complement the original engine sound of the subject vehicle (528).
In some implementations, the EHE processor 412 can be configured to determine an amount of gain or delay to be applied to the engine harmonic enhancement signal, he(t). For example, the EHE processor 412 may include an EHE gain and delay determiner which determines an amount of gain and delay to be applied to the engine harmonic enhancement signal, he(t), by an EHE overall enhancement gain. The EHE gain and delay determiner may apply a gain function (also referred to as “mapping function” or “mapping”) which includes as variables the engine load, the change in engine load, the RPM, and the rate of change in RPM to determine the EHE gain (as described in U.S. patent application Ser. No. 12/716,887). Additionally, the gain function applied by the EHE gain and delay determiner may use as variables values of other parameters, such as the gear in which the vehicle is operating; the transmission ratio or transmission ratio interval of a continuously variable transmission (CVT); and an operational mode of the engine. The EHE gain and delay determiner may smooth the gain values so that the sound variation is natural, and undistorted, similar to the sound variation in time of a mechanical system.
Alternatively or additionally, the EHE processor 412 can include a sound stage processor (as described in U.S. patent application Ser. No. 13/004,630) that processes the engine harmonic enhancement signal, he(t), through a separate equalization filter for each loudspeaker 422-1-422-4 and 424 of FIG. 4. The engine harmonic enhancement signal, he(t), can be monophonic, indicating that the same signal is provided to all loudspeakers 422-1-422-4 and 424 (FIG. 4), or may be multichannel, for example stereophonic.
Systems that implement the techniques described above can be implemented in software, firmware, in digital electronic circuitry, or in computer hardware, or in combinations of them. The systems can include a computer program product tangibly embodied in a machine-readable storage device for execution by a programmable processor, and method steps can be performed by programmable processor executing a program of instructions to perform functions by operating on input data and generating output. The systems can be implemented in one or more computer programs that are executable on a programmable system including at least one programmable processor coupled to receive data and instructions from, and to transmit data and instructions to, a data storage system, at least one input device, and at least one output device. Each computer program can be implemented in a high-level procedural or object-oriented programming language, or in assembly or machine language if desired; and in any case, the language can be a compiled or interpreted language. Suitable processors include, by way of example, both general and special purpose microprocessors. Generally, a processor will receive instructions and data from read-only memory and/or a random access memory. Generally, a computer will include one or more mass storage devices for storing data files; such devices include magnetic disks, such as internal hard disks and removable disks; magneto-optical disks; and optical disks. Storage devices suitable for tangibly embodying computer program instructions and data include all forms of non-volatile memory, including by way of example semiconductor memory devices, such as EPROM, EEPROM, and flash memory devices; magnetic disks such as internal hard disks and removable disks; magneto-optical disks; and CD-ROM disks. Any of the foregoing can be supplemented by or, or incorporated in, ASICs (application-specific integrated circuits).
A number of implementations have been described. Nevertheless, it will be understood that additional modifications may be made without departing from the spirit and scope of the inventive concepts described herein, and, accordingly, other implementations are within the scope of the following claims.

## Claims

1. A method comprising:
in a vehicle engine harmonic enhancement system,
receiving an RPM signal representing a current RPM of a vehicle engine;
based on the RPM signal, determining an instant magnitude, ck(t), for each of N harmonics of a target engine sound by interpolating the instant magnitude, ck(t), from a vector [Ck], wherein the vector [Ck] comprises samples of an RPM-dependent shape, {tilde over (c)}k(RPM), for each of the N harmonics over a subset of RPM values;
based on the RPM signal, determining an instant phase, pk(t), for each of the N harmonics by interpolating the instant phase, pk(t), from a vector [Pk], wherein the vector [Pk] comprises samples of an RPM-dependent phase {tilde over (p)}k(RPM), for each of the N harmonics over the subset of RPM values;
smoothing the instant magnitude, ck(t), to generate a smoothed harmonic magnitude, ĉk(t), for each of the N harmonics;
smoothing the instant phase, pk(t), to generate a smoothed harmonic phase, {circumflex over (p)}k(t), for each of the N harmonics;
generating a magnitude perturbation signal, dck(t), for each of the N harmonics;
generating a phase perturbation signal, dpk(t), for each of the N harmonics; and
generating an individual harmonic signal, hk(t), for each of the N harmonics according to: hk(t)=(ĉk(t)+dck(t))·sin(2πflkt+{circumflex over (p)}k(t)+dpk(t));
summing the individual harmonic signals, hk(t), and thereby generating an engine harmonic enhancement signal, he(t) to be added on top of a baseline engine sound; and
transducing the engine harmonic enhancement signal, he(t), to acoustic energy.

2. The method of claim 1, wherein generating the magnitude perturbation signal comprises generating the magnitude perturbation signal, dck(t), for each of the N harmonics according to:




  
    
      d
      ⁢
      
        
      
      ⁢
      
        
          c
          k
        
        ⁡
        
          (
          t
          )
        
      
    
    =
    
      
        
          ncf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              c
              k
            
          
        
        
          σ
          
            ncf
            k
          
        
      
    
  
  ,



where σΔck, is the standard deviation of fluctuations in the harmonic magnitude of each of the N harmonics of the target engine sound; ncfk(t) is filtered random noise used for generating the magnitude perturbation signal; and σncfk is the standard deviation for the filtered random noise used for generating the magnitude perturbation signal.

3. The method of claim 1, wherein generating the phase perturbation signal comprises generating the phase perturbation signal, dpk(t), for each of the N harmonics according to:




  
    
      
        dp
        k
      
      ⁡
      
        (
        t
        )
      
    
    =
    
      
        
          npf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              k
            
          
        
        
          σ
          
            npf
            k
          
        
      
    
  
  ,



where σΔpk, is the standard deviation of fluctuations in the harmonic phase of each of the N harmonics of the target engine sound; npfk(t) is filtered random noise used for generating the phase perturbation signal; and σnpfk is the standard deviation for the filtered random noise used for generating the phase perturbation signal.

4. An engine harmonic enhancement system, comprising:
a digital signal processor; and
a memory coupled to the digital signal processor and containing instructions that, when executed, cause the digital signal processor to:
receive an RPM signal representing a current RPM of a vehicle engine;
based on the RPM signal, determine an instant magnitude, ck(t), for each of N harmonics of a target engine sound by interpolating the instant magnitude, ck(t), from a vector [Ck], wherein the vector [Ck] comprises samples of an RPM-dependent shape, {tilde over (c)}k(RPM), for each of the N harmonics over a subset of RPM values;
based on the RPM signal, determine an instant phase, pk(t), for each of the N harmonics by interpolating the instant phase, pk(t), from a vector [Pk], wherein the vector [Pk] comprises samples of an RPM-dependent phase {tilde over (p)}k(RPM), for each of the N harmonics over the subset of RPM values;
smooth the instant magnitude, ck(t), to generate a smoothed harmonic magnitude, ĉk(t), for each of the N harmonics;
smooth the instant phase, pk(t), to generate a smoothed harmonic phase, {circumflex over (p)}k(t), for each of the N harmonics;
generate a magnitude perturbation signal, dck(t), for each of the N harmonics;
generate a phase perturbation signal, dpk(t), for each of the N harmonics; and
generate an individual harmonic signal, hk(t), for each of the N harmonics according to: hk(t)=(ĉk(t)+dck(t))·sin(2πflkt+{circumflex over (p)}k(t)+dpk(t));
sum the individual harmonic signals, hk(t), and thereby generating an engine harmonic enhancement signal, he(t) to be added on top of a baseline engine sound; and
transduce the engine harmonic enhancement signal, he(t), to acoustic energy.

5. The engine harmonic enhancement system of claim 4, wherein, when executed, the instructions cause the one or more processors to generate the magnitude perturbation signal, dck(t), for each of the N harmonics according to:




  
    
      d
      ⁢
      
        
      
      ⁢
      
        
          c
          k
        
        ⁡
        
          (
          t
          )
        
      
    
    =
    
      
        
          ncf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              c
              k
            
          
        
        
          σ
          
            ncf
            k
          
        
      
    
  
  ,



where σΔck, is the standard deviation of fluctuations in the harmonic magnitude of each of the N harmonics of the target engine sound; ncfk(t) is filtered random noise used for generating the magnitude perturbation signal; and σncfk is the standard deviation for the filtered random noise used for generating the magnitude perturbation signal.

6. The engine harmonic enhancement system of claim 4, wherein, when executed, the instructions cause the one or more processors to generate the phase perturbation signal, dpk(t), for each of the N harmonics according to:




  
    
      
        dp
        k
      
      ⁡
      
        (
        t
        )
      
    
    =
    
      
        
          npf
          k
        
        ⁡
        
          (
          t
          )
        
      
      ·
      
        
          σ
          
            Δ
            ⁢
            
              
            
            ⁢
            
              p
              k
            
          
        
        
          σ
          
            npf
            k
          
        
      
    
  
  ,



where σΔpk, is the standard deviation of fluctuations in the harmonic phase of each of the N harmonics of the target engine sound; npfk(t) is filtered random noise used for generating the phase perturbation signal; and σnpfk is the standard deviation for the filtered random noise used for generating the phase perturbation signal.