Abstract:
A method for solving a system of N linear equations in N unknown variables. The method comprising:  
     (a) storing an estimate value for each unknown variable;  
     (b) initialising each estimate value to a predetermined value;  
     (c) for each estimate value:  
     (i) determining whether a respective predetermined condition is satisfied; and  
     (ii) updating the estimate if and only if the respective predetermined condition is satisfied; and  
     repeating step (c) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.

Description:
CROSS REFERENCE TO RELATED APPLICATIONS  
       [0001]    This application is a continuation in part of International Application PCT/GB03/001568, filed Apr. 10, 2003, which claims priority to Great Britain Application No. GB 0208329.3, filed Apr. 11, 2002, the contents of each of which are incorporated herein by reference. 
     
    
     
       FIELD OF INVENTION  
         [0002]    The present invention relates to systems and methods for solving systems of linear equations.  
         BACKGROUND OF INVENTION  
         [0003]    Systems of linear equations occur frequently in many branches of science and engineering. Effective methods are needed for solving such equations. It is desirable that systems of linear equations are solved as quickly as possible.  
           [0004]    A system of linear equations typically comprises N equations in N unknown variables. For example, where N=3 an example system of equations is set out below:  
           15 x+ 5 y− 2 z =15  
           5 x+ 11 y+ 4 z= 47  
           −2 x+ 4 y+ 9 z= 51  
           [0005]    In this case, it is necessary to find values of x, y, and z which satisfy all three equations. Many methods exist for finding such values of x, y and z and in this case it can be seen that x=1, y=2 and z=5 is the unique solution to the system of equations.  
           [0006]    In general terms, existing methods for solving linear equations can be categorised as either direct methods, or iterative methods. Direct methods attempt to produce an exact solution by using a finite number of operations. Direct methods however suffer from a problem in that the number of operations required is often large which makes the method slow. Furthermore, some implementations of such methods are sensitive to truncation errors.  
           [0007]    Iterative methods solve a system of linear equations by repeated refinements of an initial approximation until a result is obtained which is acceptably close to the accurate solution. Each iteration is based on the result of the previous iteration, and, at least in theory, each iteration should improve the previous approximation. Generally, iterative methods produce an approximate solution of the desired accuracy by yielding a sequence of solutions which converge to the exact solution as the number of iterations tends to infinity.  
           [0008]    It will be appreciated, that systems of linear equations are often solved using a computer apparatus. Often, equations will be solved by executing appropriate program code on a microprocessor. In general terms, a microprocessor can represent decimal numbers in fixed point or floating point form.  
           [0009]    In fixed point form, a binary number of P bits comprises a first part of length q and a second part of length r. The first part represents the whole number part, and the second part represents the fractional part. In general terms, arithmetic operations can be relatively efficiently implemented for fixed point numbers. However, fixed point numbers suffer from problems of flexibility given that the position of the decimal point is fixed and therefore the range of numbers which can be accurately represented is relatively small given that overflow and round off errors regularly occur.  
           [0010]    Floating point numbers again in general terms comprise two parts. A first part, known as the exponent, and a second part known as the mantissa. The mantissa represents the binary number, and the exponent is used to determine the position of the decimal point within that number.  
           [0011]    For example considering an eight bit value 00101011 where the first bit represents sign, the subsequent three bits represent the exponent, and the final four bits represent the mantissa, this value is analysed as follows. The first bit represents sign, and is interpreted such that the number is positive if the sign bit is equal to ‘0’ and negative if the sign bit is equal to ‘1’. In this case, given that the first bit is ‘0’, the number is determined to be positive. The mantissa (1011) is written out and a decimal point is placed to its left side as follows:  
           [0012]    0.1011  
           [0013]    The exponent is then interpreted as an integer, that is 010 is interpreted as the integer value 2 (given that the first bit of the exponent field is a sign bit which determines the direction in which the decimal point should move). In this case, the decimal point is moved to the right two bits to give:  
           [0014]    10.11  
           [0015]    Which represents a value of two and three quarters.  
           [0016]    Although floating point numbers give considerable benefits in terms of their flexibility, arithmetic operations involving floating point numbers are inherently slower than corresponding operations on fixed point numbers. Therefore, where possible, the benefits of speed associated with fixed point numbers should be exploited.  
           [0017]    When considering the implementation of an algorithm in hardware, its efficiency is a prime concern. Many algorithms for the solution of linear equations involve computationally expensive division and/or multiplication operations. These operations should, where possible be avoided, although this is often not possible with known methods for solving linear equations.  
           [0018]    Many systems of linear equations have sparse solutions. In such cases the number of iterations required to solve the system of equations should be relatively low. However, this does not occur with some known methods.  
           [0019]    In summary there is a need for a more efficient algorithm which can be used to solve systems of linear equations.  
           [0020]    It is an object of the present invention to obviate or mitigate at least one of the problems identified above.  
         SUMMARY OF THE INVENTION  
         [0021]    According to a first aspect of the present invention, there is provided a method for solving a system of N linear equations in N unknown variables, the method comprising:  
           [0022]    (a) storing an estimate value for each unknown variable;  
           [0023]    (b) initialising each estimate value to a predetermined value;  
           [0024]    (c) for each estimate value:  
           [0025]    (i) determining whether a respective predetermined condition is satisfied; and  
           [0026]    (ii) updating the estimate if and only if the respective predetermined condition is satisfied; and  
           [0027]    (d) repeating step (c) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0028]    Thus, the invention provides a method for solving linear equations in which estimates for solutions of the equations are updated only if a predetermined condition is satisfied. The predetermined condition is preferably related to convergence of the method. Therefore such an approach offers considerable benefits in terms of efficiency, given that updates are only carried out when such updates are likely to accelerate convergence.  
           [0029]    According to a second aspect of the present invention, there is provided a method for solving a system of N linear equations in N unknown variables, the method comprising:  
           [0030]    (a) storing an estimate value for each unknown variable;  
           [0031]    (b) initialising each estimate value to a predetermined value;  
           [0032]    (c) attempting to update each estimate value using a scalar value d;  
           [0033]    (d) updating the scalar value if no updates are made in step (c); and  
           [0034]    (e) repeating step (c) and step (d) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0035]    By updating the scalar value in accordance with the second aspect of the present invention, it has been discovered that benefits of efficiency are obtained.  
           [0036]    According to a third aspect of the present invention, there is provided a method for solving a system of N linear equations in N unknown variables of the form:  
           Rh=β 
           [0037]    the method comprising:  
           [0038]    generating a quadratic function of the form:  
           J        (   h   )       =         ∑     m   =   1     N            ∑     n   =   1     N            R        (     m   ,   n     )            h        (   m   )            h        (   n   )             -     2          ∑     n   =   1     N            β        (   n   )            h        (   n   )                 ;                         
 
           [0039]    and  
           [0040]    minimising said function using co-ordinate descent optimisation; wherein R is a coefficient matrix of the system of linear equations; h is a vector of the N unknown variables; β is a vector containing the value of the right hand side of each equation; R(m,n) is an element of the matrix R; h(m) is the mth element of the matrix h; and β (n) is the nth element of the vector β.  
           [0041]    The present inventors have discovered that solving a system of linear equations by minimising a quadratic function using co-ordinate descent optimisation offers considerable and surprising efficiency benefits.  
           [0042]    According to a fourth aspect of the present invention, there is provided a computer processor configured to solve a system of N linear equations in N unknown variables, comprising:  
           [0043]    storage means for storing an estimate value for each unknown variable;  
           [0044]    storage means for storing coefficients of each unknown variable in each equation;  
           [0045]    storage means for storing a scalar value d;  
           [0046]    initialising means for initialising each estimate value;  
           [0047]    computing means configured to process each estimate value by determining whether a respective predetermined condition is satisfied, and to update the estimate if and only if the respective predetermined condition is satisfied, said computing means being configured to repeatedly process each estimate value until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0048]    According to a fifth aspect of the present invention, there is provided a computer processor configured to solve a system of N linear equations in N unknown variables, comprising:  
           [0049]    storage means for storing an estimate value for each unknown variable;  
           [0050]    storage means for storing coefficients of each unknown variable in each equation;  
           [0051]    storage means for storing a scalar value d;  
           [0052]    initialising means for initialising each estimate value;  
           [0053]    computing means configured to:  
           [0054]    (a) attempt to update each estimate value using a scalar value d,  
           [0055]    (b) update the scalar value d if no updates are made in step (a); and  
           [0056]    (c) repeat step (a) and step (b) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0057]    According to a sixth aspect of the present invention, there is provided a multiuser receiver device for obtaining data transmitted by a plurality users, the device comprising:  
           [0058]    a plurality of filters, each filter being arranged to filter out a spreading code used by a respective user;  
           [0059]    equation solving means to find a solution h of a system of linear equations of the form Rh=β where R is the cross correlation of the spreading codes used by the plurality of users, and β is a vector containing the filter output signals; and  
           [0060]    means for obtaining the transmitted data using a solution provided by the equation solving means;  
           [0061]    wherein the equation solving means:  
           [0062]    (a) stores an estimate value for each value of the solution h;  
           [0063]    (b) initialises each estimate value to a predetermined value;  
           [0064]    (c) for each estimate value:  
           [0065]    (i) determines whether a respective predetermined condition is satisfied; and  
           [0066]    (ii) updates the estimate if and only if the respective predetermined condition is satisfied; and  
           [0067]    (d) repeats step (c) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0068]    According to a seventh aspect of the present invention, there is provided a multiuser receiver device for obtaining data transmitted by a plurality of users, the device comprising:  
           [0069]    a plurality of filters, each filter being arranged to filter out a spreading code used by a respective user;  
           [0070]    equation solving means to find a solution h of a system of linear equations of the form Rh=β where R is the cross correlation of the spreading codes used by the plurality of users, and β is a vector containing the filter output signals; and  
           [0071]    means to obtain the transmitted data using a solution provided by the equation solving means;  
           [0072]    wherein the equation solving means:  
           [0073]    (a) stores an estimate value for each unknown variable;  
           [0074]    (b) initialises each estimate value to a predetermined value;  
           [0075]    (c) attempts to update each estimate value using a scalar value d;  
           [0076]    (d) updates the scalar value if no updates are made in step (c); and  
           [0077]    (e) repeats step (c) and step (d) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0078]    According to an eighth aspect of the present invention, there is provided a method for generating filter coefficients for use in an echo cancellation apparatus, the method comprising:  
           [0079]    (a) generating a cross correlation matrix R containing the cross correlation of first and second signals and;  
           [0080]    (b) generating an auto correlation vector β containing an autocorrelation of the first signal; and  
           [0081]    (c) determining a vector h for which Rh=β, said vector h containing the said filter coefficients;  
           [0082]    wherein the vector h is determined by solving the system of equations Rh=β by:  
           [0083]    (d) storing an estimate value for each element of the vector h;  
           [0084]    (e) initialising each estimate value to a predetermined value;  
           [0085]    (f) for each estimate value:  
           [0086]    (i) determining whether a respective predetermined condition is satisfied; and  
           [0087]    (ii) updating the estimate if and only if the respective predetermined condition is satisfied; and  
           [0088]    (g) repeating step (f) until each estimate value is sufficiently close to an accurate value of the respective unknown variable.  
           [0089]    According to a ninth aspect of the present invention, there is provided a method for generating filter coefficients for use in an echo cancellation apparatus, the method comprising:  
           [0090]    (a) generating a cross correlation matrix R containing the cross correlation of first and second signals;  
           [0091]    (b) generating an auto correlation vector β containing an autocorrelation of the first signal; and  
           [0092]    (c) determining a vector h for which Rh=β, said vector h containing the said filter coefficients;  
           [0093]    wherein the vector h is determined by solving the system of equations Rh=β by:  
           [0094]    (d) storing an estimate value for each unknown variable;  
           [0095]    (e) initialising each estimate value to a predetermined value;  
           [0096]    (f) attempting to update each estimate value using a scalar value d;  
           [0097]    (g) updating the scalar value if no updates are made in step (f); and  
           [0098]    (h) repeating step (g) and step (h) until each estimate value is sufficiently close to an accurate value of the respective unknown variable. 
       
    
    
     BRIEF DESCRIPTION OF DRAWINGS  
       [0099]    Embodiments of the present invention will now be described, by way of example, with reference to the accompanying drawings, in which:  
         [0100]    [0100]FIG. 1 is a flow chart of an algorithm for solving linear equations in accordance with the present invention;  
         [0101]    [0101]FIG. 2 is a graph showing how the value of the solution vector changes as the algorithm of FIG. 1 solves a system of equations;  
         [0102]    [0102]FIG. 3 is a graph showing how the error of the solution vector varies as the algorithm of FIG. 1 solves the system of equations;  
         [0103]    [0103]FIG. 4 is a graph showing how the number of passes through the system equations varies in dependence upon the bit number of the solution vector elements being considered in the algorithm of FIG. 1;  
         [0104]    [0104]FIG. 5 is a graph showing the number of updates to the solution vector carried out for each bit as the algorithm of FIG. 1 solves the system of equations;  
         [0105]    [0105]FIG. 6 is a flow chart showing a variant to the algorithm of FIG. 1;  
         [0106]    [0106]FIG. 7 is a MATLAB code fragment implementing the algorithm of FIG. 6;  
         [0107]    [0107]FIG. 8 is a flow chart showing a variant to the algorithm of FIG. 6, where values of the solution vector are constrained between upper and lower bounds.  
         [0108]    [0108]FIG. 9 is flow chart showing how the algorithm of FIG. 1 can be adapted to solve equations having complex coefficients and complex solutions;  
         [0109]    [0109]FIG. 10 is a flow chart showing a method for determining the next course of action in a step of FIG. 9;  
         [0110]    [0110]FIG. 11 is a flow chart showing a variant to the algorithm of FIG. 10;  
         [0111]    [0111]FIG. 12 is a MATLAB code fragment implementing the algorithm of FIG. 11;  
         [0112]    [0112]FIG. 13 is a schematic illustration of a device configured to solve linear equations in accordance with the invention;  
         [0113]    [0113]FIG. 14 is a schematic illustration of the equation solving microprocessor of FIG. 13;  
         [0114]    [0114]FIG. 15 is a schematic illustration showing block R of FIG. 14 in further detail;  
         [0115]    [0115]FIG. 16 is a schematic illustration showing block h of FIG. 14 in further detail;  
         [0116]    [0116]FIG. 17 is a schematic illustration showing how the block h of FIG. 16 is modified when equations having complex valued solutions are to be solved;  
         [0117]    [0117]FIG. 18 is a schematic illustration showing block Q of FIG. 14 in further detail;  
         [0118]    [0118]FIG. 19 is a schematic illustration showing the minimisation block of FIG. 14 in further detail;  
         [0119]    [0119]FIG. 19 a  is a MATLAB code fragment showing a variant to the algorithm of FIG. 7;  
         [0120]    [0120]FIG. 20 is a schematic illustration of a CDMA receiver device in which algorithms of the present invention may be applied;  
         [0121]    [0121]FIG. 21 is a schematic illustration of an echo cancellation apparatus in which the algorithms of the present invention may be applied; and  
         [0122]    [0122]FIG. 22 is a schematic illustration showing part of the apparatus of FIG. 21 in further detail. 
     
    
     DESCRIPTION OF PREFERRED EMBODIMENTS  
       [0123]    A method for a system of solving linear equations is now described. A system of linear equations can be expressed in the form:  
         Rh=β  (1)  
         [0124]    where: R is a coefficient matrix of the system of equations;  
         [0125]    h is a vector of the unknown variables; and  
         [0126]    β is a vector containing the value of the right hand side of each equations  
         [0127]    For example, the system of equations (2):  
         15 x+ 5 y− 2 z= 15  
         5 x+ 11 y+ 4 z =47  
         −2 x+ 4 y+ 9 z =51  (2)  
         [0128]    can be expressed in the form of equation (1) where:  
             R   =         [         15       5         -   2             5       11       4             -   2         4       9         ]                   h     =         [         x           y           z         ]                   β     =     [         15           47           51         ]                 (   3   )                               
 
         [0129]    To solve the system of equations, it is necessary to find values for x, y, and z of h which satisfy each of the three equations.  
         [0130]    In operation, algorithm uses the matrix R and the vectors h and β as set out above, together with an auxiliary vector Q. The vector h is initialised to a predetermined initial value (see below) and updated as the algorithm proceeds until its elements represent the solution of the equations.  
         [0131]    For a system of N equations in N unknown variables, the vector h has length N and the matrix R is of size N×N.  
         [0132]    Referring to FIG. 1, at step S 1  the vectors h and Q are initialised. The vector h is initialised such that all its elements are set to ‘0’. The vector Q is initialised to contain the negative of the equivalent position of β. That is:  
         Q=−β  (4)  
         [0133]    Therefore, when working with system of equations (2), Q is initialised in accordance with equation (5):  
             Q   =     [           -   15               -   47               -   51           ]             (   5   )                               
 
         [0134]    The algorithm maintains three counter variables p, m and it, a parameter N which represents the number of elements in the solution vector (and also the number of equations), a parameter M which represents the number of bits used to represent each element of the solution vector h, a parameter Nit which represents the maximum number of iterations through which the algorithm can pass for a particular value of m, a variable Flag which is used to indicate whether or not the solution vector has been updated, and a constant H, the purpose of which is described below.  
         [0135]    Some of these variables are initialised at step S 2  and step S 3  of FIG. 1. p, m and it are all initialised to zero. N, M, and Nit are set to the values described above which can either be chosen by a user or hard coded into the algorithm. Selection and initialisation of H is described below.  
         [0136]    Operation of the algorithm can be summarised as follows. Each bit m of all elements p of the solution vector h is considered in turn. As described below, for each bit an element of the vector Q is compared with various conditions and the result of this comparison determines whether or not further processing is applicable. This further processing comprises an appropriate update of the element p of the solution vector h corresponding to the element being considered and updates of all elements of the auxiliary vector Q.  
         [0137]    When it is determined that further processing for that element is not appropriate (for the current bit), the next element is considered. When each element has been considered for that particular bit, all elements of the solution vector are considered for the next bit in turn, and updated appropriately. This process continues until all elements have been considered for all bits. If the total number of iterations for any one bit reaches a predetermined limited the algorithm again ends. The algorithm is described in further detail below.  
         [0138]    At step S 3 , the value of m is incremented to ‘1’. Thus, the algorithm is now considering the first bit of each element in the solution vector h. it is set to 0 to indicate that no iterations have yet taken place for the current value of m.  
         [0139]    At Step S 4 , a step size parameter d is calculated according to the equation:  
           d= 2 −m   H   (6)  
         [0140]    where m and H are the parameters described above.  
         [0141]    H is a value greater than or equal to the magnitude of the maximum value which is expected for any value of the solution vector. That is, the algorithm considers only solutions lying between −H and +H.  
         [0142]    As will be described below, setting d in accordance with equation (6) allows each bit of each value of the solution vector h to be considered in turn.  
         [0143]    At Step S 5  of FIG. 1, the variable it is incremented, and the variable Flag is set to ‘0’. p (the current element of the solution vector under consideration) is incremented at Step S 6 .  
         [0144]    Having performed the necessary incrementation and initialisation, the algorithm can begin processing elements of the matrix and vectors, in an attempt to solve the equations.  
         [0145]    At step S 7 , the following operation is performed:  
               arg   =     arg                 min                   {       Q        (   p   )       ,     -     Q        (   p   )         ,       -     R        (     p   ,   p     )         ·     d   2         }              
        where           (   7   )                 arg                 min     =     {           1   ,     if        {       Q        (   p   )       &lt;       -     Q        (   p   )         ⋀     Q        (   p   )         &lt;       -     R        (     p   ,   p     )         ·     d   2         }                   2   ,     if        {       -     Q        (   p   )         &lt;       Q        (   p   )       ⋀     -     Q        (   p   )           &lt;       -     R        (     p   ,   p     )         ·     d   2         }                   3   ,     if        {         -     R        (     p   ,   p     )         ·     d   2       ≤         Q        (   p   )       ⋀     -     R        (     p   ,   p     )           ·     d   2       ≤     -     Q        (   p   )           }                         (   8   )                               
 
         [0146]    The value of arg is assessed at the decision block of step S 8 .  
         [0147]    If arg=1, the element of the solution vector under consideration, that is h(p) is set according to equation (9) at step S 9 .  
           h ( p )= h ( p )+ d   (9)  
         [0148]    The auxiliary vector Q is then updated such that all values of Q are set according to equation (10) at step S 11 .  
           Q ( r )= Q ( r )+ dR ( p,r ),∀ r: 1≦ r≦N   (10)  
         [0149]    If arg=2, the element of the solution vector under consideration, that is h(p), is set according to equation (11) at step S 11 .  
           h ( p )= h ( p )− d   (11)  
         [0150]    The auxiliary vector Q is then updated such that all values of Q are set according to equation (12) at step S 10 .  
           Q ( r )= Q ( r )− dR ( p,r ),∀ r: 1≦ r≦N   (12)  
         [0151]    If arg=1 or if arg=2, Flag is set to ‘1’ at step S 13 .  
         [0152]    If arg=3, no update is made to any element of the solution vector h or the auxiliary vector Q, and Flag is not updated.  
         [0153]    Having made the updates set out above, a decision block at step S 14  checks the condition of equation (13);  
         p=N  (13)  
         [0154]    If p is not equal to N (i.e. all elements of the solution vector h have not yet been considered), control returns to step S 6  and p is incremented. This process continues until all entries in the solution vector h have been considered, and h and Q are updated in the manner set out above.  
         [0155]    If p is equal to N (step S 14 ), a check is made to determine the value of Flag (step S 15 ).  
         [0156]    Flag is initially set to ‘0’ at step S 5 , and only updated (to be equal to ‘1’) if entries of the solution and auxiliary vectors are amended by steps S 9  and S 10 , or steps S 11  and S 12 . Thus, if Flag=1, it can be deduced that a change was made to at least one element of h (i.e. one h(p) value) and all values of Q, for the current iteration it. Therefore, assuming that the total number of iterations it has not exceeded the limit set by Nit (checked at S 16 ), p is reset to ‘0’ at step S 17 , control returns to step S 5 , and the current bit is again processed for each element p of the solution vector h. This is because further processing of each element of h for the current value of m may result in further updates. If the total number of iterations has reached the limit set by Nit (step S 16 ), the algorithm exits.  
         [0157]    If it is the case that Flag=0 (step  15 ), it can be deduced that no updates have been made to any elements of the solution vector h or the auxiliary vector Q for any value of p (that is any element of the solution vector h). Given that further iterations of steps S 5  to S 15  will result in no changes to the elements of h (given that neither d, h nor Q have changed), a check is made to determine whether or not all bits m have been considered (step S 18 ), by comparing the current value of m with the total number of bits M. If it is the case that m=M, i.e. all bits have been considered, there is no work for the algorithm to do, and the algorithm exits at step S 19 .  
         [0158]    If it is the case that all bits have not been considered, p is reset to 0 at step S 20 , and control returns to step S 3 , and the algorithm processes the next bit of the solution vector entries.  
         [0159]    In the preceding discussion, it has been explained that entries of the solution vector h are processed for each bit of the solution vector entries, starting with the most significant bit. However, it can be seen from the preceding discussion, that at all steps the entire value of an element of h is used for update. However bit wise processing is in fact achieved because following each increment of m (step S 3 ) a new value of d is created at step S 4 . Given that each increment of m will result in the value of d being divided by two (given the presence of the expression 2 −m  in the expression for d), and given that d is used to update both h and Q, after an update of d the next most significant bit is then updated.  
         [0160]    It has been described that the value of H represents a value greater than or equal to the magnitude of the maximum value of the solution vector elements. In setting H it is desirable to ensure that it is a power of two. That is, H is set according to equation (14).  
         H=2 q   (14)  
         [0161]    where q is any integer (i.e. positive, negative or zero)  
         [0162]    When H is set in this way, the expression for d set out in equation (6) becomes:  
           d= 2 −m 2 q    
         d=2 q-m   (15)  
         [0163]    Thus, when H is chosen in accordance with equation (14), the value of d can be updated without multiplication or division, simply by appropriate bit shift operations. This is particularly advantageous, because microprocessors can typically carry out bit shift operations far more efficiently than multiplication or division.  
         [0164]    The application of the algorithm described with reference to FIG. 1 to the system of linear equations (2) set out above will now be described.  
         [0165]    R and β are initialised as described above:  
             R   =     [         15       5         -   2             5       11       4             -   2         4       9         ]             (   16   )               β   =     [         15           47           51         ]             (   17   )                               
 
         [0166]    h and Q are initialised in the manner described above:  
             h   =     [         0           0           0         ]             (   18   )               Q   =     [           -   15               -   47               -   51           ]             (   19   )                               
 
         [0167]    Variables are initialised as follows at step S 2  and step S 3 :  
         m=0 M=8  
         it=0 N=3  
         p=0 Nit=0  (20)  
         [0168]    H is in this case set to 16, i.e. q=4 in equation (14).  
         [0169]    m is incremented such that m=1, and it is set to ‘0’ (step S 3 ). A value of d is computed according to equation (15) and in this case:  
         d=2 4-1    
         d=8  (21)  
         [0170]    it is incremented to ‘1’ and Flag is set to ‘0’ at step S 5 . p is incremented to 1 at step S 6 . At step S 7 , the following expression is evaluated  
               arg   =     arg                 min        {       Q        (   1   )       ,     -     Q        (   1   )         ,       -     R        (     1   ,   1     )         ·     8   2         }              
          arg              =                arg                 min        {       -   15     ,   15   ,       -   15     ·     8   2         }              
          arg   =     arg                 min        {       -   15     ,   15   ,     -   60       }                 (   22   )                               
 
         [0171]    Therefore, it can be seen from equations (8) and (22) that arg=3. The decision block at step S 8  therefore passes control to step S 14 , where the condition p=N is checked. In this case p=1 and N=3, therefore p is not equal to N, and therefore control passes to step S 6 , with no changes having been made to the elements of h or Q.  
         [0172]    Step S 6  then increments p to be equal to 2 (i.e. the second element of the solution vector is being considered)  
         [0173]    Step S 7  computes:  
               arg   =     arg                 min        {       Q        (   2   )       ,     -     Q        (   2   )         ,       -     R        (     2   ,   2     )         ·     8   2         }              
          arg              =                arg                 min        {       -   47     ,   47   ,       -   11     ·     8   2         }              
          arg   =     arg                 min        {       -   47     ,   47   ,     -   44       }                 (   23   )                               
 
         [0174]    Therefore, arg=1.  
         [0175]    At step S 8  the appropriate decision outcome is chosen, and at step S 9 , h is set as follows:  
                 h          :                [         0           0           0         ]                       →   becomes                     [         0             0   +   d             0         ]       =     [         0           8           0         ]             (   24   )                               
 
         [0176]    At step S 10 , all values of Q are set by adding the current values to the values of the second row (since p=2) of R multiplied by d:  
                       Q          :                [           -   15               -   47               -   51           ]                       →   becomes                     [             -   15     +     d   ·     R        (     2   ,   1     )                       -   47     +     d   ·     R        (     2   ,   2     )                       -   51     +     d   ·     R        (     2   ,   3     )                 ]       =     [             -   15     +     8   ·   5                   -   47     +     8   ·   11                   -   51     +     8   ·   4             ]                 =     [         25           41             -   19           ]                   (   25   )                               
 
         [0177]    At step S 113 , Flag is set to ‘1’ to show that h and Q have been updated.  
         [0178]    At step S 14  p is still not equal to N (p=2, N=3) and therefore control returns to step S 6 , where p is incremented to 3.  
         [0179]    At step S 7  the following is computed:  
               arg   =     arg                 min        {       Q        (   3   )       ,     -     Q        (   3   )         ,       -     R        (     3   ,   3     )         ·     8   2         }              
          arg   =     arg                 min        {       -   19     ,   19   ,       -   9     ·     8   2         }              
          arg   =     arg                 min        {       -   19     ,   19   ,     -   36       }                 (   26   )                               
 
         [0180]    Therefore arg=3. No updates are made, and the condition of step S 14  is checked. In this instance p=N=3, and the condition is therefore true. The value of Flag is then checked at step S 15 . Flag was set to ‘1’ at step S 13  while p was set to 2, and therefore the condition of step S 15  evaluates to false.  
         [0181]    The condition of step S 16  is then checked. Given that it=1, and Nit=10, the condition of step S 16  returns False, and control passes to step S 17  where p is reset to ‘0’ and then to step S 5 , where it is incremented and Flag is reset to ‘0’.  
         [0182]    p is incremented at step S 6 , and at step S 7 , the following expression is considered  
               arg   =     arg                 min        {       Q        (   1   )       ,     -     Q        (   1   )         ,       -     R        (     1   ,   1     )         ·     8   2         }              
          arg   =     arg                 min        {     25   ,     -   25     ,       -   15     ·     8   2         }              
          arg   =     arg                 min        {     25   ,     -   25     ,     -   60       }                 (   27   )                               
 
         [0183]    It can be seen from equation (27) that arg=3. The decision block therefore passes control to step S 14 , where the condition p=N is checked. In this case p=1 and N=3, therefore p is not equal to N, and control therefore passes to step S 6 , with no changes having been made to the elements of h or Q.  
         [0184]    p is incremented to 2 at step S 6  and the following expression is considered at step S 7 :  
               arg   =     arg                 min        {       Q        (   2   )       ,     -     Q        (   2   )         ,       -     R        (     2   ,   2     )         ·     8   2         }              
          arg   =     arg                 min        {     41   ,     -   41     ,       -   11     ·     8   2         }              
          arg   =     arg                 min        {     41   ,     -   41     ,     -   44       }                 (   28   )                               
 
         [0185]    It can be seen from equation (28) that arg=3. The decision block therefore passes control to step S 14 , where the condition p=N is again checked. p is still not equal to N, and therefore control passes to step S 6 , with no changes having been made to the elements of h or Q.  
         [0186]    The following expression is then considered:  
               arg   =     arg                 min        {       Q        (   3   )       ,     -     Q        (   3   )         ,       -     R        (     3   ,   3     )         ·     8   2         }              
          arg   =     arg                 min        {       -   19     ,   19   ,       -   9     ·     8   2         }              
          arg   =     arg                 min        {       -   19     ,   19   ,     -   36       }                 (   29   )                               
 
         [0187]    Again, arg=3, and no updates are made. In this case the condition of step S 14  returns true, and the condition of step S 15  also returns true, given that no updates where made during the last pass through all elements of h, and consequently Flag is still set to ‘0’.  
         [0188]    The condition of step S 18  is then checked to determine whether all bits of the solution vector entries have been considered. In this case m=1 and M=8, therefore step S 18  returns false. p is set to ‘0’ at step S 20 , and control passes to step S 3  where m is incremented to 2 and it is reset to ‘0’.  
         [0189]    At step S 4  d is set where m is equal to 2, and therefore d=4. Note that as discussed above, d has been halved. At step S 5  it is incremented to 1 and Flag is set to ‘0’.  
         [0190]    p is incremented to 1 at step S 6 . At step S 7 , the following expression is considered:  
               arg   =     arg                 min        {       Q        (   1   )       ,     -     Q        (   1   )         ,       -     R        (     1   ,   1     )         ·     4   2         }              
          arg   =     arg                 min        {     25   ,     -   25     ,       -   15     ·     4   2         }              
          arg   =     arg                 min        {     25   ,     -   25     ,     -   30       }                 (   30   )                               
 
         [0191]    Therefore arg=3, and no update takes place. The algorithm continues as described above, and p is incremented to 2 at step S 6 . At step S 7  the following expression is evaluated:  
               arg   =     arg                 min        {       Q        (   2   )       ,     -     Q        (   2   )         ,       -     R        (     2   ,   2     )         ·     4   2         }              
          arg   =     arg                 min        {     41   ,     -   41     ,       -   11     ·     4   2         }              
          arg   =     arg                 min        {     41   ,     -   41     ,     -   22       }                 (   31   )                               
 
         [0192]    In this case arg=2. The decision block of step S 8  therefore directs control to step S 11 . Step S 11  updates h by subtracting the current value of d (d=4) from the second element of h, as shown in equation (32):  
                 h          :            [         0           8           0         ]                         →   becomes                     [         0             8   -   d             0         ]       =     [         0           4           0         ]             (   32   )                               
 
         [0193]    Step S 12  updates Q by subtracting the row of R (given that p=2) multiplied by d from the current values of Q as set out in equation (33):  
                       Q          :            [         25           41             -   19           ]                         →   becomes                     [           25   -     d   ·     R        (     2   ,   1     )                     41   -     d   ·     R        (     2   ,   2     )                       -   19     -     d   ·     R        (     2   ,   3     )                 ]       =     [           25   -     4   ·   5                 41   -     4   ·   11                   -   19     -     4   ·   4             ]                 =     [         5             -   3               -   35           ]                   (   33   )                               
 
         [0194]    Flag is set to 1 at step S 13 .  
         [0195]    Another iteration is then carried out, wherein p is set to 3 at step S 6 , and the following expression is considered at step S 7 :  
               arg   =     arg                 min        {       Q        (   3   )       ,     -     Q        (   3   )         ,       -     R        (     3   ,   3     )         ·     4   2         }              
          arg   =     arg                 min        {       -   35     ,   35   ,       -   9     ·     4   2         }              
          arg   =     arg                 min        {       -   35     ,   35   ,     -   18       }                 (   34   )                               
 
         [0196]    Therefore, arg is set to 1, the decision block of step S 8  directs control to step S 9 , and steps S 9  and S 10  set h and Q as set out below:  
                 h          :            [         0           4           0         ]                         →   becomes                     [         0           4             0   +   d           ]       =     [         0           4           4         ]             (   35   )                         Q          :            [         5             -   3               -   35           ]                         →   becomes                     [           5   +     d   ·     R        (     3   ,   1     )                       -   3     +     d   ·     R        (     3   ,   2     )                       -   35     +     d   ·     R        (     3   ,   3     )                 ]       =     [           5   +     4   ·     -   2                     -   3     +     4   ·   4                   -   35     +     4   ·   9             ]                 =     [           -   3             13           1         ]                   (   36   )                               
 
         [0197]    Given that p=N and Flag=1, p is reset (step S 17 ) and control passes to S 5 . For each of the three values ofp steps S 6  to. S 14  are executed. On each occasion, arg=3 and no updates take place. The individual expressions considered by step S 7  during each iteration are not set out in full, although they can be readily deduced from the information presented above.  
         [0198]    Given that Flag=0, the algorithm continues for the next value of m.  
         [0199]    Execution of the algorithm then continues in the manner outlined above, for each bit of the solution vector elements in turn.  
         [0200]    The values of h after each iteration through the solution vector elements are set out below. It should be noted that iteration number referred to here is equivalent to a cumulative iteration count instead of the bit by bit iteration count it described above.  
             [           Iteration        :           0       1       2       3       4       5       6       7       8       9       10             1      st                 element        :           0       0       0       0       0       0       0       1       1       1       1             2      nd                 element        :           0       8       8       4       4       2       2       2       2       2       2             3      rd                 element        :           0       0       0       4       4       4       4       5       5       5       5         ]                           
 
         [0201]    [0201]FIG. 2 is a graph showing the value of each solution vector element after each iteration. A first line  1  represents changes to the first solution vector element, h(1), a second line  2  represents changes to the second solution vector element h(2), and a third line  3  represents changes to the third solution vector element h(3).  
         [0202]    Solving the set of equations set out at (2) above in a conventional way yields  
             h   =     [         1           2           5         ]             (   37   )                               
 
         [0203]    Thus, it can be seen that the algorithm effectively solves the system of equations after seven passes through the solution vector elements.  
         [0204]    The error in the values of h after each iteration is shown below:  
             [           Iteration        :           0       1       2       3       4       5       6       7       8       9       10             1      st                 element        :             -   1           -   1           -   1           -   1           -   1           -   1           -   1         0       0       0       0             2      nd                 element        :             -   2         6       6       2       2       0       0       0       0       0       0             3      rd                 element        :             -   5           -   5           -   5           -   1           -   1           -   1           -   1         0       0       0       0         ]                           
 
         [0205]    These values are plotted in the graph of FIG. 3, where a first line  4  represents changes to the error of the first solution vector element, h(1), a second line  5  represents changes to the error of the second solution vector element h(2), and a third line  6  represents changes to the error of the third solution vector element h(3). It can be seen that errors diminish as the algorithm proceeds.  
         [0206]    [0206]FIG. 4 is a graph showing the number of iterations carried out for each bit, i.e. the value of it when processing of each bit m has been completed. FIG. 5 shows the number of updates made to the solution vector h and auxiliary vector Q for each bit m.  
         [0207]    As described above, the auxiliary vector Q is updated as the algorithm progresses. The values of Q after each update of the vector Q and the solution vector h are set out below:  
             [         Update       0       1       2       3       4       5       6             1      st                 element           -   15         25       5         -   3           -   13         2       0             2      nd                 element           -   47         41         -   3         13         -   9           -   4         0             3      rd                 element           -   51           -   19           -   35         1         -   7           -   9         0         ]                           
 
         [0208]    As a further example, consider the system of equations set out below:  
         15 x +5 y− 2 z= 15  
         5 x +11 y +4 z =−37  
         −2 x+ 4 y +9 z =−55  (38)  
         [0209]    In this case the accurate solution of the equations is:  
             h   =     [         1             -   2               -   5           ]             (   39   )                               
 
         [0210]    It should also be noted that the value H is now equal to 256. Other parameters of the algorithm remain unchanged.  
         [0211]    The value of the solution vector after each pass of the algorithm is set out below. Again, it can be seen that the algorithm correctly solves the system of equations.  
                   [           Iteration        :           0       1       2       3       4       5       6       7       8       9       10       11       12             1      st                 element        :           0       0       0       0       0       0       0       0       0       0       0       1       1             2      nd                 element        :           0       0       0       0       0       0       0       0       0         -   2           -   2           -   2           -   2               3      rd                 element        :           0       0       0       0       0         -   8           -   8           -   8           -   6           -   6           -   6           -   5           -   5           ]           (   40   )                                 
 
         [0212]    [0212]FIG. 6 shows a variant to the algorithm described above with reference to FIG. 1. Many steps of FIG. 6 are identical to steps of FIG. 1. Such steps are identified by like reference numerals in both FIG. 1 and FIG. 6. Only steps which differ from FIG. 1 are described in further detail below.  
         [0213]    In FIG. 6, Step S 4  of FIG. 1 is replaced by step S 21 . It can be seen that in addition to setting d in accordance with equation (6):  
           d= 2 −m   H   (6)  
         [0214]    a two element array delta is established as follows:  
         delta[1]=d  (41)  
         delta[2]=−d  (42)  
         [0215]    In FIG. 1 step S 8  determines the value of arg (i.e. 1, 2 or 3) and chooses a different action in dependence upon the value. In FIG. 6, step S 8  is replaced with a single comparison at step S 22 :  
         arg&lt;3  (43)  
         [0216]    If the condition returns false, it can be determined that arg=3 (given that arg can only ever take values of ‘1’, ‘2’, and ‘3’). Therefore in accordance with the algorithm of FIG. 1, FIG. 6 shows that if the condition is false, no updates to h or Q are made, and control passes to step S 14  as in FIG. 1.  
         [0217]    If the condition of step S 22  of FIG. 6 is satisfied, it can be deduced that arg=1 or arg=2. In FIG. 1, if arg=1, steps S 9  and S 10  update h and Q using expressions including d. Similarly, if arg=2, steps S 11  and S 12  update h and Q using expressions including −d.  
         [0218]    In the variant of the algorithm shown in FIG. 6, the steps S 9  and S 11  of FIG. 1 are replaced with a single step S 23 , and similarly steps S 10  and S 12  are replaced with a single step S 24 . Both of steps S 23  and S 24  involve the array delta and more particularly the element arg of the array delta.  
         [0219]    Step S 23  updates h according to equation (44) and step S 24  updates Q according to equation (45).  
           h ( p )= h ( p )+ delta ( arg )  (44)  
           Q ( r )= Q ( r )+ delta ( arg )· R ( p,r ),∀ r:r= 1 , . . . ,N   (45)  
         [0220]    Given that the first element of delta contains d and the second element contains −d, it can be seen that equations (44) and (45) correctly correspond to equivalent operations of the algorithm of FIG. 1.  
         [0221]    [0221]FIG. 7 shows MATLAB® source code for implementing the algorithm illustrated in FIG. 6.  
         [0222]    [0222]FIG. 8 is a flow chart showing a further variant to the algorithm described above. The algorithm of FIG. 8 is intended to be used where upper and lower bounds are to be enforced upon elements of the solution vector h. For example, it may be known that the solutions lie between +5 and −5, and it may therefore speed up execution of the algorithm if only values in this range are considered.  
         [0223]    [0223]FIG. 8 illustrates an algorithm where all values h(p) of the solution vector h are constrained with the bounds of equation (46):  
           h   1   &lt;h ( p )&lt; h   2   (46)  
         [0224]    Before execution of the algorithm constants h 1  and h 2  are established. The establishment of these constants is not illustrated in the flowchart of FIG. 8. Steps of FIG. 8 which are identical to steps of FIG. 6 are identified by like reference numerals in both figures. Only steps which are different between the two Figures are described in further detail below.  
         [0225]    Having updated h(p) at step S 23 , a test is executed at a decision block of step S 25  to ensure that newly set value of h(p) lies between the bounds specified in equation (46). If the new value of h(p) satisfies equation (46), the algorithm proceeds to update the auxiliary vector Q at S 24  and to update Flag at step S 13  as described above with reference to FIG. 6.  
         [0226]    If step S 25  determines that the newly set value of h(p) does not satisfy equation (46), h(p) is again updated at step S 26 . This updating reverses the update of step S 23 , such that h(p) is equal to the value before execution of step S 23 . Therefore, h(p) has not been updated (because the attempted update did not satisfy equation (46)), Flag is not set to ‘1’, and execution of the algorithm continues at step S 14 .  
         [0227]    In yet a further embodiment of the present invention, the constants h 1  and h 2  of the embodiment of FIG. 8 are replaced by vectors h 1  and h 2  both of which have a length equal to the solution vector h. The vector h 1  contains a lower bound for each element of the solution vector, and the vector h 2  contains an upper bound for each element of the solution vector. Thus, the condition of step S 25  becomes:  
           h   1 ( p )&lt; h ( p )&lt; h   2 ( p )  (46)  
         [0228]    The embodiments of the present invention described above are concerned with the application of the algorithm to systems of equations which have real valued solutions. The present invention is also applicable to the solution of systems of equations having complex valued solutions, and the application of the invention to such systems of equations is described below.  
         [0229]    Consider the system of equations set out below:  
         ( a   1   +a   2   i ) x +( b   1   +b   2   i ) y +( c   1   +c   2   i ) z=A   1   +A   2   i    
         ( d   1   +d   2   i ) x +( e   1   +e   2   i ) y +( f   1   +f   2   i ) z=B   1   +B   2   i    
         ( g   1   +g   2   i ) x +( h   1   +h   2   i ) y +( j   1   +j   2   i ) z=C   1   +C   2   i   (47)  
         [0230]    where each of the unknown variables x, y and z is a complex number defined as follows:  
         
       x=x 
       1 
       +x 
       2 
       i  
     
         
       y=Y 
       1 
       +y 
       2 
       i  
     
           z=z   1   +z   2   i   (48)  
         [0231]    and  
         i={square root}{square root over (− 1 )} 
         [0232]    From equation (1) above, the system of equations (47) can be expressed as follows:  
             R   =     [             a   1     +       a   2        i               b   1     +       b   2        i               c   1     +       c   2        i                   d   1     +       d   2        i               e   1     +       e   2        i               f   1     +       f   2        i                   g   1     +       g   2        i               h   1     +       h   2        i               j   1     +       j   2        i             ]             (   49   )               h   =     [             x   1     +       x   2        i                   y   1     +       y   2        i                   z   1     +       z   2        i             ]             (   50   )               β   =     [             A   1     +       A   2        i                   B   1     +       B   2        i                   C   1     +       C   2        i             ]             (   51   )                               
 
         [0233]    To solve the system of equations, a matrix A, and vectors b, and c are created from the data set out above. A is a 2N by 2N real-valued coefficient matrix, b is a real-valued solution vector of length 2N, and c is a real-valued right hand side vector of length 2N, where N is the number of unknown variables (i.e. N=3 in this example).  
         A, b, and c are set as described below:  
           A (2 m− 1,2 n− 1)= Re{R ( m,n )} 
           A (2 m ,2 n )= Re{R ( m,n )} 
           A (2 m −1,2 n )= Im{R ( m,n )} 
           A (2 m ,2 n− 1)=− Im{R ( m,n )}  (52)  
           b (2 n− 1)= Re{h ( n )} 
           b (2 n )= Im{h ( n )}  (53)  
           c (2 n− 1)= Re{β ( n )} 
           c (2 n )= Im{β ( n )}  (54)  
         [0234]    Where Re{ } is a function returning the real coefficient of a complex number, and Im{ } is a function returning the imaginary coefficient of a complex number.  
         [0235]    Thus, in the case of equations (47), A, is set as follows:  
             A   =     [           a   1           a   2           b   1           b   2           c   1           c   2               -     a   2             a   1           -     b   2             b   1           -     c   2             c   1               d   1           d   2           e   1           e   2           f   1           f   2               -     d   2             d   1           -     e   2             e   1           -     f   2             f   1               g   1           g   2           h   1           h   2           j   1           j   2               -     g   2             g   1           -     h   2             h   1           -     j   2             j   1           ]             (   55   )                               
 
         [0236]    The present invention is often used to solve normal equations. Where normal equations are solved, the matrix A is such that:  
         b 1 =d 1  a 1 &gt;0  
         b 2 =−d 2  e 1 &gt;0  
         f 1 =h 1  j 1 &gt;0  
         f 2 =−h 2  a 2 =0  
         c 1 =g 1  e 2 =0  
         c 2 =−g 2  j 2 =0  (55a)  
         [0237]    The vectors b and c are set as follows:  
             b   =     [           x   1               x   2               y   1               y   2               z   1               z   2           ]             (   56   )               c   =     [           A   1               A   2               B   1               B   2               C   1               C   2           ]             (   57   )                               
 
         [0238]    Having created the vectors a and b and the matrix A set out above, the methods for solving real valued equations set out above with reference to FIGS. 1, 6 and  8  can be used where:  
         R=A  
         h=b  
         β=c  (58)  
         [0239]    Values of the solution vector h set by the algorithm can then be used to determine both the real and imaginary components of the complex numbers x, y and z, and to create the complex valued solutions.  
         [0240]    In alternative embodiments of the present invention, the algorithms described above are modified such that the algorithm operates directly on R, h and β as set out in equations (49) to (51) above. These modifications are now described with reference to the flow chart of FIG. 9.  
         [0241]    The algorithm used to solve equations involving complex numbers is based upon that of FIG. 1, but steps S 7  to S 12  are replaced by the steps shown in FIG. 9. Step S 27  of FIG. 9 replaces step S 7  of FIG. 1, step S 28  of FIG. 9 replaces step S 8  of FIG. 1, and steps S 29  to S 36  replace steps S 9  to S 12  of FIG. 1. Steps S 13  and S 14  are identical to the eqivilent steps of FIG. 1 and are shown in dotted lines. They are included in FIG. 9 for the sake of clarity.  
         [0242]    Referring to FIG. 9, at step S 27 , the following expression is evaluated:  
               arg   =     argmin        {       Re        {     Q        (   p   )       }       ,       -   Re          {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       R        (     p   ,   p     )       ·     d   2         }              
        where           (   59   )               argmin   =     {           1   ,           if        {       min        (       Re        {     Q        (   p   )       }       ,       -   Re          {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       -   Im          {     Q        (   p   )       }       ,       -     R        (     p   ,   p     )         ·     d   2         )       =     Re        {     Q        (   p   )       }         }                 2   ,           if                   {       min        (       Re        {     Q        (   p   )       }       ,       -   Re          {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       -   Im          {     Q        (   p   )       }       ,       -     R        (     p   ,   p     )         ·     d   2         )       =       -   Re          {     Q        (   p   )       }         }                 3   ,           if                   {       min        (       Re        {     Q        (   p   )       }       ,       -   Re          {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       -   Im          {     Q        (   p   )       }       ,       -     R        (     p   ,   p     )         ·     d   2         )       =     Im        {     Q        (   p   )       }         }                 4   ,           if                   {       min        (       Re        {     Q        (   p   )       }       ,       -   Re          {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       -   Im          {     Q        (   p   )       }       ,       -     R        (     p   ,   p     )         ·     d   2         )       =       -   Im          {     Q        (   p   )       }         }                 5   ,           if                   {       min        (       Re        {     Q        (   p   )       }       ,       -   Re          {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       -   Im          {     Q        (   p   )       }       ,       -     R        (     p   ,   p     )         ·     d   2         )       =       -   R          {     Q        (   p   )       }         }                       (   60   )                               
 
         [0243]    and Re{ } and Im{ } are as defined above.  
         [0244]    The value of arg set at step S 27  is checked at the decision block of S 28 . This determines the update (if any) to be made to the elements of the solution vector h and the auxiliary vector Q.  
         [0245]    If arg is 1, then the element p of the solution vector h currently under consideration is updated as follows at step S 29 :  
           h ( p )= h ( p )+ d   (61)  
         [0246]    Equation (61) means that the real part of h(p) is increased by d.  
         [0247]    All elements of the auxiliary vector Q are also updated, in accordance with equation (62) at step S 30 .  
           Q ( r )= Q ( r )+ dR ( p,r ),∀ r: 1≦ r≦N   (62)  
         [0248]    If arg is 2, then the element p of the solution vector h currently under consideration is updated as follows at step S 31 :  
           h ( p )= h ( p )− d   (63)  
         [0249]    Equation (64) means that the real part of h(p) is decreased by d.  
         [0250]    All elements of the auxiliary vector Q are also updated, in accordance with equation (62) at step S 32 .  
           Q ( r )= Q ( r )− dR ( p,r ),∀ r: 1≦ r≦N   (64)  
         [0251]    If arg is 3, then the element p of the solution vector h currently under consideration is updated as follows at step S 33 :  
           h ( p )= h ( p )+ i·d   (65)  
         [0252]    Equation (65) means that the imaginary part of h(p) is increased by d.  
         [0253]    All elements of the auxiliary vector Q are also updated, in accordance with equation (66) at step S 34 .  
           Q ( r )= Q ( r )+ d·i·R ( p,r ),∀ r: 1≦ r≦N   (66)  
         [0254]    If arg is 4, then the element p of the solution vector h currently under consideration is updated as follows at step S 35 :  
           h ( p )= h ( p )− i·d   (65a)  
         [0255]    Equation (65a) means that the imaginary part of h(p) is decreased by d.  
         [0256]    All elements of the auxiliary vector Q are also updated at step S 36 , in accordance with equation (66a).  
           Q ( r )= Q ( r )− d·i·R ( p,r ),∀ r: 1≦ r≦N   (66a)  
         [0257]    If arg=1, 2, 3, or 4, Flag is set to ‘1’ at step S 13 , control then passes to step S 14  and the algorithm proceeds as described above with reference to FIG. 1.  
         [0258]    If arg=5, no updates are made to h or Q and control passes to step S 14 , and the algorithm proceeds as described above.  
         [0259]    Step S 27  described above, requires the minimum of five values to be identified so as to determine the next course of action:  
       min                   {       Re        {     Q        (   p   )       }       ,       -   Re                     {     Q        (   p   )       }       ,     Im        {     Q        (   p   )       }       ,       -   Im          {     Q        (   p   )       }       ,       -     R        (     p   ,   p     )         ·     d   2         }                           
 
         [0260]    a method for finding this minimum, and efficiently identifying the correct action (at step S 28  of FIG. 9), is now described with reference to the flowchart of FIG. 10  
         [0261]    The algorithm takes as input two values a and b. a is input at a step S 101  and is set to be equal to the real part of Q(p), and b is input at a step S 102  and is set to be equal to the imaginary part of Q(p).  
         [0262]    A decision block S 103  determines whether or not a is greater than 0. If a is positive, a is set to be equal to the negative of its current value at a step S 104 , and a variable Ja is set to ‘1’ at step S 105 . If a is not positive, the variable Ja is set to ‘0’ at step S 106 .  
         [0263]    b is processed in a similar manner. Step S 107  checks whether b is positive. If b is positive, it is set to the negative of its current value at step S 108 , and a variable Jb is set to ‘1’ at step S 109 . If b is not positive, Jb is set to ‘0’ at step S 110 .  
         [0264]    Thus after execution of steps S 101  to S 110 , Ja is set to ‘1’ if the input value of a was positive, and ‘0’ if the input value of a was not positive. Jb is similarly set. Given the action of steps S 104  and S 108 , it can be seen that both a and b will now not be positive.  
         [0265]    At step S 111  a check is made to determine whether a is greater than b. If this condition is true a variable J 1  is set to be equal to Jb, a variable J 2  is set to ‘0’ and a variable c is set to be equal to b. This is accomplished by step S 112 . If the condition of step S 111  is false, a variable J 1  is set to be equal to Ja, a variable J 2  is set to ‘0’ and c is set to be equal to a. This is accomplished by step S 113 .  
         [0266]    At step S 114 , a check is made to determined whether or not c is greater than −R(p,p)d/2. If this condition returns true, a variable J 3  is set to be equal to ‘1’ at step S 115 . If this condition is false, the variable J 3  is set to be equal to ‘0’ at the step S 116 .  
         [0267]    The method of FIG. 10 is such that after execution of all steps, the variables J 1 , J 2  and J 3  are set as follows:  
         [0268]    J 3 =1         no update necessary  
         [0269]    J 3 =0         update necessary  
         [0270]    J 2 =1         update to imaginary part necessary  
         [0271]    J 2 =0         update to real part necessary  
         [0272]    J 1 =1         update should comprise subtraction  
         [0273]    J 1 =0         update should comprise addition  
         [0274]    It will be appreciated that the modifications made to the algorithm of FIG. 1, as illustrated in FIGS. 6 and 8 can similarly made to the algorithm of FIG. 9. For example, the algorithm of FIG. 9 can be implemented using an array delta, similar to that used in FIG. 6, where:  
         delta[1]=d;  
         delta[2]=−d,  
           delta[ 3 ]=d·i;    
           delta[ 4 ]=−d·i;    
         [0275]    Having established delta in this way, the algorithm can proceed by simply adding the appropriate element of delta to the appropriate element or elements of h or Q as described above. FIG. 11 shows this variant of the algorithm of FIG. 9 where steps S 28  to S 36  are replaced by steps S 37  to S 39 . It will be appreciated that additionally the array delta must be initialised and updated after each update of d in accordance with equation (67). This is not shown in FIG. 11. FIG. 12 shows a MATLAB program implementing the algorithm illustrated in FIG. 11.  
         [0276]    The description presented above has illustrated various implementations of algorithms in accordance with the invention. It will be appreciated that various amendments can be made to these algorithms without departing from the invention.  
         [0277]    For example, in the description set out above execution ends when the number of iterations it for any particular bit m reaches a predetermined limit Nit. It will be appreciated that execution need not end in this circumstance. Instead, a timer t may be set to ‘0’ each time m is updated, and execution can end if this timer exceeds a predetermined time threshold.  
         [0278]    The vectors h and Q need not necessarily be initialised as indicated above. Indeed, the initial value for h should usually be substantially centrally positioned in the range −H to H in which solutions are being sought, so as to obtain quick convergence. In some embodiments of the present invention the auxiliary vector need not be created. Instead the vector β is used directly, and is updated in a manner similar to the manner described above for vector Q, although it will be appreciated that different updates will be required. Suitable updates for such embodiments of the invention are set out in the derivation of the algorithm presented below.  
         [0279]    It has been described above that d is updated by dividing the previous value of d by two. This is preferred because considerable benefits are achieved because computations involving multiplication of d can be carried out using efficient bit shift operations in place of relatively inefficient multiplication operations. However, it will be appreciated that alternative methods for updating d may be used in some embodiments of the invention. For example, if d is updated by division by a power of two such as four or eight, computations can still be efficiently implemented by carrying out two or three bit shifts instead of the single bit shifts required when d is updated by division by two.  
         [0280]    In preferred embodiments of the present invention, each value of the solution vector h is represented by a fixed point binary word. This is particularly beneficial given that mathematical operations can typically be carried out more efficiently using fixed point arithmetic. Furthermore, a fixed point representation is likely to be acceptable because the different unknown variables are likely to have an approximately equal magnitude.  
         [0281]    In circumstances where a fixed point representation is inappropriate for the solution vector values, a conventional floating point representation can be used. Although the algorithm will operate more slowly with floating point values than with fixed point values, the algorithm still offers very favourable performance as compared with other methods for solving linear equations.  
         [0282]    It will be appreciated that the algorithm described above can be implemented in software or hardware. A software implementation will typically comprise appropriate source code executing on an appropriate microprocessor. For example, as shown in FIGS. 7 and 12, code can be created in MATLAB which is compiled to create object code which is specified in the instruction set of a microprocessor. Although MATLAB provides a convenient implementation for the algorithms of the invention, it will be appreciated that the algorithms could instead be coded in any one of the large number of widely available computer programming languages.  
         [0283]    A hardware implementation of the algorithm can either be provided by configuration of appropriate reconfigurable hardware, such as an application specific integrated circuit (ASIC), field programmable gate array (FPGA), or a configurable computer (CC), or alternatively by a bespoke microprocessor built to implement the algorithm. FIGS.  13  to  18  illustrate a device and microprocessor configured to solve linear equations.  
         [0284]    [0284]FIG. 13 schematically illustrates the general architecture of a device configured to solve linear equations. The device comprises a host processor  7  which is a general purpose microprocessor of known form configured to control operation of the device by running appropriate program code. This program code can be stored on a non volatile storage device  8  (e.g ROM or FLASH memory) and read by the general purpose microprocessor  7  when it is to be executed. Alternatively the program code can be copied to a volatile memory  9  (e.g RAM) prior to execution. In order to solve linear equations using the method of the invention, the device comprises an equation solving microprocessor  10 , operation of which is controlled by the host processor  7 . The aforementioned components of the device are connected together by a host bus  11 .  
         [0285]    [0285]FIG. 14 illustrates the architecture of the equation solving microprocessor  10  in further detail. The equation solving microprocessor  10  comprises a controller  12 , the function of which is to control operation of the equation solving microprocessor  10 . The equation solving microprocessor  10  further comprises a h-block  13 , which stores and updates the solution vector h of the algorithm, a Q-block  14  which stores and updates the auxiliary vector Q of the algorithm, and a R-block  15  which stores the coeffient matrix R used by the algorithm. The equation solving microprocessor  10  also comprises a minimisation block  16  which finds the minimum of a number of values. Components of the equation solving microprocessor  10  are connected together by an internal bus  17 , along which control instructions and data can pass.  
         [0286]    [0286]FIG. 15 shows the structure of R-block  15  in further detail. The R-block  15  comprises a storage element  18  which stores the elements of the coefficient matrix R. The R-block also comprises multiplexers  19 ,  20  and  21  which generate a column address signal  22 , and a row address signal  23  for reading data from and writing data to the storage element  18  in the manner described below. Finally, the R-block comprises a bit-shift element  24  which can left-shift a value presented at its input by a value determined by the current value of d.  
         [0287]    The multiplexers  19 ,  20 ,  21  generate address signals as follows. The row multiplexer and the column multiplexer  21  each have selection inputs connected to a common input line  25  which can carry an initialisation signal Init which is typically received from the controller  12  (FIG. 14). If an initialisation signal is received, this indicates that data is to be written to the storage element  18  representing the coefficients of the equations which are to be solved. In this case the row multiplexer should generate a row address  23  which is equal to the input row address  26 , and similarly the column multiplexer  21  should generate a column address  22  which is equal to the input column address  27 . During initialisation, the row address  26  and column address  27  will typically count through elements of the matrix R writing data provided to the storage element  18  on input line  28  to the appropriate element of the storage element  18 .  
         [0288]    When initialisation has been completed, the Init signal will no longer be received by the selection input of the multiplexers  20 ,  21 , and therefore the row address is determined by the input  29  to the row multiplexer  20  and the column address is determined by the input  30  to the column multiplexer  21 . The input  30  is the value p of the algorithm as described above. The input  29  is connected to the update multiplexer  19  whose output is dependent upon whether values of R are being read for purposes of analysis (e.g step S 7  of FIG. 1) or update of Q (e.g. steps S 10  and S 12  of FIG. 1). The update multiplexer  19  has a selection input  31  which indicates either update or analysis. This can conveniently be achieved, for example, by a single bit input where ‘0’ indicates update and ‘1’ indicates analysis.  
         [0289]    If R is being read for purposes of analysis (i.e. the selection input  31  is set to ‘1’), then the element R(p,p) is required (see step S 7  of FIG. 1). Therefore, update multiplexer output  29  is set to be equal to the input  32  which is equal to p, and the row multiplexer  20  subsequently generates a row address  23  equal to p.  
         [0290]    If R is being read for purposes of update of the auxiliary vector Q (i.e. selection input is set to ‘0’) R(p,r) is required (see steps S 10  and S 12  of FIG. 1). Therefore the update multiplexer output  29  is equal to the input  33  which is set to be equal to r, and the row multiplexer  20  subsequently generates a row address  23  equal to r.  
         [0291]    In either case, the column multiplexer  21  generates a column address  22  which is equal to p.  
         [0292]    From the preceding description, it can be seen that the multiplexers  19 ,  20  and  21 , act together to ensure that the correct row address  23  and column address  22  are sent to the storage element  18 . Having read the appropriate element of R from the storage element  18 , the output needs to be multiplied by d regardless of whether it is being used for update or analysis (see steps S 7 , S 10  and S 12  of FIG. 1). As described above, this can be achieved by a bit shift (assuming that d is a power of 2), the number of bits to shift being determined by the expression:  
         log 2 d  (68)  
         [0293]    which is passed to an input  34  of the bit shift block  24 . The output  35  from the bit shift block  24  is then correctly set as follows:  
         d·R(p,r), if update  (68a)  
         d·R(p,p), if analysis  (68b)  
         [0294]    [0294]FIG. 16 illustrates the h-block  13  of FIG. 14 in further detail. The h-block comprises a storage element  36  for storing elements of the solution vector h, an update/reading multiplexer  37 , and an adder-subtractor  38 .  
         [0295]    To initialise elements of the solution vector, an initialisation signal  39  is input to the storage element  36 . Upon receipt of this signal, all elements of the solution vector h are initialised to ‘0’.  
         [0296]    The update/reading multiplexer  37  sends an appropriate address  39   a  to the storage element  36 . The address  39   a  sent to the storage element  36  is determined by the signal provided to a selection input  40  of the multiplexer  37 . If the selection input  40  is set to update, the address p provided at input  41  becomes the address  39   a.    
         [0297]    If the selection input  40  is set to read (i.e. the values of h are being read to obtain the solution of the equations), a read address  42  input to the multiplxer  37  becomes the address  39 . In the case of a read operation, the input read address will count through all elements of h in turn, and each element will be transmitted to the internal bus  17  of the equation solving microprocessor  10  in turn.  
         [0298]    In the case of an update operation, a single value of p is provided to the multiplexer  37 , and an input signal I 1  is provided to the adder/subtractor  38 . The signal I 1  indicates whether the element of the solution vector h(p) is to be updated by adding d or subtracting d. Upon receipt of the address  39 , the storage element  36  outputs the appropriate element to the adder/subtractor  38  along a line  44 . By analysing the signal I 1  the adder/subtractor  38  can determine whether its output should be set to h(p)+d or h(p)−d. The appropriate expression is calculated, and written back to the storage element  36  along a line  45 .  
         [0299]    [0299]FIG. 17 shows how the h-block  16  can be implemented when the algorithm is used to solve equations having complex values (for example using the algorithm of FIGS. 10, 11 or  12 ). It can be seen that in addition to the components illustrated in FIG. 16, the update/read multiplexer  37  has an additional input  46  which carries a signal I 2 . This signal is sent to the storage element  36  along a line  47  along with the address p when the selection input is set to update. The signal I 2  indicates whether the update should be made to the real part or the imaginary part of h(p).  
         [0300]    [0300]FIG. 18 illustrates the Q-block  14  of FIG. 14 in further detail. It can be seen that the Q-block comprises a storage element  48  for storing the vector Q, an address multiplexer  49 , a data multiplexer  50  and an adder/subtractor  51 .  
         [0301]    The address multiplexer  49  has a selection input  52  which selects update, analysis, or initialisation mode. The address multiplexer  49  also has three data inputs, a first data input  53  carries the value r used in the algorithm, a second data input  54  carries the value p used in the algorithm, and a third data input  55  carries initialisation address data.  
         [0302]    When the selection input  52  of the address multiplexer  49  is set to initialise, an output  56  of the address multiplexer carries the initialisation address supplied at the third input  55  of the address multiplexer. When the selection input  52  is set to update, the output  56  carries the value r supplied at the first input  53  of the address multiplexer  49 . When the selection input  52  is set to analysis, the output  56  of the address multiplexer  49  is set to the value p supplied at the second input  54  to the address multiplexer  49 . The output  56  of the address multiplexer is passed to the storage element  48  as an address for reading or writing data.  
         [0303]    Data is written to the storage element  48  only in the update and initialisation modes. Therefore, the data multiplexer  50 , has a selection input  57  having two recognised values, update and initialise (it will be appreciated that in practice this selection input may be common with the selection input  52  of the address multiplexer  49 , although the behaviour of the data multiplexer is not well defined when the selection input is set to analysis).  
         [0304]    The data multiplexer  50  comprises a first data input  58  carrying initialisation data (elements of the vector β) and a second data input  59  carrying update data. When the selection input  57  of the data multiplexer  50  is set to initialise, data from the first data input  58  is provided at the output  60  of the data multiplexer  50 . When the selection input  57  is set to update, data from the second data input  59  is provided at the output  60 .  
         [0305]    In operation, when the selection inputs of both multiplexers  49 ,  50  are set to initialise, a sequence of addresses (counting through all elements of the vector Q) are provided at the third data input  55  of the address multiplexer  49 . In synchronisation with these addresses, the appropriate data is provided at the first data input  58  of the data multiplexer  50 . The data and addresses are provided to the storage element  48 , and the storage element is then set to contain the appropriate elements of the vector Q.  
         [0306]    When the selection inputs of both multiplexers  49 ,  50  are set to update, the address r is provided to the storage element  48  by the multiplexer  49 . The appropriate entry from the storage element  48 , Q(r) is provided at an output  61  of the storage element, and passed to a first input  62  of the adder/subtractor  51 . A second input  63  of the adder/subtractor is set to d·R(p,r) (provided by the R-block  15  of FIG. 14). The adder/subtractor also includes an input  64  carrying a signal I 1  which indicates whether the adder/subractor  51  should perform an addition or subtraction operation. On receipt of appropriate data at its inputs an output  65  of the adder/subtractor is set to the output of the operation preformed, and this data is fed to the data multiplexer  50 . Given that the selection input  57  is set to update, the data provided at the input  59  is provided to the output  60  of the data multiplexer, and from there to the storage element  48  where it is written to the element of Q specified by the address  56 .  
         [0307]    Where the algorithm is being implemented to solve complex valued equations, the adder/subtractor comprises a further input  66  carring an signal I 2  which indicates whether the update should be applied to the real part or the imaginary part of the appropriate element of the vector Q.  
         [0308]    In analysis mode, the Q-block is required to provide a current value of Q(p), with no update functions being needed. The data multiplexer  50  is therefore not involved with the analysis operation. The address multiplexer  49  provides the address p provided at its second input  54  to the storage element  48 . The appropriate value of Q(p) is read from the storage element  48  and provided at the output  61  of the storage element  48 .  
         [0309]    [0309]FIG. 19 provides an implementation for the minimisation block  16  of FIG. 14. This implementation corresponds to that illustrated in FIG. 10, and is therefore configured for use in algorithms solving equations having complex valued solutions. The minimisation block comprises first and second converter blocks  67 ,  68 , first and second comparators  69 ,  70  and first and second multiplexers  71 ,  72 .  
         [0310]    Input data is provided to the first and second converter blocks  67 ,  68 . The first converter block  67  takes as input the real part of a value Q(p), and the second converter block  68  takes as input the imaginary part of the value Q(p). The converter blocks  67  and  68  provide two outputs, a first output  67   a ,  68   a  indicates the sign of the input data, and second output  67   b ,  68   b  indicates the modulus of the input data. The modulus outputs  67   b ,  68   b  are input to the first comparator  69 , which generates a signal I 2  which at an output  73 . The signal I 2  is defined as being ‘0’ if the value provided at output  67   b  is greater than that provided at  68   b , and ‘1’ if the value provided at output  68   b  is greater than or equal to that provided at output  67   b .  
         [0311]    The values provided at outpus  67   b  and  68   b  are also provided to the first multiplexer  71 , along with the output value  73  created by the first comparator  69 . The output  73  of the first comparator acts as a selection input to the first multiplexer  71 . The first multiplexer  71  then acts to provide the larger of the values  67   b ,  68   b  at its output  74 .  
         [0312]    The output  74  of the multiplexer  71  is input to the second comparator  70 , along with the value d·R(p,p) at an input  75 , which is computed by the R-block  16  of FIG. 15. The second comparator  70  then produces an output  76  which is a signal I 3  which is set such that 13 is set to ‘0’ if the input  74  greater than the input  75 , and ‘1’ if the imputer  75  is greater than or equal to the input  74 .  
         [0313]    The second multiplexer  72  is a bit multiplexer taking as input the signs of the input data generated by the converter blocks  67 ,  68 . In dependence upon the value of the signal I 2  output from the first comparator  69 , the multiplexer generates an output  77  which is a signal I 1 .  
         [0314]    The signals I 1 , I 2  and I 3  can together be used to determine what update (if any) should be made to the elements of Q and h. I 3  indicates whether or not the current iteration is successful. If I 3  is equal to ‘0’, the element h(p) is updated, and all elements of Q are updated. If I 3  is equal to ‘1’, no updates are necessary.  
         [0315]    I 2  indicates whether the real or the imaginary part of the appropriate elements should be updated. If I 2  is equal to ‘0’ the real part of the appropriate values is updated, while if I 2  is equal to ‘1’ the imaginary part of the appropriate values is updated. I 1  indicates whether the update should comprise addition or subtraction. If I 1  is set to ‘0’, addition is used, and if I 1  is set to ‘1’, subtraction is used.  
         [0316]    Having described the structure and function of the individual components of the equation solving microprocessor  10  (FIG. 14), operation of the processor is now described.  
         [0317]    To perform initialisation, the equation solving microprocessor receives a signal to cause initialisation, for example from the host processor (FIG. 13). The controller  12  of the equation solving microprocessor then generates an appropriate internal initialisation signal which is communicated to all blocks of the equation solving microprocessor, via the internal bus  17 . Upon receipt of this signal the h-block  13 , the Q-block  14  and the R-block  15  perform initialisation as described above. The data required for initialisation may be located in a predefined read only memory, or an address where the data is to be found may be provided to the controller  12  for onward transmission to the appropriate block. Parameters used by the algorithm (e.g. N, M and Nit) can either be provided to the controller, or alternatively specified within the microprocessor  10 .  
         [0318]    When initialisation is complete, the controller begins executing an algorithm in accordance with the invention using the blocks of the microprocessor  10  to update values of h and Q as appropriate. In some embodiments of the invention data that is required to be passed between blocks of the microprocessor is passed directly between blocks, under the control of the sending and/or receiving block. In alternative embodiments of the invention all data is passed between blocks as directed by the controller  12 . When the controller determines that the equations have been solved, such that the vector h contains the solution of the equations, the vector h can then be copied from the storage element  36  (FIGS. 16 and 17) to an appropriate location within the device, where the solutions can be used as required.  
         [0319]    It will be appreciated that although a specific hardware implementation of algorithms of the invention has been described above, numerous modifications could be made to the implementation while remaining within the scope and spirit of the invention. Such modifications will be readily apparent to those of ordinary skill in the art.  
         [0320]    The preceding description has described algorithms in accordance with the invention, and hardware suitable for implementing such algorithms. The mathematical basis of algorithms in accordance with the invention is now described.  
         [0321]    To aid understanding of the invention, the known co-ordinate descent optimisation method for minimising a function of many variables is presented. The co-ordinate descent optimisation method seeks to find:  
         min{J(h)}  (69)  
         [0322]    where J is a function of many variables; and  
         [0323]    h is an N-dimensional vector containing the variables.  
         [0324]    Thus we want to find the value of h when the function J is a minimum.  
         [0325]    In many circumstances the elements of h have a maximum amplitude H. Therefore, the minimisation problem is considered for a subset of values of h defined by equations (70) and (71).  
         hεU  (70)  
         [0326]    where:  
         U={ h ( m ):1 ≦m≦NΛ|h ( m )|&lt; H};   (71)  
         [0327]    where:  
         [0328]    h(m) are elements of the vector h and H is a known number such that H&gt;0.  
         [0329]    Define a vector ei where the ith coordinate is ‘1’, and all other coordinates are ‘0’. Let h 0  be an initial value of the solution vector h and let α 0  be an initial value of the step size parameter (that is d in the algorithms described above).  
         [0330]    h k  is defined to be the value of the solution vector h after some number of iterations k, where k is a positive integer. α k  is the step size parameter after k iterations.  
         [0331]    A vector p k  is also defined, according to the equation:  
         p k =e i     k     (72)  
         [0332]    where:  
           i   k   =k−N ( div ( K,N ))+1  (73)  
         [0333]    where div (A,B) is a function whose result is defined as the integer part of result of dividing A by B.  
         [0334]    It can be seen from equation (73) that i k  will repeatedly count from 1 to N as k is increased. This results in values of p cycling through e 1  to e N : p, e 1 ,p 2 =e 2 , . . . , p N =e N ,p N+1 =e 1 ,p N+2 =e 2 , . . . ,p 2N =e N   (74)  
         [0335]    The value of the function J is calculated at the point h=h k +α k p k  and two conditions are checked:  
         h k +α k p k εU  (75)  
           J ( h   k +α k   p   k )&lt; J ( h   k )  (76)  
         [0336]    If the conditions of equations (75) and (76) are satisfied, then the solution vector h is pdated as follows:  
           h   k+1   =h   k +α k   p   k .  (77)  
         [0337]    Given that that update involves a scalar multiple of the vector p k  and given also that p k  is a vector containing only a single non-zero element, it can be seen that equation (77) means that h k +1 will be identical to hk save for a single element.  
         [0338]    The step size parameter is not updated at this stage, as indicated by equation (78)  
         α k+1 =α k   (78)  
         [0339]    If the conditions of equations of (75) and (76) are not satisfied, the value of the function J is calculated at the point h=h k −α k p k  and two conditions are again checked:  
         h k −α k p k εU  (79)  
           J ( h   k −α k p k )&lt; J ( h   k )  (80)  
         [0340]    If the conditions of equations (79) and (80) are satisfied, then the solution vector h is updated as follows:  
           h   k+1   =h   k −α k p k .  (81)  
         [0341]    Again, given that that update involves a scalar multiple of the vector p k  and given also that p k  is a vector containing only a single non-zero element, it can be seen that equation (81) means that h k+1  will be identical to h k  save for a single element.  
         [0342]    The step size parameter is not updated at this stage, as indicated by equation (82)  
         α k+1 a k   (82)  
         [0343]    The kth iteration is considered to be successful if either the conditions of equations (75) and (76) or equations (79) and (80) are satisfied. If neither of these conditions are satisfied, the iteration is considered to be unsuccessful.  
         [0344]    α k  is updated according to equation (83):  
               α     k   +   1       =     {             λ                   α   k       ,       if                   i   k       =       N   ⋀     h   k       =     h     k   -   N   +   1                         α   k     ,   otherwise                     (   83   )                               
 
         [0345]    where λ is a parameter of the algorithm, and is such that 0&lt;λ&lt;1  
         [0346]    The condition of equation (83) means that if the last N iterations involve no successful updates (i.e. the value of h has not changed), the step size parameter is updated by multiplication by λ. If however there is at least one update during the previous N iterations, the step size parameter is not updated. It should be recalled that N is defined to be the number of elements in the solution vector h, and it can therefore be seen that h is updated only when all elements have been processed and no update has occurred.  
         [0347]    The method described above is generally known as co-ordinate descent optimisation.  
         [0348]    It is known that for some arbitrary function J, the method converges to find the value of h for which the function is a minimum, providing that the function J is convex and differentiable on U, providing that α 0 &gt;0, and 0&lt;λ&lt;1 and providing that h 0  εU. This result is shown, for example, in Vasiliev, F. P.: “numerical methods for solutions of optimisation problems”, Nauka, Moscow 1988 (published in Russian), the contents of which is herein incorporated by reference.  
         [0349]    It is often necessary to solve the linear least squares problem. The linear least squares roblem is concerned with the minimisation of the function J specified in equation (84):  
           J ( h )=| Zh−d|   2 =( Zh−d ) T ( Zh−d )  (84)  
         [0350]    with respect to an unknown vector h, where Z is a known M×N matrix, d is a known vector of length N, and denotes the transpose of a matrix.  
         [0351]    This is discussed in Sayed, A. H., and Kailath, K.: “Recursive least-squares adaptive filters”, The Digital Signal Processing Handbook, CRC Press, IEEE Press  1998 , pages 21.1 to 21.37, which discussion is herein incorporated by reference.  
         [0352]    It can be shown that the minimisation of the function of equation (84) is equivalent to minimisation of a quadratic function.  
         [0353]    Equation (84) can be rewritten as:  
           J ( h )=h T   Z   T   Zh−h   T   Z   T   d−d   T   Zh+d   T   d   (85)  
         [0354]    by multiplying out the bracketed expressions of equation (84).  
         [0355]    Given that the purpose of the method is to minimise J with respect to h it can be concluded that the term d T d of equation (85) will not effect the minimisation process, and therefore can be removed from the expression without affecting the minimum value. Therefore minimisation of equation (85) is equivalent to minimisation of equation (86):  
           J ( h )= h   T   Z   T   Zh−h   T   Z   T   d−d   T   Zh   (86)  
         [0356]    A matrix R is defined according to equation (87):  
         R=Z T Z  (87)  
         [0357]    A matrix β is defined according to equation (88):  
         β=Z T d  (88)  
         [0358]    Equation (86) can then be rewritten using R and β as shown in equation (89):  
           J ( h )=h T   Rh−h   T β−β T   h   (89)  
         [0359]    Simplifying this expression yields:  
           J ( h )= h   T   Rh− 2 h   T β  (90)  
         [0360]    The expression h T Rh can be rewritten in terms of the elements of h, and R as follows, using the definitions of matrix multiplication, and the matrix transpose operation:  
                 h   T        Rh     =       ∑     m   =   1     N            ∑     n   =   1     N            R        (     m   ,   n     )            h        (   m   )            h        (   n   )                     (   91   )                               
 
         [0361]    Similarly the expression h T β can be written in form of equation (92):  
                 h   T        β     =       ∑     n   =   1     N            β        (   n   )            h        (   n   )                   (   92   )                               
 
         [0362]    Substituting equations (91) and (92) into equation (90) yields:  
               J        (   h   )       =         ∑     m   =   1     N            ∑     n   =   1     N            R        (     m   ,   n     )            h        (   m   )            h        (   n   )             -     2          ∑     n   =   1     N            β        (   n   )            h        (   n   )                       (   93   )                               
 
         [0363]    It can be seen that equation (93) is a quadratic function of h. Thus it can be seen that solving the linear least squares problem is equivalent to minimisation of the function of equation (93). Furthermore, it is also known that solving a system of linear equations of the form of equation (1):  
         Rh=β  (1)  
         [0364]    is equivalent to minimisation of the function of equation (93), for any set of normal linear equations. This is explained in, for example, Moon, Todd K., and Stirling, Wynn C.: “Mathematical methods and algorithms for signal processing”, Prentice Hall, 2000, section 3.4. “Matrix representations of least-squares problems”, pages 138-139. This explanation is incorporated herein by reference.  
         [0365]    Given that many sets of linear equations occurring the electronics and physics are normal linear equations, minimisation of the function of equation (93) has wide applicability in solving linear equations.  
         [0366]    The explanation presented above has set out a method for minimising a function J using the co-ordinate descent optimisation method. The material presented above has also set out the relationship between the minimisation of equation (93) and a set of linear equations of the form of equation (1).  
         [0367]    The present inventors have surprisingly discovered that applying the known co- ordinate descent method to the minimisation of equation (93) provides a particularly efficient method for solving linear equations.  
         [0368]    Minimisation of the function of equation (94) is considered:  
               J        (   h   )       =         1   2            ∑     m   =   1     N            ∑     n   =   1     N            R        (     m   ,   n     )            h        (   m   )            h        (   n   )               -       ∑     n   =   1     N            β        (   n   )            h        (   n   )                     (   94   )                               
 
         [0369]    This minimisation process finds values for the elements of h which minimise the function J(h). The matrix R and the vector β are known. It is known that the function of equation (94) is convex and differentiable. This is shown in Vasiliev, F. P.: “Numerical methods for solutions of optimisation problems”, Nauka, Moscow 1988 (published in Russian), page 345, which explanation is incorporated herein by reference. Therefore, as explained above, the co-ordinate descent optimisation method can be used to find the minimum value of the function J.  
         [0370]    During operation of the co-ordinate descent optimisation method, the following expressions are computed:  
         Δ J ( h   k )= J ( h   k +α k   e   i     k   )− J ( h   k )  (95)  
         [0371]    and  
         Δ J ( h   k )= J ( h   k α k   e   i     k   )− J ( h   k )  (96)  
         [0372]    It can be seen that equation (95) relates to the condition of equation (76) set out above, while equation (96) relates to the condition of equation (80) set out above.  
         [0373]    The inequality of equation (97) is also considered:  
         Δ J ( h   k )&lt;0  (97)  
         [0374]    It can be recalled that the values of the vector e i     k    has elements defined by:  
                 [          i   k       ]     i     =     {             1                 if                 i                -                i   k                 0                 otherwise                     (   98   )                               
 
         [0375]    Also, it is known that the matrix R is symmetric, given that the system of equations is normal. Therefore, substituting equation (94) into equation (95) yields:  
               Δ                   J        (     h   k     )         =         α   k     2          [         -   2          β        (   i   )         +     2          ∑     m   =   1     N              h     (   k   )            (   m   )            R        (     m   ,   i     )             +       α   k          R        (     i   ,   i     )           ]               (   99   )                               
 
         [0376]    where h (k) (m) are elements of the vector h k .  
         [0377]    If a vector Q is defined as:  
               Q        (   i   )       =         -   2        β                   (   i   )       +     2          ∑     m   =   1     N              h     (   k   )            (   m   )            R        (     m   ,   i     )                       (   100   )                               
 
         [0378]    then equation (99) can be written as:  
               Δ                   J        (     h   k     )         =         α   k     2          [       Q        (   i   )       +       α   k          R        (     i   ,   i     )           ]               (   101   )                               
 
         [0379]    Given that α k &gt;0, from equation (101), equation (95) can be rewritten as:  
         α k   R ( i,i )&lt;− Q ( i )  (102)  
         [0380]    Similarly, equation (96) can be rewritten as:  
         α k   R ( i,i )&gt; Q ( i )  (103)  
         [0381]    An auxiliary vector Q k  is defined as the vector Q after the kth iteration. If the (k+1)th iteration is not successful, then elements of h and Q can be updated as follows:  
         h k+ 1=h k   (104)  
         [0382]    and  
         Q k+ 1=Q k   (105)  
         [0383]    the (k+1)th iteration is successful, then elements of h and Q are updated as ollows:  
           h   (k+1) ( i )= h   (k) ( i )±α k   (106)  
           Q   (k+1) ( n )= Q   (k) ( n )±2α k   R ( n,i ), n= 1 , . . . ,N   (107)  
         [0384]    The vector h can be initialised to a vector h 0  where  
           h   0 ( n )=0 , n= 1 , . . . ,N   (108)  
         [0385]    Then, from equation (100), elements of Q are initialised as follows:  
           Q   0 ( n )=−2β( n ),  n= 1 , . . . ,N   (109)  
         [0386]    That is, each element of Q is set to be the negative of the corresponding element of β multiplied by two.  
         [0387]    Thus, the solution vector h is initialised to contain all ‘0’ values, while the auxiliary vector Q is initialised to be the negative of the vector β.  
         [0388]    As described above, multiplication operations may be avoided by setting H in accordance with equation (110):  
           H= 2 P+M     b     (110)  
         [0389]    where M b  is a positive integer, and P is any integer. α 0  is initialised to be  
       H   2                         
 
         [0390]    and λ is initialised to be  
         1   2     .                         
 
         [0391]    The multiplications described above can then be replaced by bit shift operations.  
         [0392]    The algorithms described thus far have used an auxiliary vector Q which is intialised n accordance with equation (4):  
         Q=−β  (4)  
         [0393]    However, some embodiments of the invention use β itself as an auxiliary vector. In such embodiments of the invention, the auxiliary vector update rule of equation (107) above becomes:  
         β (k+1) ( n )=β (k) ( n )±α k   R ( n,i ), n= 1 , . . . , N   (111)  
         [0394]    Similarly, the inequalities of equations (102) and (103) become:  
         α k   R ( i,i )&lt;2β( i )  (112)  
         [0395]    and:  
         α k   R ( i,i )&gt;−2β( i )  (113)  
         [0396]    [0396]FIG. 19 a  shows MATLAB source code for an algorithm which uses β in place of Q. This algorithm is based upon that of FIG. 7, but the code has been amended as set out above.  
         [0397]    The step size parameter α k  is updated if N consecutive iterations are not successful. At every update of α k , α k  is decreased by a factor of two. The algorithm described above is therefore referred to as the dichotomous coordinate descent algorithm for the solution of linear equations.  
         [0398]    From the description set out above, it can be observed that the algorithms of the invention solve linear equations by minimisation of an appropriate quadratic function. It has been explained above that it is known that such minimisation can be mployed to solve normal linear equations. However, it should be noted that the resent invention is not limited simply to normal linear equations, but is instead pplicable to a wider class of linear equations.  
         [0399]    In the methods described above, the elements of the solution vector h are analysed in predetermined order (i.e. from element ‘1’ to element N). However, it will be appreciated that elements of the solution vector h can be analysed in any convenient manner. For example, the values of h can be sorted on the basis of some corresponding auxiliary value (e.g. a corresponding element of the vector Q), and elements of the solution vector h can then be processed in that order. For some applications, ordering elements of the vector h in this way will provide a more rapid convergence, although this must of course be balanced against the computational cost of sorting the elements of h.  
         [0400]    It has been explained above that the present invention can be usefully applied in any application in which it is necessary to solve linear equations. Two such applications are now described.  
         [0401]    In a multiuser Code Division Multiple Access (CDMA) communications system, a plurality of users transmit data using a common collection of frequencies. A narrow band data signal which a user is to transmit is multiplied by a relatively broad band spreading code. Data is then transmitted using this broad band of frequencies. Each user is allocated a unique spreading code.  
         [0402]    A receiver needs to be able to receive data transmitted by a plurality of users simultaneously, each user using his/her respective spreading code. The receiver therefore needs to implement functions which allow the spreading code to be removed from the received data to yield the originally transmitted data. Typically filters are used to extract the spreading code to obtain the transmitted data. It should be noted that the process is complicated by interfering signals from multiple users, and also from different propagation paths which may be followed by different signals.  
         [0403]    [0403]FIG. 20 illustrates a receiver  80  suitable for use in a CDMA communications system. The receiver comprises a receiver circuit  81  including an antenna (not shown), an analog to digital convertor  82 , a bank of filters  83   a ,  83   b  and  83   c , an equation solving circuit  84  and a decision circuit  85 .  
         [0404]    Spread spectrum signals are received by the receiver circuit  81 , and converted into digital data by the analog to digital converter  82 . Digital data is then passed to all of the filters  83   a ,  83   b ,  83   c . Each filter of the bank of filters  83   a ,  83   b ,  83   c  relates to a unique user, and has filter coefficients selected to match the spreading code of that user. It will therefore be appreciated that in practical embodiments a receiver will include more than three filters as is illustrated in FIG. 20.  
         [0405]    If a single signal is transmitted at any one time, and this signal travels between a sender and the receiver  80  by a direct path, the output of the filters alone should provide the data which the sender intended to transmit. However, in the more likely and more complicated situation where interference between signals occurs, the outputs of the filters alone will be inconclusive. However, it is known that solving an equation:  
         Rh=β  (114)  
         [0406]    where R is the cross correlation matrix of the spreading sequences of all users;  
         [0407]    β is a vector containing the filter output values; and  
         [0408]    h is a solution vector  
         [0409]    will allow the originally transmitted data to be obtained.  
         [0410]    In general, for a system involving N users, there will be N filters, and the vector β will therefore have length N, and the matrix R will have size N×N.  
         [0411]    The cross correlation matrix R can be defined as  
         R=S T S  (115)  
         [0412]    where the matrix S contains the spreading codes, and is defined as follows:  
             S   =     [             s   1          (   1   )               s   2          (   1   )           ⋯           s   k          (   1   )                   s   1          (   2   )               s   2          (   2   )           ⋯           s   k          (   2   )               ⋮       ⋮       ⋰       ⋮               s   1          (   m   )               s   2          (   m   )           ⋯           s   k          (   m   )             ]             (   116   )                               
 
         [0413]    where s j (i) denotes the ith element of the spreading code for user j.  
         [0414]    As has been described above linear equations of the form shown in equation (114) can be solved using an algorithm in accordance with the invention. Therefore, the invention provides a novel multi user receiver apparatus, in which the equation (114) is solved as described above, thereby achieving the considerable performance benefits provided by solving equations in accordance with the invention.  
         [0415]    The equation solver  84  of FIG. 20 can either be implemented by means of a computer program carrying out the method of the invention executing on an appropriate microprocessor, or alternatively by means of hardware, such as that described above with reference with FIGS.  14  to  19 .  
         [0416]    The equation solver provides a vector h as output, and this is input to the decision circuit  85 , which then determines the nature of the transmitted data.  
         [0417]    It will be appreciated that the cross correlation matrix described with reference to equations (115) and (116) is merely exemplarily. Cross correlation matricies can be created in a variety of different ways which will be known to one skilled in the art. Regardless of how the cross correlation matrix is formed, a system of equations (114) is created which can be solved using methods in accordance with the present invention.  
         [0418]    It will also be appreciated that in addition to the components illustrated in FIG. 20 and described above, a CDMA receiver may require other components to function properly. The selection and use of these components will be readily apparent to those of ordinary skill in the art.  
         [0419]    The algorithms of the invention can also be employed in adaptive filtering applications such as echo cancellation in a hands free communications system.  
         [0420]    A system of interest in illustrated in FIG. 21. A signal travels along input line  86  and is output through a loudspeaker  87 . A further signal such a human voice (not shown) is passed to an input of a microphone  88 . It is desirable that the signal at the output  89  of the microphone  88  contains only the human voice signal, and none of the signal output by the loudspeaker  87 . However, in practice, some of the signal output by the loudspeaker  87  will be received by the microphone  88  such that the microphone output  89  comprises a combination of the human voice signal and part of the loudspeaker output signal (referred to as “the echo signal”). It is desirable to remove the echo signal present in the microphone output  89 .  
         [0421]    As shown in FIG. 21, the echo cancellation apparatus comprises a filter  90 , which is configured to provide an estimate  91  of the echo signal. This estimate  91  is subtracted from the microphone output signal  89  by a subtractor  92 . Therefore, if the echo is accurately estimated, an output  93  of the subtractor will be equal to the human voice signal.  
         [0422]    The echo cancellation apparatus comprises a filter coefficient setting circuit  94 , which takes as inputs a signal  86  which is input to the loudspeaker and the signal  89  which is output from the microphone. The circuit  94  should generate coefficients to allow the filter  90  to accurately model the echo signal.  
         [0423]    [0423]FIG. 22 shows the filter coefficient setting circuit  94  in further detail. It can be seen that the circuit comprises an auto correlation element  95 , a cross correlation element  965  and an equation solver  97 . The auto correlation element  95  finds the auto correlation of the signal  86 . The cross correlation element  96  finds the cross correlation between the signal  86  and the signal  89 . It is known that when a auto correlation matrix R and a cross correlation vector β have been generated, the optimal filter coefficients h can be found by solving the system of equations:  
         Rh=β  (117)  
         [0424]    where the auto correlation matrix R is generated according to the equation:  
                 R        (     m   ,   n     )       =       ∑     t   =   1       T   -          n   -   m                     x        (   t   )            x        (     t   +          n   -   m            )             ,     
          ∀   m     ,       n        :                   1     &lt;   m   &lt;     N   ⋀   1     &lt;   n   &lt;   N             (   118   )                               
 
         [0425]    where t=1, . . . , T are discrete time moments;  
         [0426]    and the cross correlation vector β is generated according to the equation:  
         β        (   n   )       =       ∑     t   =   1       T   -   n   -   1              x        (   t   )            y        (     t   +   n   -   1     )                                 
 
         [0427]    where x is the loudspeaker input signal  86  and y is the microphone output signal  89 .  
         [0428]    An echo cancellation system operating in the manner described above is described in US5062102 (Taguchi).  
         [0429]    Having generated a system of equations of the form of equation (117), algorithms in accordance with the present invention can be used to solve linear equations to determine a solution vector h containing optimal filter coefficients. Therefore, referring back to FIG. 22, the equation solver  97  can either be a microprocessor executing code in accordance with one of the algorithms described above, or alternatively hardware suitably configured to implement a suitable algorithm.  
         [0430]    It will be appreciated that although this application of the algorithm has been described with reference to echo cancellation, it is widely applicable in all cases where an adaptive filter is required, and where solving a system of linear equations yields appropriate filter coefficients. A suitable example system in which the invention could be beneficially employed is described in WO 00/38319 (Heping).  
         [0431]    Applications of the invention to CDMA receivers, and echo cancellers have been described above. However, it will be appreciated that many other applications exist which can benefit by the improved efficiency with which linear equations can be solved in accordance with the invention. For example, the invention can be used in tomographic imaging systems, where a large system of linear equations is solved to generate an image.  
         [0432]    Although the present invention has been described above with reference to various preferred embodiments, it will be apparent to a skilled person that modifications lie within the scope and spirit of the present invention.