Abstract:
A device to perform DFT calculations, for example in a GNSS receiver, including two banks of multipliers by constant integer value, the values representing real and imaginary part of twiddle factors in the DFT. A control unit selectively routes the data through the appropriate multipliers to obtain the desired DFT terms. Unused multipliers are tied to constant input values, in order to minimize dynamic power.

Description:
FIELD OF THE INVENTION 
     The present invention concerns a DFT (Discrete Fourier Transform) processor adapted in particular, but not exclusively, for integration in a receiver for a Global Navigation Satellite Systems (GNSS). Embodiments of the present invention relate to a power and area optimized architecture based on selectively activation of single DFT lines. 
     RELATED ART 
     The Fourier Transform (FT) is a function that converts a signal from the time domain into the frequency domain. In the case of discrete signals of finite duration {x 0 , . . . , x N-1 } the Fourier Transform is often referred to as Discrete Fourier Transform (DFT). 
     Fourier Transforms are applied, among many other applications, in demodulation and processing of GPS, Galileo, GLONASS, and other GNSS signals. In these applications, the Fourier Transform applied to the received data allows to process several carrier frequencies in a parallel fashion, with an important reduction of hardware complexity. Modern GNSS receiver or processors include in most cases a “DFT engine”, that is a section that is especially dedicated to the calculation of DFTs. 
     The theory and details of FT-based GNSS signal processing are known in the art, and will not be discussed at length in the present specification. Exhaustive information can be found in the available literature, for example in the book edited by E. D: Kaplan and C. Hegarty “Understanding GPS and its applications”, 2 nd  edition, Published by Artech House, London (December 2005), which is hereby incorporated by reference. 
     The definition of a DFT point X k  for the mentioned finite sampled signal {x N } is: 
                       X   k     =         ∑     n   =   0       N   -   1       ⁢         x   n     ·     ⅇ       -   j     ⁢       2   ⁢   π     N     ⁢   kn         ⁢           ⁢   k       =   0       ,   1   ,   …   ⁢           ,     N   -   1             (   1   )                 X   k     =       ∑     n   =   0       N   -   1       ⁢       x   n     ·     W   n   kn                 (   2   )               
