Abstract:
A method for adaptively learning a sparse impulse response ( 100 ) of a continuous channel to which an input signal (x (t)) is applied and which delivers an output signal (y(t)), comprising the following steps:
       low-pass filtering the input signal and the output signal and obtain a filtered input signal (x F (t)) and a filtered output signal (y F (t))   sampling the filtered input signal and the filtered output signal with a sampling rate below the Nyquist rate and obtaining a sampled input signal (x S (t)) and a sampled output signal (y S (t))   retrieving from the sampled input signal (x S (t)) and the sampled output signal (y S (t)) an estimate ( 400 ) of the sparse impulse response ( 100 ) of the continuous channel.       
 
     This method can be applied in CDMA channels, in acoustic room context, in ultra-wideband ranging and line echo cancellation problems, in transmission systems for optical fibres, in body scan devices, to name a few.

Description:
RELATED APPLICATIONS 
       [0001]    This application claims the benefit of U.S. Provisional Application No. 61/255,765, filed on Oct. 28, 2009, which is hereby incorporated by reference in its entirety. 
     
    
     TECHNICAL FIELD 
       [0002]    This disclosure relates to a method and an apparatus for adaptively learning a sparse impulse response of a continuous channel. 
         [0003]    This method and apparatus can be applied to various situations where an adaptive learning of a sparse impulse response of a continuous channel is required, such as without restrictions in some wireless communication channels, such as e.g. CDMA channels, in acoustic room context, in ultra-wideband ranging and line echo cancellation problems, in transmission systems for optical fibers, in body scan devices, to name a few. 
       BACKGROUND 
       [0004]    In this context the adjective “continuous” when applied to a channel and in particular to its impulse response, means continuous in the time-domain, or continuous-time, i.e. a varying quantity whose domain, in this case the time, is a continuum, e.g. a connected interval of real numbers. In other words, the function&#39;s domain is an uncountable set. The impulse response function itself need not be continuous. To contrast, the impulse response of a discrete-time channel has a countable domain, like natural numbers. The impulse response of a discrete-time channel can be obtained by sampling the impulse response of a continuous-time channel. 
         [0005]    In each of the above cited applications, the continuous channel is unknown although it is known that it is sparse in time, i.e. it consists of few well localised signals, for example due to echoes in an acoustic room. An example of a continuous sparse channel is a channel whose impulse response h(t) can be modeled as a linear combination of several Diracs, i.e. 
         [0000]    
       
         
           
             
               
                 
                   
                     h 
                      
                     
                       ( 
                       t 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       K 
                     
                      
                     
                         
                     
                      
                     
                       
                         c 
                         k 
                       
                        
                       
                         δ 
                          
                         
                           ( 
                           
                             t 
                             - 
                             
                               t 
                               k 
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where {c k } k=1   K  and {t k } k=1   K  are some unknown parameters. K is the sparsity of the channel. 
         [0006]    In the context of the one embodiment of the invention the function h(t) is the impulse response of the continuous-time channel. Since a channel can be modeled as a filter, in this context we also refer h(t) as the filter that models the continuous-time channel. Similar considerations are valid for the other cited channel, e.g. the discrete-time channel. An example of an impulse response of the continuous-time channel h(t) is illustrated in  FIG. 3 , which will be discussed later. 
         [0007]    It is often wanted to learn this unknown continuous sparse channel by building an estimate of this filter h(t). Since the channel is changing in time, it is often necessary to modify this channel iteratively in order to adapt this estimate to the variations of the continuous sparse channel in time. 
         [0008]    A possible system comprising an automatic and adaptive channel estimate is illustrated in  FIG. 1 , which represents a scenario involving hands-free telephony in an acoustic and close room  1 . In this figure, s 1  represents a known voice signal coming from a party at a remote side of the conversation and played back with a loudspeaker  2 . A speaker  4  is speaking by using the microphone  6 ; s 2  represents the acoustic signal which is picked up by the microphone  6 . However the microphone  6  picks up also a distorted copy of s 1 : the distortion is due to the room reverberations and can be modeled as the convolution between s 1  and an unknown filter H, indicated with the reference  100  and corresponding to the acoustic channel between the source s 1  and the microphone  4 . In other words the microphone  6  measures both the speaker voice s 2  and the result of above convolution: 
         [0000]      s 2 +s 1 *H  (2)
 
         [0009]    The goal of the system is to transmit to the party at the remote side only the signal s 2 , so that this party will not hear the echoes s 1 *H of her/his own voice. This task is also known as “acoustic echo cancellation” and it is achieved by using an adaptive filter Ĥ, indicated with the reference  400 , which represents an estimate of the unknown filter H  100 . 
         [0010]    The signal s 1  is also the input of the filter Ĥ  400  and its output is indicated by ŝ 1 . The closer the filter Ĥ to the unknown filter H, the more ŝ 1  is equal to s 1 . For this reason ŝ 1  is subtracted (reference  50 ) from the output of the microphone given by the above expression (2): in the ideal case, i.e. in the case wherein the filter Ĥ exactly reproduced the unknown filter H 
         [0000]      ( s   2   +s   1   *H )− ŝ   1   =s   2   +s   1   *H−s   1   *Ĥ   H=Ĥ   =s   2   (3)
 
         [0011]    The arrow of the filter Ĥ  400  in  FIG. 1  indicates that it can be iteratively modified in order to adapt Ĥ to temporal variations of H. 
         [0012]    The example of  FIG. 1  has been introduced for focusing the problem. However the scope of one embodiment of the invention is not limited to the described example: as already mentioned, the man skilled in the art can imagine other scenarios in which it is necessary to adaptively learn a sparse impulse response of a continuous, almost shift-invariant channel. 
         [0013]    Generally speaking, a linear channel that is “almost” shift-invariant, i.e. that is slowly varying in time, is an example of a continuous—and unknown—channel for which it is useful to learn its sparse impulse response. How “slow” a channel varies is relative to the tracking speed of the channel identification methods. In general, the channel should remain roughly constant within N taps, where N is the number of steps taken by an adaptive method for converging to a solution. 
         [0014]    Computing a filter estimate Ĥ is based on linear system identification, which is a classical and well-studied subject. A known method for learning a general sparse impulse response h(t) of a continuous channel can be summarized as follows. x(t) and y(t) are the input and output signal of the unknown continuous channel, respectively. The input x(t) is supposed to be band-limited, so that x(t) and y(t) can be converted into the discrete domain by sampling them at or above the Nyquist rate. 
         [0015]    Since the sampling process is alias-free, the resulting discrete signals—denoted by x[n] and y[n] respectively—are also related by a linear shift-invariant system, i.e. 
         [0000]        y[n]= ( x*h   d ) [n]   (4)
 
         [0016]    The discrete filter h d  [n] above either completely determines the original continuous-time impulse response h(t) if the latter is band-limited and has a bandwidth below that of x(t), or is a low-pass approximation of h(t). In practice, h d [n] is usually approximated by a finite impulse response (FIR) filter with N taps, where N is chosen to be large enough to retain most of the energy of h d [n]. In other words after sampling, the discrete filter h d  [n] is approximated by an FIR model, whose length, N, must be large enough to capture the longest possible delay in the system. More specifically: 
         [0000]    
       
         
           
             
               
                 
                   N 
                   ≈ 
                   
                     
                       
                         max 
                         k 
                       
                        
                       
                         t 
                         k 
                       
                     
                     T 
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where T is the sampling interval. N is generally much larger than the number of delays K that make up the system itself. 
         [0017]    This finite-order system can then be learned from x[n] and y[n] by using standard adaptive filtering techniques, such as the least mean squares (LMS) algorithm. 
         [0018]    Despite its widespread and successful use in practice for learning general linear time-invariant systems, the traditional approach described above is computationally expensive, since all the N filter taps have to be determined. Moreover the bandwidth of the input signal x(t) needs to be limited, so that x(t) and y(t) can be converted into the discrete domain by sampling them at or above the Nyquist rate. However there are wide classes of very common signals, including stream of Dirac pulses, bi-level signals, piecewise polynomial signals, etc. which are not band-limited and which cannot be used for estimating the channel with this known method. An example of a system using such no band-limited signals is the CDMA. 
         [0019]    Even if the signals are band-limited, their bandwidth is sometime very large, so that the current methods requires fast, expensive and power-consuming A/D converter, fast digital circuits and a waste storage place for storing the digitized signals while fulfilling the Nyquist condition. Examples of signals with a large bandwidth that makes them impracticable for the above mentioned methods include UWB signals. 
         [0020]    Moreover the known standard adaptive filtering techniques are not efficient, since when either the input signal x(t) or the channel noise is not white, there is no possibility to select the frequency band of the input signal x(t) with the highest SNR (Signal to Noise Ratio). 
         [0021]    Known algorithm which adaptively obtains the impulse response of the filter Ĥ  400  include the Frequency Least Mean Square (FLMS), which operates in the frequency-domain, the Recursive Least Square (RLS) and the reweighted zero attracting method (RZA), both operating in the time-domain. An example of RZA method is given by Y. Chen, Y. Gu, and A. O. Hero, “ Sparse LMS for system identification,” in Proc. IEEE Int. Conf Acoust., Speech, and Signal Proc., Taipei,  2009. Another LMS sparsity-exploiting method operating in the time domain is given by R. K. Martin, W. A. Sethares, R. C. Williamson, and C. R. Johnson, “ Exploiting sparsity in adaptive filters,” IEEE Trans. Signal Process ., vol. 50, no. 8, pp. 1883-1894, August 2002. 
         [0022]    If for comparison purposes the parameters of each algorithm are made equal when possible, the known FLMS, compared to the RLS and the RZA, has the highest Mean-Squared output Error (MSE). The RLS is not efficient to track a change of the unknown filter taps. The RZA is efficient in tracking these changes, it presents a MSE better than a classical FLMS but it is computationally complex and requires a sampling rate above the Nyquist system, i.e. the convergence rate is quite low. Moreover RZA does not handle well the modelling of the sparsity off of the sampling grid since it assume that the discrete channel is exactly K-sparse. However, if the continuous dirac delta (t−t_k) falls on a fractional number (i.e. off the discrete grid), one can see a sinc function in the discrete domain, which is not exactly sparse. Finally RLS and RZA methods are computationally less efficient than frequency domain methods, if block-based processing is used. 
         [0023]    WO10083321 describes a method for reconstructing a pair of signals, wherein the first signal is linked to the second signal by an unknown but sparse filter. In this case the sparseness information is used for reconstructing the two signals, but not the filter itself Moreover the method assumes that the unknown channels are constant over time, i.e. it cannot track slowly-varying channels. 
         [0024]    A method and an apparatus for adaptively learning a sparse impulse response of a continuous channel which does not necessarily require a band-limited input signal x(t) and with a lower computational complexity than known methods are needed. 
         [0025]    A method and an apparatus for adaptively learning a sparse impulse response of a continuous channel which are more efficient in the case when either the input signal x(t) or the channel noise is not white are also needed. 
         [0026]    A new FLMS method which is efficient to track a change of the unknown filter taps, which presents a MSE better than a classical FLMS, which is more versatile and less computationally complex than other known method is also needed. 
       SUMMARY 
       [0027]    In general, this disclosure describes techniques for adaptively learning a sparse impulse response of a continuous channel. 
         [0028]    The approach described above does not exploit the special sparse structure of the unknown continuous channel. Intuitively, since the impulse response in (1) can be specified by only a small number of parameters, i.e. 2K, one should expect a much more efficient scheme in learning h(t). 
         [0029]    By exploiting the sparse structure of the continuous channel, according to one embodiment of the invention it is possible to reduce the sampling rate below the Nyquist rate, and still recover the entire impulse response of h(t). Reducing the sampling rate below the Nyquist rate allows to use less expensive and power-consuming digital circuits and to reduce the storage place for storing the digitized signals. Moreover it is not necessary that the input signal is band-limited, so that the new method can also be used to non band-limited signals, such as CDMA signals. 
         [0030]    The sparse system according to one embodiment of the invention can be learned by an estimate, for example an FIR filter, which has roughly the same order as the sparsity level K, where K&lt;&lt;N, e.g. K=12 and N=1024. 
         [0031]    The new adaptive FLMS method according to the one embodiment of invention uses a sub-Nyquist sampling step, i.e., a step of sampling a signal at a rate which is below the sampling rate given by Nyquist. The new method is based on the observation that the Fourier transform of the sparse impulse response of a continuous-time channel is a sum of complex exponentials. By borrowing tools from spectral analysis and array signal processing, it is possible to use only a small band of the Fourier transform to recover its frequencies, phases, and amplitude, the collection of which parameterizes the entire channel. 
         [0032]    In other words the proposed method allows to select a frequency band of the input signal x(t) with the highest signal to noise ratio (SNR), which is useful in the case when either the input signal x(t) or the channel noise is not white. 
         [0033]    According to one example of the invention, the original continuous-time channel can be determined by considering a sampled discrete-time system, whose sampling rate can vary freely below the Nyquist rate. Moreover, with sufficient sub-sampling, this equivalent discrete system can be learned by running LMS for an FIR filter with only 2K taps. Let T be the sampling interval, i.e. 1/T the sampling rate, L the maximum span, i.e. the delay spread in the time domain, of the original continuous filter and N the number of taps of the discrete filter, then N=L/T. By increasing T, i.e. by reducing the sampling rate, it is possible to reduce N until N is roughly equals 2K. 
         [0034]    By exploiting the prior knowledge of the complex exponential structure of the continuous channel in the frequency domain, a new frequency-domain LMS (FLMS) method with significantly improved convergence performance can be built. The parameters of the original continuous-time channel can then be recovered by a new method based on different tools known in the field of harmonic retrieval. 
         [0035]    In one embodiment, the continuous-time channel is a wireless RF channel and its input output signals are electromagnetic signals. 
         [0036]    In one example, a method for adaptively learning a sparse impulse response of a continuous channel to which an input signal is applied and which delivers an output signal includes:
       low-pass filtering said input signal and said output signal and obtaining a filtered input signal and a filtered output signal   sampling said filtered input signal and said filtered output signal with a sampling rate below the Nyquist rate and obtaining a sampled input signal and a sampled output signal   retrieving from said sampled input signal and said sampled output signal an estimate of said sparse impulse response of said continuous channel.       
 
         [0040]    In another example a computer-readable medium, such as a computer-readable storage medium, for adaptively learning a sparse impulse response of a continuous channel to which an input signal is applied and which delivers an output signal, is encoded with instructions that cause a programmable processor to
       low-pass filter said input signal and said output signal and obtain a filtered input signal and a filtered output signal   sample said filtered input signal and said filtered output signal with a sampling rate below the Nyquist rate and obtain a sampled input signal and a sampled output signal   retrieve from said sampled input signal and said sampled output signal an estimate of said sparse impulse response of said continuous channel.       
 
         [0044]    In another example, an apparatus for adaptive learning a sparse impulse response of a continuous channel having an input signal and an output signal includes
       means for low-pass filtering said input signal and said output signal and obtaining a filtered input signal and a filtered output signal   means for sampling said filtered input signal and said filtered output signal with a sampling rate below the Nyquist rate and obtaining a sampled input signal and a sampled output signal   means for retrieving from said sampled input signal and said sampled output signal an estimate of said sparse impulse response of said continuous channel.       
 
         [0048]    In another example, an apparatus for adaptive learning a sparse impulse response of a continuous channel having an input signal and an output signal includes
       a low-pass filter arranged to low-pass filter said input signal and said output signal and obtain a filtered input signal and a filtered output signal   a sampler arranged to sample said filtered input signal and said filtered output signal with a sampling rate below the Nyquist rate and obtain a sampled input signal and a sampled output signal   a calculator configured to retrieve from said sampled input signal and said sampled output signal an estimate of said sparse impulse response of said continuous channel.       
 
         [0052]    In one embodiment the apparatus can be a radio-transmitter comprising also circuits configured to transmit said input signal and receive said output signal. 
         [0053]    In one embodiment, the radio-transmitter is a mobile phone. 
         [0054]    In another embodiment the apparatus can be an acoustic echo canceller, comprising also circuits configured to cancel the acoustic echo from the continuous-time channel. 
         [0055]    In another embodiment the apparatus can be a line echo canceller, comprising also circuits configured to cancel the line echo from the continuous-time channel, e.g. an optical fiber. 
         [0056]    The details of one or more examples are set forth in the accompanying drawings and the description below. Other features, objects, and advantages will be apparent from the description and drawings, and from the claims. 
     
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         [0057]      FIG. 1  is a block diagram illustrating an acoustic system. 
           [0058]      FIG. 2  is a block diagram illustrating a system for adaptive learning a sparse impulse response of a continuous channel. 
           [0059]      FIG. 3  is a chart illustrating a continuous-time sparse channel. 
       
    
    
     DETAILED DESCRIPTION 
       [0060]      FIG. 1  is a block diagram illustrating an acoustic system. 
         [0061]    An example of such impulse response is illustrated in  FIG. 3 : the first Dirac  102  could correspond for example to the first reflection of an acoustic signal against a first wall of a close acoustic room, the second as the second reflection against a second wall, etc. The last part of the  FIG. 3 , denoted by the reference  104  and comprising signals having a small amplitude with regard to the Diracs  102 , changes fast, while the first part constituted by Diracs  102  slowly changes. In the context of the one embodiment of the invention, only the first part of this impulse response, constituted by a number of Diracs equal to the sparsity K of the channel, is taken into account in the proposed method, and the last part  104  of  FIG. 3  is disregarded. 
         [0062]    In  FIG. 1  a frequency modulation  10  is applied with a frequency modulator to the input signal x(t) and to the output signal y(t): in the illustrated example these input and output signals are multiplied by e −jΩ     0     t . A frequency modulated input signal x M (t) and a frequency modulated output signal y M (t) are then respectively obtained. 
         [0063]    In one embodiment these frequency modulated signals are then passed through a low-pass filter  20 , defined as 
         [0000]        L (Ω)=1 [−π/T,π/T] (Ω)  (6)
 
         [0000]    where 1         (•) denotes the indicator function defined on an interval F. The specification of this filter guarantees that there is not aliasing in the next sampling process. The filter  20  can be an electrical or software filter. 
         [0064]    Subsequently, the two filtered signals are sampled at a sampling rate 1/T below the Nyquist rate (reference  30  in  FIG. 1 ) in order to obtain a sampled input signal x[n] and a sampled output signal y[n] according to the following definitions 
         [0000]        x   T   [n]   def   =x   F ( nT ) and  y   T   [n]   def   =y   F ( nT )  (7)
 
         [0000]    where x F  and y F  are respectively the filtered input and output signal of in  FIG. 2 . Those methods can be regarded as a special case of the scheme of  FIG. 1 , when the modulation frequency Ω 0  is set to 0, the low-pass filter is always present and the sampling rate 1/T is equal to the Nyquist rate of the filtered input signal x F (t). 
         [0065]    According to one aspect of the invention, the discrete signals x T [n] (sampled input signal) and y T [n] (sampled output signal) are related by a discrete-time filtering operation: 
         [0000]        y   T   [n]= ( x   T   *h   T ) [n]   (8)
 
         [0000]    where the filter h T [n], illustrated by a dashed line in  FIG. 2  (reference  200 ), has a discrete-time Fourier transform H T (e jω ) given by 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       H 
                       T 
                     
                      
                     
                       ( 
                       
                          
                         jω 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       K 
                     
                      
                     
                       
                         c 
                         k 
                       
                        
                       
                          
                         
                           
                             - 
                             
                               jΩ 
                               0 
                             
                           
                            
                           
                             t 
                             k 
                           
                         
                       
                        
                       
                          
                         
                           
                             - 
                             jω 
                           
                            
                           
                               
                           
                            
                           
                             
                               t 
                               k 
                             
                             / 
                             T 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
         [0066]    In other words, the estimate in the frequency domain of the unknown continuous channel can be modeled as a filter comprising a sum of complex exponentials. See the appendix A for further details. 
         [0067]    The sampled discrete-time channel h T [n] specified in (8) and (9) can completely determine the original continuous-time sparse channel h(t). See the appendix B for a demonstration. 
         [0068]    The proposition of the previous paragraph is valid for an arbitrary choice of T and Ω. Thus, by increasing T, the sampling rate for the equivalent discrete-time channel h T [n] can be set to well below the Nyquist rate. As outlined, reducing the sampling rate below the Nyquist rate allows to use less expensive and power-consuming digital circuits and to reduce the storage place for storing the digitized signals. Moreover it is not necessary that the input signal is band-limited, so that the method can be applied no non band-limited signals, such as CDMA for example. 
         [0069]    Additionally, the modulation frequency Ω 0  can be used to select a frequency band of the input signal x(t) with the highest possible signal to noise ratio (SNR). As discussed, this is useful when either the input signal x(t) or the noise is not white. 
         [0070]    In  FIG. 1  the noise v[n] has been introduced for describing a realistic and no-ideal situation: it is added at element  40  to the output sampled signal y[n]. 
         [0000]    Learning the Discrete-Time Channel h T [n] 
         [0071]    In order to learn the discrete-time channel h T [n] that, as demonstrated in Appendix B, is equivalent to the continuous-time channel h(t), a new frequency-domain LMS (FLMS) method is used. The idea is to approximate h T [n], which is of infinite length, by a truncated FIR filter {tilde over (h)} T [n]. It is possible to show that the discrete Fourier transform (DFT) of this truncated filter {tilde over (h)} T [n] an be approximated by a discretized version of the discrete-time Fourier transform (9) of the discrete-time channel h T [n]. 
         [0072]    From (9) it follows that the discrete-time channel h T [n] has the following form: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       h 
                       T 
                     
                      
                     
                       [ 
                       n 
                       ] 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       K 
                     
                      
                     
                       
                         c 
                         k 
                       
                        
                       
                          
                         
                           
                             - 
                             
                               jΩ 
                               0 
                             
                           
                            
                           
                             t 
                             k 
                           
                         
                       
                        
                       
                         
                           sin 
                            
                           
                             ( 
                             
                               π 
                                
                               
                                 ( 
                                 
                                   n 
                                   - 
                                   
                                     
                                       t 
                                       k 
                                     
                                     / 
                                     T 
                                   
                                 
                                 ) 
                               
                             
                             ) 
                           
                         
                         
                           π 
                            
                           
                             ( 
                             
                               n 
                               - 
                               
                                 
                                   t 
                                   k 
                                 
                                 / 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
         [0073]    Unless the t k  parameters are all integer multiples of T, i.e. the sampling interval, the discrete-time channel h T [n] of (10) can be modeled as a filter of infinite length. It is then approximated by an FIR filter of length N, where 
         [0000]    
       
         
           
             
               
                 
                   N 
                   ≥ 
                   
                     
                       ⌈ 
                       
                         
                           
                             max 
                             k 
                           
                            
                           
                             ( 
                             
                               t 
                               k 
                             
                             ) 
                           
                         
                         T 
                       
                       ⌉ 
                     
                     + 
                     c 
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
         [0000]    wherein ┌ ┐ indicates the ceiling function and c some constant. Even for relatively small value of c, the above truncated FIR filter, denoted by {tilde over (h)} T [n], can capture most of energy of h T [n]. For example c=5 include at least 99% of the energy for the largest t k  in the worst case. 
         [0074]    As mentioned the DFT of the truncated filter {tilde over (h)} T [n], denoted by H[m], can be approximated by a discretized version of the discrete-time Fourier transform H T (e jω ) given in (9). An example of this discretized version is 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       H 
                       m 
                     
                     = 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                        
                       
                         
                           c 
                           k 
                         
                          
                         
                            
                           
                             
                               - 
                               
                                 jΩ 
                                 0 
                               
                             
                              
                             
                               t 
                               k 
                             
                           
                         
                          
                         
                            
                           
                             
                               - 
                               j2π 
                             
                              
                             
                                 
                             
                              
                             
                               
                                 
                                   t 
                                   k 
                                 
                                  
                                 m 
                               
                               TN 
                             
                           
                         
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       m 
                       = 
                       0 
                     
                     , 
                     1 
                     , 
                     … 
                      
                     
                         
                     
                     , 
                     
                       N 
                       - 
                       1 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
         [0075]    The filter of the formula (12) has a parametric structure, since it depends from the parameters c k  and t k . By increasing the sampling interval T, not only the sampling rate 1/T becomes lower and below the Nyquist rate, but also the model order N can be decreased, according to the constraint in (11). As will be discussed, a possible lower bound of N can be 2K+2. 
         [0076]    For learning the truncated FIR filter {tilde over (h)} T [n], or equivalently its approximated DFT {H m } m=0   N−1  given in (11) many type of adaptive filters can be used. The frequency-domain LMS (FLMS) in particular can learn the approximated DFT {H m } m=0   N−1  directly: at each iteration it finds the filter coefficients that relate to producing the least mean squares of the error signal, e.g. the difference performed by the element  50  in  FIG. 1 . It is a stochastic gradient descent method in that the filter is only adapted based on the error at the current time. In other words at each iteration, an estimate of {H m } m=0   N−1  is made by following a gradient descent step. This estimate is indicated by 
         [0000]      Ĥ=[Ĥ 0 , . . . , Ĥ N−1 ] T   (13)
 
         [0077]    The initial value of this estimate can be a null-vector or any other vectors. The new FLMS according to the one embodiment of the invention exploits the special parametric structure of the filter given by (12). In a first step a Toeplitz matrix A, i.e. a matrix in which each descending diagonal from left to right is constant, is formed according to 
         [0000]    
       
         
           
             
               
                 
                   A 
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               H 
                               L 
                             
                           
                           
                             
                               H 
                               
                                 L 
                                 - 
                                 1 
                               
                             
                           
                           
                             … 
                           
                           
                             
                               H 
                               0 
                             
                           
                         
                         
                           
                             
                               H 
                               
                                 L 
                                 + 
                                 1 
                               
                             
                           
                           
                             
                               H 
                               L 
                             
                           
                           
                             … 
                           
                           
                             
                               H 
                               1 
                             
                           
                         
                         
                           
                             … 
                           
                           
                             … 
                           
                           
                             … 
                           
                           
                             … 
                           
                         
                         
                           
                             
                               H 
                               N 
                             
                           
                           
                             
                               H 
                               
                                 N 
                                 - 
                                 1 
                               
                             
                           
                           
                             … 
                           
                           
                             
                               H 
                               
                                 N 
                                 - 
                                 L 
                               
                             
                           
                         
                       
                       ] 
                     
                     
                       
                         ( 
                         
                           N 
                           - 
                           L 
                           + 
                           1 
                         
                         ) 
                       
                       × 
                       
                         ( 
                         
                           L 
                           + 
                           1 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where L is the floor function of N/2, i.e. 
         [0000]    
       
         
           
             L 
             = 
             
               
                 ⌊ 
                 
                   N 
                   2 
                 
                 ⌋ 
               
               . 
             
           
         
       
     
         [0078]    From J. A. Cadzow, “ Signal enhancement—A composite property mapping algorithm,” IEEE Trans. Acoust., Speech, and Signal Process ., vol. 36, no. 1, pp. 49-62, January 1988 it is known that in the noiseless case the rank of A must be equal to K. This known property can be used for accelerating the convergence of the new FLMS method, which can be called Cadzow-FLMS or C-FLMS, as summarized in the following steps. 
         [0079]    For each iteration of the new FLMS method or C-FLMS, the current estimate of the approximated DFT of the truncated FIR filter {tilde over (h)} T [n] which approximates the discrete-time filter, given in (12), is denoted by Ĥ n . The new FLMS method or C-FLMS for each iteration is composed by the following steps:
       Finding a standard update vector of the C-FLMS, indicated by S, by using for example the method described in Ali H. Sayed,  Adaptive Filters , John Wiley &amp; Sons, 2008, p. 445 and adding it to a current estimate of the Fourier coefficient Ĥ n , in order to find the next estimate of the Fourier coefficient Ĥ n+1  according to:       
 
         [0000]        Ĥ   n+1   =Ĥ   n   +S   (15)
       Building a Toeplitz matrix A from the next estimate Ĥ n+1      In order to enforce the rank condition of A (rank A =K) computing the Singular Value Decomposition (SVD) of the Toeplitz matrix A, which is a particular kind of factorization based on the use of eigen-values and eigen-vectors, used for generating an approximation of the matrix A having a minor rank according to       
 
         [0000]      A=UΛV T   (16)
           where
               U is an unitary matrix (N−L+1)×(N−L+1)   Λ is a diagonal matrix (N−L+1)×(L+1)   V T  is the conjugate transpose matrix of a unitary matrix (L+1)×(L+1)   
                   Building a diagonal matrix Λ′ from the diagonal matrix Λ the SVD by zerogoing all term less than the K most significative, where K is the sparsity of the continuous channel.   Producing a new matrix A′ having rank K by using this built diagonal matrix Λ′ according to       
 
         [0000]      A′=UΛ′V T   (17)
           This step allows to build a matrix with a rank equal to the sparsity of the continuous channel and it is important in the context of the one embodiment of the invention, since the method at this step promotes and exploits the sparsity of the continuous channel.       Computing a denoised version C, C being a vector, of the next estimate Ĥ n+1  by averaging the diagonals of the new matrix A′.
           Intuitively, this “SVD denoising step” brings the next estimate closer to the model assumption made in the formula (12), i.e. the approximation of the DFT of the truncated filter {tilde over (h)} T [n], denoted by H[m], by a discretized version of the discrete-time Fourier transform H T [e jω ] given in (9).   
           Updating the next estimate Ĥ n+1  taking into account the update vector S and the denoised version C according to       
 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       H 
                       ^ 
                     
                     
                       n 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       
                         H 
                         ^ 
                       
                       n 
                     
                     + 
                     S 
                     + 
                     
                       
                         ( 
                         
                           I 
                           - 
                           
                             
                               SS 
                               H 
                             
                             
                               
                                 S 
                                 H 
                               
                                
                               S 
                             
                           
                         
                         ) 
                       
                        
                       
                         ( 
                         
                           C 
                           - 
                           
                             
                               H 
                               ^ 
                             
                             n 
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 where S H  represents the Hermitian conjugate of the update vector S. 
                 The matrix 
               
             
           
         
       
     
         [0000]    
       
         
           
             ( 
             
               I 
               - 
               
                 
                   SS 
                   H 
                 
                 
                   
                     S 
                     H 
                   
                    
                   S 
                 
               
             
             ) 
           
         
       
       
         
           
             
               
                 ensures that the inner product between S and the denoised version C is always positive, which is important in guaranteeing the convergence of the method. 
               
             
           
         
       
     
         [0096]    Numerical experiments have shown that the proposed method described above can dramatically improve the convergence rate and the steady state error of the FLMS. 
       Recovering the Original Continuous-Time Sparse Channel 
       [0097]    After finding the approximation H[m] of the DFT of the truncated filter {tilde over (h)} T [n], given in (12), according to the previous steps, the final step is to estimate the parameters c k  and t k , which can be used for reconstructing the original continuous-time sparse impulse response h(t). 
         [0098]    Different methods can be used for this parameter estimation task. A possible choice can be the following, which is called the Prony&#39;s method:
       From the Fourier coefficient {Ĥ m } m=0   N−1  learned by the previous steps forming the Toeplitz matrix A as in (14). In this case the parameter L has to be set equal to K.   Performing the Singular Value Decomposition SVD of the built Toeplitz matrix A according to the previous formula (17).   Taking the smallest possible right singular vector a=[a 0 , a 1 , . . . , a k ] T , i.e. the column of V with the smallest possible singular value   Forming a K-th-order polynomial whose coefficients are elements of a. The K zeros or roots of this polynomial lie directly on the complex exponential frequency of (12).   Forming a list of 2π-normalized frequencies by taking the complex angles of these roots.       
 
         [0000]    
       
         
           
             NT 
             
               2 
                
               π 
             
           
         
       
       
         
           
             Obtaining the values {t k } k=1   K  by multiplying this list by 
             Obtaining the values {c k } k=1   K  by the least-square solution of (12). 
           
         
       
     
         [0106]    In one embodiment the means for low-pass filtering comprise a hardware-implemented low-pass filter or a software-implemented low-pass filter. 
         [0107]    In one embodiment the means for sampling comprise a hardware-implemented sampler or a software-implemented sampler. 
         [0108]    In one embodiment the means for retrieving comprise at least one processor, such as one or more digital signal processors (DSPs), general purpose microprocessors, application specific integrated circuits (ASICs), field programmable logic arrays (FPGAs), or other equivalent integrated or discrete logic circuitry. 
         [0109]    In one or more examples, the functions described may be implemented in hardware, software, firmware, or any combination thereof. If implemented in software, the functions may be stored on or transmitted over as one or more instructions or code on a computer-readable medium. Computer-readable media may includes computer data storage media or communication media including any medium that facilitates transfer of a computer program from one place to another. Data storage media may be any available media that can be accessed by one or more computers or one or more processors to retrieve instructions, code and/or data structures for implementation of the techniques described in this disclosure. By way of example, and not limitation, such computer-readable media can comprise RAM, ROM, EEPROM, CD-ROM or other optical disk storage, magnetic disk storage or other magnetic storage devices, or any other medium that can be used to carry or store desired program code in the form of instructions or data structures and that can be accessed by a computer. Also, any connection is properly termed a computer-readable medium. For example, if the software is transmitted from a website, server, or other remote source using a coaxial cable, fiber optic cable, twisted pair, digital subscriber line (DSL), or wireless technologies such as infrared, radio, and microwave, then the coaxial cable, fiber optic cable, twisted pair, DSL, or wireless technologies such as infrared, radio, and microwave are included in the definition of medium. Disk and disc, as used herein, includes compact disc (CD), laser disc, optical disc, digital versatile disc (DVD), floppy disk and Blu-ray disc where disks usually reproduce data magnetically, while discs reproduce data optically with lasers. Combinations of the above should also be included within the scope of computer-readable media. 
         [0110]    The code may be executed by one or more processors, such as one or more digital signal processors (DSPs), general purpose microprocessors, application specific integrated circuits (ASICs), field programmable logic arrays (FPGAs), or other equivalent integrated or discrete logic circuitry. Accordingly, the term “processor,” as used herein may refer to any of the foregoing structure or any other structure suitable for implementation of the techniques described herein. In addition, in some aspects, the functionality described herein may be provided within dedicated hardware and/or software modules configured for encoding and decoding, or incorporated in a combined codec. Also, the techniques could be fully implemented in one or more circuits or logic elements. 
         [0111]    The techniques of this disclosure may be implemented in a wide variety of devices or apparatuses, including a wireless handset, an integrated circuit (IC) or a set of ICs (i.e., a chip set). Various components, modules or units are described in this disclosure to emphasize functional aspects of devices configured to perform the disclosed techniques, but do not necessarily require realization by different hardware units. Rather, as described above, various units may be combined in a codec hardware unit or provided by a collection of interoperative hardware units, including one or more processors as described above, in conjunction with suitable software and/or firmware. 
         [0112]    Various examples have been described. These and other examples are within the scope of the following claims. 
         [0113]    It is to be understood that the claims are not limited to the precise configuration and components illustrated above. Various modifications, changes and variations may be made in the arrangement, operation and details of the methods and apparatus described above without departing from the scope of the claims. 
       APPENDIX A 
       [0114]    Proof of the Form of the Discrete-Time Fourier Transform of the Discrete-Time Filtering Operation Between the Input Sampled Signal and the Output Sampled Signal 
         [0115]    The specification of the low-pass filter in (6) guarantees that there is no aliasing in the sampling process. Consequently, the discrete-time Fourier transform of x T [n] can be obtained as 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       X 
                       T 
                     
                      
                     
                       ( 
                       
                          
                         jω 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         1 
                         T 
                       
                        
                       
                         X 
                          
                         
                           ( 
                           
                             
                               ω 
                               T 
                             
                             + 
                             
                               Ω 
                               0 
                             
                           
                           ) 
                         
                       
                        
                       
                          
                         ω 
                          
                       
                     
                     ≤ 
                     π 
                   
                 
               
               
                 
                   ( 
                   A1 
                   ) 
                 
               
             
           
         
       
     
         [0116]    where X(Ω) denotes the continuous-time Fourier transform of x(t). Similarly for y T [n] 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         Y 
                         T 
                       
                        
                       
                         ( 
                         
                            
                           jω 
                         
                         ) 
                       
                     
                     = 
                     
                       
                         
                           1 
                           T 
                         
                          
                         
                           Y 
                            
                           
                             ( 
                             
                               
                                 ω 
                                 T 
                               
                               + 
                               
                                 Ω 
                                 0 
                               
                             
                             ) 
                           
                         
                       
                       = 
                       
                         
                           1 
                           T 
                         
                          
                         
                           X 
                            
                           
                             ( 
                             
                               
                                 ω 
                                 T 
                               
                               + 
                               
                                 Ω 
                                 0 
                               
                             
                             ) 
                           
                         
                          
                         
                           H 
                            
                           
                             ( 
                             
                               
                                 ω 
                                 T 
                               
                               + 
                               
                                 Ω 
                                 0 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                        
                       ω 
                        
                     
                     ≤ 
                     π 
                   
                 
               
               
                 
                   ( 
                   A2 
                   ) 
                 
               
             
           
         
       
     
         [0117]    where H(Ω) denotes the Fourier transform of h(t). From the definition of h(t) in (1) 
         [0000]    
       
         
           
             
               
                 
                   
                     H 
                      
                     
                       ( 
                       Ω 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       K 
                     
                      
                     
                       
                         c 
                         k 
                       
                        
                       
                          
                         
                           
                             - 
                             jΩ 
                           
                            
                           
                               
                           
                            
                           t 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   A3 
                   ) 
                 
               
             
           
         
       
     
         [0118]    Substituting the above equality (A 3 ) into (A 2 ) and comparing with (A 1 ) it is possible to verify that 
         [0000]        Y   T ( e   jω )= X   T ( e   jω ) H   T ( e   jω ) 
         [0119]    where H T (e jω ) is as defined in (9). 
       APPENDIX B 
       [0120]    Proof of the Complete Determination of the Original Continuous-Time Sparse Channel h(t) by Using the Sampled Discrete-Time Channel h T [n] as Specified in (9) 
         [0121]    On one hand the parameters {c k } k=1   K  and {t k } k=1   K  from (9) clearly determine the original continuous-time sparse channel h(t). On the other hand, it is known from harmonic retrieval that these parameters can be uniquely determined from at least 2K values of H T (e jω ), with a variety of methods available for that purpose.