where W N   kn  are the N th  order complex roots of the unit, also called “twiddle factors” or, writing explicitly the real and imaginary parts of the sum terms in (1):
 
     
       
         
           
             
               
                 
                   
                     
                       x 
                       n 
                     
                     = 
                     
                       I 
                       + 
                       
                         j 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         Q 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       W 
                       N 
                       kn 
                     
                     = 
                     
                       
                         
                           cos 
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   2 
                                   ⁢ 
                                   π 
                                 
                                 N 
                               
                               ⁢ 
                               kn 
                             
                             ) 
                           
                         
                         + 
                         
                           j 
                           · 
                           
                             sin 
                             ⁡ 
                             
                               ( 
                               
                                 
                                   
                                     2 
                                     ⁢ 
                                     π 
                                   
                                   N 
                                 
                                 ⁢ 
                                 kn 
                               
                               ) 
                             
                           
                         
                       
                       := 
                       
                         
                           C 
                           I 
                         
                         + 
                         
                           j 
                           · 
                           
                             C 
                             Q 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       x 
                       n 
                     
                     · 
                     
                       W 
                       N 
                       kn 
                     
                   
                   = 
                   
                     
                       ( 
                       
                         
                           I 
                           · 
                           
                             C 
                             I 
                           
                         
                         - 
                         
                           Q 
                           · 
                           
                             C 
                             Q 
                           
                         
                       
                       ) 
                     
                     + 
                     
                       j 
                       ⁡ 
                       
                         ( 
                         
                           
                             I 
                             · 
                             
                               C 
                               Q 
                             
                           
                           + 
                           
                             Q 
                             · 
                             
                               C 
                               I 
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     Equation (1) can thus be expressed as a sum of DFT terms x n ·W N   kn . The number of DFT terms needed to compute all the X k  points is N 2  but each will require 4 multiplications and 2 additions. 
     An efficient algorithm to calculate the DFT is the well known Fast Fourier Transform (FFT) that is based on a split and conquers approach. If the N (power of 2) samples data stream is halved and processed in parallel, the computation order is reduced to N 2 /2 complex additions and (N 2 /2+N) complex multiplications. As the number of possible splits is equal to log 2 (N) then it follows that the computation order is given by N·log 2 (N) complex additions and N·log 2 (N)/2 complex multiplications. 
     The use of the FFT algorithm is generally considered the most efficient way of calculating N DFT points from N samples. However, there are certain DFT configurations where the FFT algorithm is not optimal. For instance in those applications where only a reduced set of M DFT lines is required (with M≦N), as it is often the case in signal processing and in particular in GNSS processors, the computation order of the FFT architecture is not optimized. 
     Moreover, the FFT algorithm takes its simplest and most efficient form only if N is a power of two. Variant FFT algorithm for an arbitrary N exist, but they are in general less efficient. 
     In a GNSS signal processor, the DFT computation reflects directly on cost, silicon area and power consumption of the receiver. There is therefore a need to provide a DFT algorithm having the lowest possible computation load. 
     The Fourier transform can be regarded as a spectral factorization of a function in the time domain over an orthonormal base of sine and cosine functions. Many other discrete integral transforms are relevant in signal processing techniques, corresponding to different orthonormal bases. These transforms include, for example, the Cosine transform (DCT) and its various modifications (MDCT), the discrete Hartley transform (DHT) and many others. The foregoing specification refers, for simplicity&#39;s sake, to the DFT transform only. It must be understood, however, that the present invention is not limited to this particular case, but include all the discrete integral transforms to which it applies. 
     It is an aim of the present invention to provide a more efficient algorithm to compute a set of DFT lines in a signal processor for processing GNSS signals. The present invention aims moreover to provide a low-power GNSS receiver. 
     BRIEF SUMMARY OF THE INVENTION 
     The goals of the present invention are achieved by the object of the appended independent claims, the variants of the dependent claims incorporating important, but not strictly essential features. The invention will be better understood referring to the detailed description of some embodiments, which is provided by way of example only, and to the drawings that illustrate schematically: 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1 : a known implementation of a DFT computing. 
         FIG. 2 : a pipelined variant of the known device of  FIG. 1 . 
         FIG. 3 : the structure of a DFT engine according to an aspect of the present invention. 
         FIG. 4 : a half butterfly architecture used in the structure of  FIG. 3 . 
         FIG. 5 : a detail of the even coefficient bank, with (K+1)/2 multiplier-by-constant modules. 
         FIG. 6 : a detail of the odd coefficient bank, with (K+1)/2 multiplier-by-constant modules. 
         FIG. 7 : a representation of the search space in GNSS satellite acquisition. 
         FIG. 8 : a GNSS receiver architecture. 
         FIG. 9 : the position, in a particular case, of the twiddle factors of the DFT in the complex plane. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present invention is based on the observation that GNSS processing requires the calculation, in general of a limited number of DFT lines. Twiddle factors used in GNSS DFT engine are a finite set and, moreover, have special symmetry properties in their real and imaginary parts. The invention thus proposes an optimised DFT algorithm that is able to calculate a reduced set M of DFT lines from N samples (with M≦N) being N not necessarily a power of 2, which takes advantage from these symmetries. 
     For these purposes data format considerations, dynamical disabling of non-used DFT lines and routing-effective design partitioning have to be taken into account. 
     A direct DFT implementation, as illustrated on  FIG. 1 , is based on 4 multipliers and 2 adders needed to calculate the real and imaginary part of equation (4). This is normally referred as a full butterfly circuit or algorithm. The coefficient C I , C Q  can be fetched from a pre-calculated table, not shown in the figure, for example. 
       FIG. 2  illustrates a variant of the full butterfly circuit, or Half Butterfly, which adopts a pipeline approach. The multiplexer  201  and the de-multiplexer  206  are used to compute separately the real part  207  I·C I −Q·C Q  and the imaginary part  208  I·C Q +Q·C I  in the arithmetic block  205  that has only two multipliers  202 ,  203  and one adder  204 . The results are stored in complex register  209 . 
     Conventional DFT implementations that make use of real general purpose multipliers are sub-optimal because they do not take into considerations the fact that the possible twiddle factors W N   kn  and the C I , C Q  coefficients can not take any value, but are necessarily included in a predefined finite set. The twiddle factors W N   kn  are the N th  order complex roots of the unit and, therefore, their real and imaginary are parts have certain symmetries that are not exploited in conventional DFT implementations. 
     Preferably, the circuit of the present invention uses a set of multiplier-by-constant modules whose number is equal to the number K of different possible coefficients, disregarding their signs, plus 1 (the extra module is needed to take into account roots of the unity of the form±√{square root over (2)}/2±j·√{square root over (2)}/2, that have the same coefficient, in absolute value, for the real and imaginary part). 
     A simple formula for an approximate computation of the frequency span of each DFT is given by: 
                     Δ   ⁢           ⁢   f     =       1     T   s       ⁢     (       M   N     +     1     N   eff         )               (   6   )               
where T S  is the sample period, N is the length of the input signal for the DFT (expressed in number of points), M is the number of calculated DFT lines (with M≦N), and N eff  is the effective DFT length (how many points are effectively integrated).
 
     The number K of different coefficients depends on the number N of the DFT points. Supposing that N is an integer multiple of 8 and exploiting the symmetry of the N th  order complex roots of the unit, the number K of different coefficients is equal to N/4+1. This formula gives for N=8, 16, 24, 32, 40, 48 the following values of K=3, 5, 7, 9, 11, 15. It is straightforward to extend this to a DFT engine supporting other values of N, or any finite choice of values of K. If for instance N is intended to be configurable in the range of values {8, 16, 24, 32}, the number of needed coefficients K is (K(24)+K(32)−K C ), being K C  the number of common coefficients between N=24 and N=32 (those having complex argument 0° or 45°. 
     An embodiment of the invention will now be described with reference to  FIG. 3 . The Half butterfly arithmetic block  205 , is functionally equivalent to the block designated under the same number in  FIG. 2 . Its structure is however different, as it will be explained in the following. 
     A control logic unit  303  drives the operation of the DFT engine. It makes use of incoming configuration parameters  301  to select the data to be processed  201  from an input buffer  302 . Configuration parameters  301  and input data in buffer  302  can have various sources. In the case of a GNSS receiver, for example, configuration may be provided by a non represented CPU, and that the input data could be generated by a correlation unit. Other configuration are however possible. The control logic unit  303  also controls the selection of the current coefficients pair  304  and determines if the partial terms should be added or subtracted in the half butterfly unit  205 , depending whether the real or imaginary part of equation (4) is being processed. 
     The Half Butterfly  205  calculates first the real part  207  (I·C I −Q·C Q ) and the demultiplexer  206  drives the result to the I part of the complex result register  209 . Afterwards, the Half Butterfly calculates the imaginary part  208  (I·C Q +Q·C I ) and stores the result in the Q part of the complex result register  209 . 
     The complex result  209  is optionally multiplied by a constant factor by scaling unit  305  and accumulated, together with previous DFT terms in the accumulator comprising the adder  306  and memory  308 , for example a RAM. 
     The DFT of the invention includes preferably a saturation stage  307 , to detect overflow in the accumulation RAM  308 . If saturation occurs a flag  310  is activated. The DFT control logic  303  or an external control unit implement an algorithm to calculate a proper scaling factor to avoid saturation. 
     The DFT processing unit of the present invention is preferably arranged to disable dynamically any unused DFT line as indicated by configuration parameters  301 . The input data selector  201  can be stuck to ‘0’, reducing the toggling of the Half Butterfly  205  and scaler  305  for the discarded DFT line. Additionally, the control logic  303  can be used to update the address of the RAM  308 , but avoid the unnecessary and power consuming memory read/write accesses. 
     According to a preferred variant of the invention, the DFT processor includes a Half Butterfly  205  having the structure illustrated in  FIG. 4 , which calculates equation (4) in 2 steps. First the real part (I·C I −Q·C Q ) is calculated then the imaginary one (I·C Q +Q·C I ). 
     The structure of  FIG. 4  takes advantage of the fact that, with the sole exception of the diagonal coefficient C 45 , having a complex argument multiple of 45° it is always C i ≠C q ; the set of coefficients can be split into two separate groups of even constant coefficients {C 0 , C 2 , . . . , C k-1 } and odd constant coefficients {C 1 , C 3 , . . . , C k-2 }. The subdivision is done in a way that it never happens that two coefficient for the same group are needed to calculate the Half Butterfly terms. Block  403  contains a plurality of multipliers specifically arranged to multiply by one of the constant even coefficients {C 0 , C 2 , . . . , C k-1 }, whereas block  406  contains a plurality of multipliers specifically arranged to multiply by one of the constant odd coefficients {C 1 , C 3 , . . . , C k-2 }. Both blocks contain a multiplier by the C 45  coefficient. 
     This feature of the invention is exemplified by  FIG. 9  that illustrates the position of the twiddle factors W 32   kn  appearing in equation (4) for the special case N=32. The W 32   kn  are distributed along the unity circle in the complex plane, and are symmetrically placed about the 45° dashed line. Neglecting the sign, that can be computed trivially, the coefficient C Q  and C I  must necessarily take one of the values C 0 ,-C 8  shown. With the sole exception of the 45° twiddle factor, all the terms of equation (4) involve one coefficient from the finite set of constant numbers {C 0 , C 2 , C 4 , C 6 , C 8 } and another coefficient from the finite set of constant numbers {C 1 , C 3 , C 5 , C 7 , C 8 }. The structure of the Half Butterfly unit  205  is designed to take advantage of this symmetry. 
     Reverting now to  FIG. 4 , input samples I, Q and the K coefficients are preferably coded as sign-magnitude. Eventually additional conversion logic can be added in front of the I, Q data with a small overhead in terms of area if the incoming data are not in sign-magnitude format. Magnitudes of the data and of the selected coefficients are multiplied to obtain the magnitude of the product in blocks  403  and  406 . Due to the fact that the multipliers have to deal with absolute magnitudes, internal toggling and dynamic power consumption are much reduced with respect to traditional multipliers arranged to tread numbers in two&#39;s complement format. 
     The multiplier by constant values contained in banks  403  and  406  are preferably implemented in integer arithmetic and are highly optimized. Trivial coefficients having 0° or 90° complex argument are implemented as shift and truncation operation to further minimize area and power. If appropriate, some coefficients may deviate slightly from the theoretical value, in order to simplify the structure and reduce the power consumption of the multiplier. 
     In the presented example the output of the Half Butterfly is in two&#39;s complement format. Since the products calculated in  403  and  406  are encoded in sign-magnitude format, conversion blocks  409  carry out the necessary conversion, before the data are combined by adder  410 . The configuration port  407  selects if the Half Butterfly is calculating the real (I·C I −Q·C Q ) or imaginary part (I·C Q +Q·C I ) of the complex DFT result  209 . It is also necessary to know the sign of the data and the sign of the coefficients before converting the data to two&#39;s complement format. It is also possible to bypass the conversion block  409  to have a sign-magnitude coding of the outputs. The output of the adder  410  is the real or the imaginary part of one of the DFT of equation (4). According to a variant, the values I, Q are represented as unsigned integers, and the multipliers in banks  403  and  406  operate in unsigned mode. The sign of the result, computed separately, is set by acting on the two&#39;s complement units  408  and  409 . 
     Each of the multiplier-by-constant modules contained in banks  403  and  406  is dynamically activated by the configuration bus  402  and  405  from the control unit  303 . Thereby only the part of the circuitry really needed for the current DFT line calculation is active at any given moment, and the multipliers in banks  403  and  406  are in an inactive quiescent state for most of the time. Preferably the order of sum of terms in equation (4) can be rearranged (scrambling). 
       FIG. 5  shows a possible structure for the even coefficient bank  403 . The input signal  501  is the result of selecting between I and Q performed by multiplexer  401  ( FIG. 4 ). This signal  501  is common to all the (K+1)/2 multiplier-by-constant modules  503 . The possibility of activating only the part of the circuitry really needed for the current DFT line calculation is performed using the control signals from  402  and listed as (Sel_I_ 0 , Sel_Q_ 0 , . . . , Sel_I_K−1, Sel_Q_K−1). 
     A similar approach is shown in  FIG. 6  for the odd coefficient bank. The input signal  601  is the result of selecting between I and Q performed by multiplexer  404 . This signal  601  is common to all the (K+1)/2 multiplier-by-constant modules within  406 . The possibility of activating only the part of the circuitry really needed for the current DFT line calculation is performed using the control signals from  405  and listed as (Sel_I_ 1 , Sel_Q_ 1 , . . . , Sel_I_K−2, Sel_Q_K−2, Sel_I — 45°, Sel_Q — 45°. A multiplication unit  605  with the coefficient for 45° is also available in  403 , but has to be duplicated in  406 . What coefficient from  403  is also marked as  605  within  406  depends on the coding of coefficients adopted. 
     Only 2 coefficients are needed hence only 2 multiply-by-constant units are activated simultaneously (one from each coefficient bank  403  and  406 ). The other multiplier-by-constant blocks have their inputs tied to 0 (no consumption due to combinatorial logic toggling). 
     This approach reduces the toggling activity of about 30% if compared to a standard multiplier approach. Moreover the architecture is totally combinatorial and no pipeline stages are present inside it. Typically a pipeline in digital circuitry is implemented using simple flip-flop based registers that are not optimized in term of area and power consumption. Avoiding them area and power are minimized. 
     The order of the DFT lines being calculated is managed inside the control logic blocks  303  and  304 . The DFT lines can be calculated with a programmable order so that the post-processing computation load for a CPU is reduced. 
     With reference to  FIGS. 5 and 6  and supposing that C I =C 0  and C Q =C 1  the calculation of equation (4) can be performed in 2 steps described below.
         Step 1: calculation of real part of equation (4): Sel_I_ 0 =1, Sel_Q_ 1 =1, Sel_I_x=0 if x≠0, Sel_Q_x=0 if x≠1, Sel_I=1, Sel_Q=1, DFT_Re=1. Then A=I·C 0 , B=−Q·C 1      Step 2: calculation of imaginary part of equation (4): Sel_I_ 1 =1, Sel_Q_ 0 =1, Sel_I_x=0 if x≠1, Sel_Q_x=0 if x≠0, Sel_I=0, Sel_Q=0, DFT_Re=0. Then D=Q·C 0 , C=I·C 1  
 
GNSS Receiver Embodiment
       

     DFT algorithms are generally known in the art and described in the technical literature. In the following only the aspects specific to GNSS implementation will be discussed with reference to  FIGS. 7 and 8 . 
     With reference to  FIG. 7 , the acquisition and tracking of a GPS space vehicle (SV) requires the determination of the frequency/code bin  705 . For this purpose specific resources are needed to determine the code phase offset bin  703  and the Doppler bin  704 . In the case that there is no estimation of one of more of the above cited parameters, a full search over the entire frequency/code search space should be performed. 
     A serial search approach, where a frequencies sweep over all possible Doppler of the incoming GPS signal and a code phases sweep over 1023 possible values for a GPS PRN (Pseudo Random Noise) code, is a widely used method for the acquisition step in a GNSS system. 
     A parallel frequency approach may be used to speed up the acquisition process. The receiver architecture illustrated in  FIG. 8  carries out a parallel search in the frequency domain calculating the DFT  806  of the signal generated from the correlation  805  between the processed GNSS signals  802  and a locally generated replica  804  of the PRN for a given SV. 
     The line-of-sight velocity of the satellite referred to the receiver cause a Doppler effect in the order of +/−10 KHz. A step frequency step of 150 Hz is the minimum required for a low level GNSS signal scenario. 
     The gain of the correlated signal  805  has the format sinc(x):=sin(x)·x −1  where x=πfT. Applying this sinc envelope to the DFT transfer function it becomes evident that all the DFT lines will be affected by an amplitude loss with the exception of the centre frequency line f (that correspond to the Doppler frequency) under the condition that the PRN code is perfectly aligned. Otherwise no peak is present. 
     The aforementioned properties justify the use of a reduced number M of DFT lines from N data samples out of the correlator  805 . 
     The present invention further concerns a DFT processor for a reduced number of spectrum lines to reduce hardware complexity and power consumption. 
     By careful application of appropriate design constraints specific to the SV navigation and analysis of the DFT algorithm an optimised hardware architecture can be realised for embedding frequency-domain analysis efficiently into a GNSS chipset.