WO2024258418A1 - Mitigation of multipath-induced error effects in a navigation satellite signal receiver - Google Patents

Mitigation of multipath-induced error effects in a navigation satellite signal receiver Download PDF

Info

Publication number
WO2024258418A1
WO2024258418A1 PCT/US2023/034191 US2023034191W WO2024258418A1 WO 2024258418 A1 WO2024258418 A1 WO 2024258418A1 US 2023034191 W US2023034191 W US 2023034191W WO 2024258418 A1 WO2024258418 A1 WO 2024258418A1
Authority
WO
WIPO (PCT)
Prior art keywords
accordance
determination
samples
iterative process
estimate
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
PCT/US2023/034191
Other languages
French (fr)
Inventor
Gang Xie
Brian C. GOODRICH
Zhe LIU
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Deere and Co
Original Assignee
Deere and Co
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from US18/373,913 external-priority patent/US20240411029A1/en
Application filed by Deere and Co filed Critical Deere and Co
Priority to AU2023456885A priority Critical patent/AU2023456885A1/en
Priority to CN202380099371.7A priority patent/CN121488173A/en
Priority to EP23895569.4A priority patent/EP4695629A1/en
Publication of WO2024258418A1 publication Critical patent/WO2024258418A1/en
Anticipated expiration legal-status Critical
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/428—Determining position using multipath or indirect path propagation signals in position determination
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/22—Multipath-related issues

Definitions

  • GNSS global navigation satellite systems
  • GNSS Global Navigation Satellite Systems
  • GPS Global Positioning System
  • GLONASS Global Orbiting Navigation Satellite System
  • BDS Chinese BeiDou System
  • WAAS Wide Area Augmentation System
  • QZSS Japanese Quasi-Zenith Satellite System
  • PN code pseudorandom code
  • PRN code pseudorandom noise code
  • the PN code is, for example, a series of ones and zeroes that looks random, but in fact uniquely identifies the satellite.
  • the GNSS receivers acquire and track the GNSS satellite signals. From its code and carrier tracking loops, a GNSS receiver can then generate range or range rate measurements, such as pseudorange measurements, carrier phase measurements and doppler shift measurements. In some instances, using other data, such as ephemeris data and almanac data, that is encoded in the received signal, “on top of” the pseudorandom code, using various modulation schemes (e.g., BPSK, BOC), the receiver can determine the time at which the signal was sent by the satellite. Using signals received from four or more satellites, the receiver can determine its position (e.g., on Earth).
  • Multipath is one of the major error sources in GNSS.
  • a signal from a GNSS satellite to a GNSS receiver is reflected by surroundings (e.g., buildings, mountains, lakes, other objects, and/or the ground)
  • this reflected signal called a multipath
  • a GNSS receiver can receive a composite of a line-of-sight direct-path signal and zero, one or more multiple signals from each satellite.
  • multipath signals are also internally processed by GNSS receivers, thus introducing errors to code and carrier phase measurements and reducing GNSS positioning accuracy.
  • a method of mitigating the effect of a multipath-induced error in a GNSS is performed at a respective GNSS signal receiver.
  • the method includes receiving a band-limited composite signal corresponding to a respective satellite in the GNSS.
  • the band-limited composite signal includes a band-limited direct-path signal and a band-limited multipath signal.
  • the direct-path signal and the multipath signal are modulated with a pseudorandom (PN) code.
  • the method includes obtaining, for a code chip edge transition of the PN code, n pairs of correlation samples of the composite signal during a corresponding n time instances, each having a different time offset from the code chip edge transition.
  • the code chip edge transition has a predetermined filter step response function SR(t).
  • Each respective pair of correlation samples of the n pairs of correlation samples consists of a respective in-phase component sample I(ti) and a respective quadrature component sample Q(ti), wherein i is a positive integer from one to n.
  • the n pairs of correlation samples consist of n in-phase component samples of the composite signal (I(t) samples) and n quadrature component samples of the composite signal (Q(t) samples).
  • the method includes obtaining, for each respective pair of correlation samples of the n pairs of correlation samples, a corresponding sample of the filter step response function SR(ti) corresponding in time with the respective pair of correlation samples, thereby obtaining n samples of the filter step response function (predetermined SR(t) samples) during the corresponding n time instances.
  • the method includes, in accordance with a determination that a computed similarity between the I(t) samples and the SR(t) samples satisfies a first threshold value, determining that a phase ⁇ ⁇ is ⁇ 90° within a first predefined margin, wherein the phase ⁇ ⁇ is equal to 180° minus a sum of (i) a carrier phase multipath error ⁇ ⁇ and (ii) a phase difference ⁇ ⁇ between the direct-path signal and the multipath signal.
  • the method includes, in accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples does not satisfy the first threshold value: solving a first set of matrix equations, thereby obtaining a solution for tan ⁇ ⁇ and for ⁇ ⁇ tan ⁇ ⁇ , wherein ⁇ ⁇ is a response time error of the code chip edge transition due to the multipath signal.
  • the method includes, in accordance with a determination that
  • the method includes, in accordance with a determination that
  • the method includes adjusting a pseudorange measurement for the respective satellite in accordance with the determined ⁇ ⁇ ; and/or adjusting a carrier phase measurement for the respective satellite in accordance with a parameter corresponding to the determined ⁇ ⁇ .
  • the method includes performing a navigation function using the adjusted pseudorange and/or the adjusted carrier phase measurement for the respective satellite. [0007] In some embodiments, the method further includes determining a sign of the carrier phase multipath error ⁇ ⁇ .
  • the method incudes, in accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples satisfies the first threshold value: determining that ⁇ ⁇ is +90° within the first predefined margin when the sign of the carrier phase multipath error ⁇ ⁇ is positive; and determining that ⁇ ⁇ is -90° within the first predefined margin when the sign of the carrier phase multipath error ⁇ ⁇ is negative.
  • the method further includes, in accordance with a determination that the phase ⁇ ⁇ is 0° or 180° within the second predefined margin: determining a position of a peak of the I(t) samples relative to a position of a peak of the predetermined SR(t) samples.
  • the method includes, in accordance with a determination that the position of the peak of the I(t) samples occurs later in time than the position of the peak of the SR(t) samples, determining that ⁇ ⁇ is 0° within the second predefined margin; and in accordance with a determination that the position of the peak of the I(t) samples occurs earlier in time than the position of the peak of the predetermined SR(t) samples, determining that ⁇ ⁇ is 180° within the second predefined margin.
  • the computed similarity between the I(t) samples and the SR(t) samples is obtained using a mean squared error (MSE), defined by ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ( ⁇ ) ⁇ ⁇ ⁇ ⁇ ( ⁇ ) ⁇ , wherein ⁇ ⁇ ⁇ cos ⁇ , ⁇ is a magnitude of the direct-path signal, and ⁇ ⁇ is the carrier phase multipath error.
  • MSE mean squared error
  • solving the first set of matrix equations includes obtaining (i) an initial estimate ⁇ ⁇ ⁇ , ⁇ for the phase ⁇ , (ii) an initial estimate ⁇ ⁇ ⁇ for the response time error ⁇ ⁇ , and (iii) an initial estimate ⁇ ⁇ ⁇ for a magnitude ⁇ of the composite signal; and using the initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , and ⁇ ⁇ ⁇ to solve the first set of matrix equations via a Least Squares fitting process to obtain the solution for tan ⁇ ⁇ and ⁇ ⁇ tan ⁇ ⁇ .
  • the method further includes, in accordance with the determination that
  • the method includes, in accordance with a determination that a number of iterations for which the first iterative process has looped satisfies a third threshold value: outputting information that no valid solution for ⁇ ⁇ has been determined; and terminating the first iterative process without producing a solution for ⁇ ⁇ .
  • the method further includes, in accordance with a determination that a number of iterations for which the first iterative process has looped does not satisfy the third threshold value: in accordance with a determination that a magnitude based on the updated estimate ⁇ ⁇ ⁇ or based on the updated estimate tan ⁇ ⁇ ⁇ , ⁇ satisfies a fourth threshold value: determining that the Least Squares fitting process has converged; outputting the updated estimates tan ⁇ , ⁇ and the updated estimate ⁇ ⁇ as the solutions for tan ⁇ and ⁇ ⁇ , respectively; calculating the phase ⁇ ⁇ from the updated estimate tan ⁇ ⁇ ⁇ , ⁇ ; and terminating the first iterative process.
  • the method also includes, in accordance with a determination that (i) the magnitude based on the updated estimate ⁇ ⁇ ⁇ or based on the updated estimate tan ⁇ ⁇ ⁇ , ⁇ does not satisfy the fourth threshold value, performing a next iteration of the first iterative process.
  • the method further includes, in accordance with a determination that a valid solution for the phase ⁇ ⁇ has been determined: determining (i) an initial estimate ⁇ ⁇ ⁇ , ⁇ for ⁇ , wherein ⁇ ⁇ ⁇ cos ⁇ , ⁇ is a magnitude of the direct-path signal, and ⁇ ⁇ is the carrier phase multipath error; (ii) an initial estimate ⁇ ⁇ ⁇ , ⁇ for a magnitude ⁇ , wherein ⁇ ⁇ ⁇ sin ⁇ ; (iii) an initial estimate ⁇ ⁇ ⁇ , ⁇ for a magnitude ⁇ of the multipath signal, and (iv) an initial estimate ⁇ ⁇ ⁇ for a time delay ⁇ of the multipath signal relative to the direct-path signal, or determine a subset of those initial estimates.
  • the method includes using at least a subset of (i) the initial estimate ⁇ ⁇ ⁇ , ⁇ , (ii) the initial estimate ⁇ ⁇ ⁇ , ⁇ , (iii) the initial estimate ⁇ ⁇ ⁇ , ⁇ , and (iv) the initial estimate ⁇ ⁇ as starting values, performing one or more iterative processes to obtain solutions for at least a subset of ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ ; computing the carrier phase multipath error ⁇ ⁇ in accordance with the obtained solutions for at least the subset of ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ ; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • the method further includes, in accordance with a determination that a valid solution for the phase ⁇ ⁇ has been determined: determining (i) an initial estimate ⁇ ⁇ ⁇ ⁇ , ⁇ for ⁇ ⁇ , wherein ⁇ ⁇ ⁇ ⁇ ⁇ cos ⁇ ⁇ , ⁇ ⁇ is a magnitude of the direct-path signal, and ⁇ ⁇ is the carrier phase multipath error; (ii) an initial estimate ⁇ ⁇ ⁇ , ⁇ for a magnitude ⁇ ⁇ , wherein ⁇ ⁇ ⁇ ⁇ sin ⁇ ⁇ ; (iii) an initial estimate ⁇ ⁇ ⁇ ⁇ , ⁇ for a magnitude ⁇ ⁇ of the multipath signal, and (iv) an initial estimate ⁇ ⁇ ⁇ for a time delay ⁇ of the multipath signal relative to the direct-path signal.
  • the initial estimate ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a last in-phase component sample I(tn) of the I(t) samples; and the initial estimate ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a last quadrature component Q(tn) of the Q(t) samples.
  • the initial estimate ⁇ ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a first in-phase component sample I(t 1 ) of the I(t) samples; and the initial estimate ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a first quadrature component Q(t1) of the Q(t) samples.
  • the method includes, in accordance with the determination that the phase ⁇ ⁇ corresponds to ⁇ 90° within the first predefined margin: assigning the obtained solution as the initial estimate ⁇ ⁇ ⁇ ⁇ ,0.
  • the method includes, in accordance with the determination that the phase ⁇ ⁇ is 0° or 180° within the second predefined margin: assigning the initial estimate ⁇ ⁇ ⁇ , ⁇ as zero. [0018] In some embodiments, the initial estimate ⁇ ⁇ ⁇ , ⁇ is 0.5.
  • the method further includes, in accordance with the determination that the phase ⁇ ⁇ does not correspond to ⁇ 90° within the first predefined margin: using the initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ as starting values to obtain a solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ via a second iterative process, wherein k is a positive integer representing a number of iterations of the second iterative process.
  • a k th iteration of the second iterative process includes: obtaining n estimated in-phase component values ( ⁇ ( ⁇ ) values) for the I(t) samples based on estimates ⁇ , , an ⁇ th ⁇ , ⁇ d ⁇ from a previous (k-1) iteration of the second iterative process; determining, for each estimated in-phase component value ⁇ ⁇ ⁇ ⁇ , ⁇ of the ⁇ ⁇ ⁇ ( ⁇ ) values, a corresponding in-phase component estimation error between (i) the respective estimated in-phase component value ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ and (ii) the respective in-phase component sample I(ti) of the I(t) samples, thereby obtaining n in-phase component estimation errors ( ⁇ ⁇ ( ⁇ ⁇ ) estimation errors); solving a second set of matrix equations based on the ⁇ ( ⁇ ) estimation errors, thereby obtaining updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇
  • the method includes, in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ , ⁇ is within a corresponding first valid range and (ii) the updated estimate ⁇ ⁇ ⁇ is within a corresponding second valid range: in accordance with a determination that the number of iterations ⁇ satisfies a fifth threshold value: in accordance with a determination that a change in value of one or more predetermined first parameters between the k th iteration and the (k-1) th iteration satisfies a respective corresponding threshold: determining that the second iterative process has converged; outputting the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ , ⁇ , and ⁇ , respectively; and terminating the second iterative process.
  • the method includes, in accordance with the determination that the change in value of the one or more predetermined parameters between the k th iteration and the (k-1) th iteration does not satisfy the respective corresponding threshold, performing a next iteration of the second iterative process.
  • the method includes, in accordance with the determination that (i) the number of iterations ⁇ does not satisfy the fifth threshold value, terminating the second iterative process without producing a solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ .
  • the method includes, in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ ⁇ , ⁇ is not within the corresponding first valid range or (ii) the updated estimate ⁇ ⁇ ⁇ is not within the corresponding second valid range: adjusting the initial estimate for the time delay ⁇ from ⁇ ⁇ ⁇ to ⁇ ⁇ ⁇ ′ ; and executing the second iterative process using the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the adjusted initial estimate ⁇ ⁇ ′ as the starting values.
  • the method includes, in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ ⁇ , ⁇ is not within the corresponding first valid range or (ii) the updated estimate ⁇ ⁇ ⁇ is not within the corresponding second valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, repeating the steps of adjusting the initial estimate for the time delay ⁇ and executing the second iterative process; and otherwise, outputting information that no valid solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ has been determined.
  • the method includes, after outputting the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ , ⁇ , and ⁇ , respectively, and in accordance with the determination that the phase ⁇ ⁇ is 0° or 180° within the second predefined margin: assigning a value of zero as the solution for ⁇ ⁇ .
  • the method includes determining the carrier phase multipath error ⁇ ⁇ based on the solutions for ⁇ ⁇ and ⁇ ⁇ ; and correcting a pseudorandom phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • the method includes, after outputting (e.g., at an end of the second iterative process) the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ , ⁇ ⁇ , and ⁇ , respectively, and in accordance with the determination that the phase ⁇ ⁇ is not 0° or 180° within the second predefined margin: starting with the initial estimate for ⁇ ⁇ ⁇ , ⁇ , obtaining a solution for ⁇ ⁇ via a third iterative process, wherein k is a positive integer representing a number of iterations of the third iterative process (e.g., k is set or reset to 1 at beginning of the third iterative process).
  • a k th iteration of the third iterative process includes: obtaining n estimated quadrature component values ( ⁇ ⁇ ⁇ ( ⁇ ) values) for the Q(t) samples based on an initial estimate from a previous (k-1) th iteration of the third iterative process; determining, for each estimated quadrature component value ⁇ ⁇ ⁇ of the ⁇ ⁇ ⁇ ( ⁇ ⁇ ) values, a corresponding quadrature component estimation error ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ between (i) the respective estimated quadrature component value ⁇ ⁇ ⁇ ⁇ , ⁇ and (ii) the respective quadrature component sample Q(t i ), thereby obtaining n quadrature component estimation errors estimation errors); solving a third set of matrix equations based on the estimation errors, thereby obtaining an updated estimate ⁇ ⁇ ⁇ , ⁇ for the k th iteration.
  • the method includes, in accordance with a determination that a change in value of one or more predetermined second parameters between the k th iteration and the (k-1) th iteration satisfies a respective corresponding threshold: determining that the third iterative process has converged; outputting the updated estimate the solution for ⁇ ; and terminating the third iterative process.
  • the method includes, in accordance with a determination that the change in value of the one or more predetermined second parameters between the k th iteration and the (k-1) th iteration does not satisfy a respective corresponding threshold, performing a next iteration of the third iterative process.
  • the method includes, in accordance with the determination that (i) the number of iterations ⁇ does not satisfy the sixth threshold value, terminating the third iterative process without producing a solution for ⁇ ⁇ .
  • the method includes computing the carrier phase multipath error ⁇ ⁇ based on the solutions for ⁇ ⁇ and ⁇ ⁇ ; correcting the pseudorange measurement for the respective satellite in accordance with ⁇ ⁇ ; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • the method includes, in accordance with the determination that the phase ⁇ ⁇ corresponds to ⁇ 90° within the first predefined margin: using the initial estimates as starting values, obtaining a solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ via a fourth iterative process, wherein k is a positive integer representing a number of iterations of the fourth iterative process.
  • a k th iteration of the fourth iterative process includes: obtaining n estimated quadrature component values ( ⁇ ⁇ ⁇ ( ⁇ ) values) for the Q(t) samples based on initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ from a previous (k-1) th iteration of the fourth iterative process; determining, for each estimated quadrature component value ⁇ ⁇ ⁇ ⁇ , ⁇ of the ⁇ ⁇ ⁇ ( ⁇ ) values, a respective quadrature component estimation error ⁇ ⁇ , ⁇ between (i) the respective estimated quadrature component value and (ii) the respective quadrature component sample Q(ti), thereby obtaining n quadrature component estimation errors ( ⁇ ⁇ ( ⁇ ⁇ ) estimation errors); and solving a fourth set of matrix equations based on the estimation errors, thereby obtaining updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇
  • the method includes, in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ ⁇ , ⁇ is within a corresponding third valid range and (ii) the updated estimate ⁇ ⁇ ⁇ ⁇ is within a corresponding fourth valid range: in accordance with a determination that the number of iterations ⁇ satisfies a seventh threshold value: in accordance with a determination that a change in value of one or more predetermined second parameters between the k th iteration and the (k-1) th iteration satisfies a respective corresponding threshold: determining that the fourth iterative process has converged; outputting the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ , respectively; and terminating the fourth iterative process.
  • the method includes, in accordance with the determination that the change in value of the one or more predetermined parameters between the k th iteration and the (k-1) th iteration does not satisfy the respective corresponding threshold, performing a next iteration of the fourth iterative process.
  • the method includes, in accordance with the determination that (i) the number of iterations ⁇ does not satisfy the seventh threshold value, terminating the fourth iterative process without producing a solution for ⁇ ⁇ .
  • the method includes, in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ ⁇ , ⁇ is not within a corresponding third valid range or (ii) the updated estimate ⁇ ⁇ ⁇ ⁇ is not within a corresponding fourth valid range: updating the initial estimate for the time delay ⁇ from ⁇ ⁇ ⁇ to ⁇ ⁇ ⁇ ′ ; and executing the fourth iterative process using the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the updated initial estimate ⁇ ⁇ ⁇ ′ as the starting values.
  • the method includes using the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the adjusted initial estimate ⁇ ⁇ ⁇ ′ as the starting values.
  • the method includes, in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ , ⁇ is not within the corresponding third valid range or (ii) the updated estimate ⁇ ⁇ ⁇ ⁇ is not within the corresponding fourth valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, repeating the steps of updating the initial estimate for the time delay ⁇ and executing the fourth iterative process; and otherwise, outputting information that no valid solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ has been determined.
  • the method includes, in conjunction with outputting the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ , respectively: obtaining a solution for ⁇ ⁇ by solving a matrix equation [0030]
  • the method includes computing the carrier phase multipath error ⁇ ⁇ based on the solution for ⁇ ⁇ and the solution for ⁇ ⁇ ; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • a navigation module, or a global navigation satellite system (GNSS) receiver, or a device for mitigating the effect of a multipath-induced error in a global navigation satellite system (GNSS) includes one or more processors and memory.
  • the memory stores instructions that, when executed by the one or more processors, cause the navigation module, or GNSS receiver, or device, to perform any of the methods disclosed herein.
  • a computer-readable storage medium stores instructions which, when executed by a navigation module, or a global navigation satellite system (GNSS) receiver, or a device, that includes one or more processors, cause the navigation module, or the GNSS receiver, or the device, to perform the any of the methods described herein.
  • GNSS global navigation satellite system
  • Figure 1B illustrates a block diagram of a GNSS receiver in accordance with some embodiments.
  • Figure 1C is a block diagram of a computer system, such as a computer system that is part of a moveable object’s navigation system, according to some embodiments.
  • Figure 2 illustrates a block diagram of a baseband DSP channel in accordance with some embodiments.
  • Figure 3 illustrates an example where both a direct-path signal and a multipath signal from a satellite are received by an antenna of a GNSS receiver, in accordance with some embodiments.
  • Figure 4A illustrates a phasor diagram of the in-phase and quadrature components of a direct-path signal, in accordance with some embodiments.
  • Figure 4B illustrates a phasor diagram of the in-phase and quadrature components of a composite signal, which is composed of a direct-path signal and a multipath signal, in accordance with some embodiments.
  • Figure 5 illustrates the magnitude of a filter step response function ⁇ ⁇ ( ⁇ ) versus time, in accordance with some embodiments.
  • Figure 6 illustrates a correlator arrangement in which 16 pairs of samples are arranged along a code chip transition, and a delayed version of the filter step response function ⁇ ⁇ ( ⁇ ⁇ ⁇ ), in accordance with some embodiments.
  • Figures 7A and 7B are phasor diagrams illustrating multipath scenarios where angle ⁇ ⁇ is equal to +90° and -90°, respectively, in accordance with some embodiments.
  • Figures 8A and 8B are phasor diagrams illustrating multipath scenarios where angle ⁇ ⁇ is equal to 0° and 180°, respectively, in accordance with some embodiments.
  • Figure 9 is a workflow for Procedure 1 of the phase multipath mitigation algorithm, in accordance with some embodiments
  • Figures 10A and 10B show an example of ⁇ ( ⁇ ⁇ ) measurement samples when angle ⁇ ⁇ is equal to 0° and 180°, respectively, in accordance with some embodiments.
  • Figure 11 illustrates a workflow 1100 for Procedure 2 of the phase multipath mitigation algorithm, in accordance with some embodiments.
  • Figure 12 illustrates a workflow for solving the unknown parameters, ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ when ⁇ ⁇ is not ⁇ 90°, in accordance with some embodiments.
  • Figure 13 illustrates a workflow for solving for the unknown parameters ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ when ⁇ ⁇ is equal to ⁇ 90°, in accordance with some embodiments.
  • Figure 14 illustrates a series of different initial estimates ⁇ ⁇ ⁇ ⁇ that can be tried sequentially, in accordance with some embodiments.
  • Figure 15A illustrates a phasor diagram for a very weak multipath scenario with ⁇ ⁇ ⁇ 30°.
  • Figure 15B illustrates a correct solution to angle ⁇ ⁇ even when this multipath scenario is detected as a case of ⁇ ⁇ ⁇ ⁇ 90°.
  • Figures 16A and 16B illustrate a flowchart diagram for a method of mitigating the effect of a multipath-induced error in a GNSS, in accordance with some embodiments.
  • DESCRIPTION OF EMBODIMENTS [0052]
  • Figure 1A is a block diagram illustrating a navigation system 100, according to some embodiments.
  • Navigation system 100 enables a moveable object 110 (e.g., a phone, specialized GNSS receiver, a boat, a truck or other vehicle, a farming appliance, mining appliance, drilling system, etc.) to determine, at any point of time, its current position 112 with respect to a global coordinate system (e.g., a coordinate system for Earth 114).
  • Moveable object 110 is equipped with a satellite receiver (navigation signal receiver 120), typically including or coupled to one or more satellite antennas 140, to receive satellite navigation signals from at least four satellites 115 that are orbiting Earth.
  • the satellite navigation signals received by moveable object 110 are typically global navigation satellite system (GNSS) signals.
  • GNSS global navigation satellite system
  • Moveable object 110 optionally receives satellite orbit correction information and satellite clock correction information (sometimes collectively called “correction information”) for the plurality of satellites.
  • the correction information is typically broadcast by and received from one or more satellites 118 distinct from the GNSS satellites 115, using an antenna 142 and signal receiver 152 (see Figure 1C) distinct from antenna 140 and navigation signal receiver 120 used to receive the satellite navigation signals. However, in some embodiments, the same antenna and receiver are used to receive both satellite navigation signals and correction information [0054]
  • Moveable object 110 determines a position of moveable object 110, using the received satellite navigation signals and, optionally, the received satellite orbit correction information and satellite clock correction information for the plurality of satellites.
  • FIG. 1B illustrates a block diagram of an example GNSS receiver 170 (e.g., navigation signal receiver 120) in accordance with some embodiments.
  • the GNSS receiver 170 receives signals from one or more GNSS satellites (e.g., GNSS satellite 171), processes the signals, and delivers Position, Velocity and Time (PVT) results 180.
  • GNSS satellites e.g., GNSS satellite 171
  • PVT Position, Velocity and Time
  • the GNSS receiver 170 includes an antenna 172, a radio frequency (RF) front end module 174, a baseband digital signal processing (DSP) module 176, and a positioning and navigation module 178.
  • the antenna 172 is the first stage of the GNSS receiver 170 and receives RF signals from all visible GNSS satellites.
  • the antenna 172 is equipped with a built-in amplifier.
  • the GNSS receiver 170 is equipped with multiple antennas or an antenna array. The received RF signals are filtered, amplified, down-converted, and digitized by the RF front end module 174.
  • the RF front end module 174 typically includes one or more stages of down-conversions and filters, to convert the RF signals to different intermediate frequency (IF) bands with designated bandwidths while excluding out-of-band interferences.
  • the final filter in high-precision GNSS receivers typically has a wide bandwidth (e.g., 30 MHz or more).
  • One important characteristic of the final filter of the RF front end module 174 is its step response, which will be described later.
  • Analog-to-digital (A/D) converters in the RF front end module convert received analog signals, e.g., from the final filter, into digital Intermediate Frequency (IF) signals.
  • the RF front end electronics 174 are integrated into a RF Application-Specific Integrated Circuit (ASIC) chip.
  • ASIC Application-Specific Integrated Circuit
  • the RF front end module 174 outputs digital IF signals, which are processed by the baseband DSP module 176.
  • the baseband DSP module 176 acquires and tracks each individual visible satellite (e.g., each satellite for which there is a direct line of sight to the receiver 170), generates range measurements, including code phase measurements (and therefore pseudorange measurements), carrier phase measurements and doppler measurements, and decodes navigation data bits encoded in the signal(s) received from each satellite.
  • phase multipath mitigation is performed by the baseband DSP module 176.
  • the baseband DSP module 176 comprises both hardware and firmware, but the division of these two portions is often not clearly defined and depends on a variety of design considerations.
  • the baseband DSP electronics are integrated into a baseband ASIC chip.
  • the RF front-end electronics and the baseband DSP electronics are integrated into a single-die System-On-Chip (SOC) chip.
  • SOC System-On-Chip
  • the positioning and navigation module 178 is implemented in software executed by a microprocessor, while the firmware portion of the baseband DSP module 176 is implemented in either the same microprocessor as the microprocessor of the positioning and navigation module 178, or in a separate microprocessor.
  • the GNSS receiver 170 includes more than one microprocessor.
  • Figure 1C is a block diagram of computer system 130 in, and used by, a moveable object (e.g., moveable object 110, Figure 1A) to determine the position of the moveable object, according to some embodiments.
  • Computer system 130 typically includes one or more processors (sometimes called CPUs, or hardware processors) 202 for executing programs or instructions; memory 210; one or more communications interfaces 206; and one or more communication buses 205 for interconnecting these components.
  • Computer system 130 optionally includes a user interface 209 comprising a display device 211 and one or more input devices 213 (e.g., one or more of a keyboard, mouse, touch screen, keypad, etc.) coupled to other components of computer system 130 by the one or more communication buses 205.
  • Navigation signal receiver 120 or GNSS receiver 170
  • supplemental receiver(s) 152 are also coupled to other components of computer system 130 by the one or more communication buses 205.
  • the one or more communication buses 205 may include circuitry (sometimes called a chipset) that interconnects and controls communications between system components.
  • Communication interface 206 e.g., a receiver or transceiver
  • Memory 210 includes high-speed random access memory, such as DRAM, SRAM, DDR RAM or other random access solid state memory devices; and may include non-volatile memory, such as one or more magnetic disk storage devices, optical disk storage devices, flash memory devices, or other non-volatile solid state storage devices.
  • Memory 210 optionally includes one or more storage devices remotely located from the CPU(s) 202.
  • Memory 210 or alternately the non-volatile memory device(s) within memory 210, comprises a computer readable storage medium (e.g., a non-transitory computer readable storage medium).
  • memory 210 or the computer readable storage medium of memory 210 stores the following programs, modules and data structures, or a subset thereof: • an operating system 212 that includes procedures for handling various basic system services and for performing hardware dependent tasks; • a communications module 214 that operates in conjunction with communication interface 206 (e.g., a receiver and/or transceiver) to handle communications between moveable object 110 and external systems 160 ( Figure 1); the connection between computer system 130 and external systems 160 may include a communication network 162, such as the internet or a public or proprietary wireless network; • optionally, a user interface module 216 for receiving information from one or more input device 213 of user interface 209, and to convey information to a user of moveable object 110 via one or more display or output devices 211; • a navigation module 218 for determining a position of the moveable object and performing one or more navigation functions (e.g., routing of the mobile object; displaying on a map one or more suggested routes for moving the mobile object from a current location to a specified or target location
  • tracking module 222 samples a GNSS satellite signal using a tracking frequency and code shift initially determined by the acquisition module 220, and then adjusted by the tracking module itself as the tracking module determines code phase corrections for respective channels of the received satellite navigation signals.
  • Operating system 212 and each of the above identified modules and applications correspond to a set of instructions for performing a function described above.
  • the set of instructions can be executed by the one or more processors 202 of computer system 130.
  • the above identified modules, applications or programs i.e., sets of instructions
  • memory 210 stores a subset of the modules and data structures identified above. Furthermore, memory 210 optionally stores additional modules and data structures not described above.
  • Figure 1C is intended more as a functional description of the various features which may be present in a computer system 130 of a moveable object 110 than as a structural schematic of the embodiments described herein. In practice, and as recognized by those of ordinary skill in the art, items shown separately could be combined and some items could be separated. For example, some items shown separately in Figure 1C could be combined into a single module or component, and single items could be implemented using two or more modules or components. The actual number of modules and components, and how features are allocated among them will vary from one implementation to another.
  • CPUs 202 include specialized hardware for performing these and other tasks.
  • some of the operations described above with respect to computer system 130 are performed by navigation signal receiver 120 and/or GNSS receiver 170 rather than computer system 130.
  • the baseband DSP module 176 processes digital IF signals channel by channel (e.g., one channel per satellite being tracked).
  • Figure 2 illustrates a block diagram of a baseband DSP channel 250 in accordance with some embodiments.
  • Each channel may include hardware and software components configured to acquire and track one code signal of one satellite.
  • multi-constellation, multi-frequency receivers there can be tens or hundreds of nearly identical channels running and sharing the baseband DSP hardware resources, for example using time multiplexing.
  • the baseband DSP module 176 receives (e.g., fetches) digital in-phase (I) IF samples 252 and quadrature (Q) IF samples 254 from the RF front-end module 174 of the GNSS receiver 170.
  • the digital I IF samples 252 and Q IF samples 254 are first mixed with the replicas of the carrier signal in carrier removal mixers 256, so the carrier frequency of the received signals is stripped off, and the digital IF signals are down-converted to baseband I and Q samples, respectively.
  • the receiver generated carrier signal is synchronized with the frequency of the received IF signals.
  • Maintaining carrier frequency synchronization is a primary goal of the carrier tracking loop in each channel of the receiver. Furthermore, when the carrier tracking loop of a DSP channel 250 is configured to align the phase of the receiver generated carrier signal with the carrier phase of the received IF signals, the carrier tracking loop is called a phase lock loop or phase-locked loop (PLL).
  • PLL phase-locked loop
  • the I and Q signals output by the carrier removal mixers 256 are then correlated with n replicas of code (e.g., the pseudorandom code corresponding to the satellite for which IF signals are being processed, provided by code phase management unit 272) in n pairs of correlators 258 (e.g., correlator pair 258-1, correlator pair 258-2, ... and correlator pair 258-n).
  • each of the n replicas of the code has a different phase delay.
  • the code phase of one of the n replicas is better aligned with the code phase of the received signals than any of the other replicas, and the correlation of that replica with the received signals has a peak value relative to the other correlations.
  • the strongest correlation is identified and used by the code tracking loop’s discriminator and filter 264, and the carrier tracking loop’s discriminator and filter 266, to adjust the code tracking loop and carrier tracking loop, and to remove the code (sometimes called a spreading spectrum code) from the received I and Q samples.
  • each of the n replicas of the code has a different phase delay.
  • the code phase of one of the n replicas is aligned with the code phase of the received signals, therefore their correlation can achieve a peak value.
  • the goal of Code Tracking Loop in the receiver is to maintain this code phase synchronization or code phase lock.
  • All 2n correlation results are coherently integrated for a pre-detection integration time in integrate and dump circuits 260, and then dumped for further signal detection tracking and/or multipath mitigation (e.g., via phase multipath mitigation unit 262), etc.
  • the length of the pre-detection integration (PDI) time depends on the received signal’s navigation data bit rate or symbol rate and other design considerations. A typical value of the PDI is 1 ms.
  • the n pairs of coherent integration I and Q samples are denoted as ⁇ ( ⁇ ⁇ ) 259-1I and ⁇ ( ⁇ ⁇ ) 259-1Q, ⁇ ( ⁇ ⁇ ) 259-2I and ⁇ ( ⁇ ⁇ ) 259-2Q, ..., and ⁇ ( ⁇ ⁇ ) 259-nI and ⁇ ( ⁇ ⁇ ) 259-nQ.
  • the coherent integration samples can be further coherently or non-coherently integrated for a longer duration in a receiver to achieve a higher signal to noise ratio (SNR).
  • ACF auto- correlation function
  • the code tracking loop’s discriminator and filter 264 can identify the code phase difference between the received signals and the Prompt replica of the code, filter it, and use it to adjust and drive code numerically controlled oscillator (NCO) 268.
  • the output of the code NCO 268 is used by code generator or loader 270 to control the frequency and phase of the Prompt code replica such that it is kept aligned with those of the received signals.
  • the code tracking loop’s discriminator and filter 264 can generate code phase measurements of the received signals. Using time information obtained by combining information obtained from the satellite signals of four or more satellites, the code phase measurements are used to generate pseudorange measurements, which are one type of key satellite range measurement for GNSS positioning (e.g., determining the position of a mobile object, such as mobile object 110, Figure 1A).
  • the carrier tracking loop’s discriminator and filter 266 typically utilizes only the pair of Prompt correlations, discriminates the frequency difference and/or phase difference between the received signals and the receiver generated carrier signal, filters the difference(s), and drives carrier NCO 274 so as to maintain synchronization between the received signals (e.g., from the satellite being tracked by the baseband DSP channel 250) and the carrier signal replica produced by carrier NCO 274 and carrier phase summation unit 276.
  • the carrier tracking loop’s discriminator and filter 266 can also output raw navigation data bits or symbols embedded in the satellite signal being processed by the baseband DSP channel 250.
  • Figure 3 illustrates an example where both a direct-path signal 302 and a multipath signal 304 from a satellite 170 are received by an antenna 172 of a GNSS receiver.
  • LOS line-of-sight
  • multipath multipath signal
  • t represents the time that the composite signal is sampled or measured
  • ⁇ ⁇ is an amplitude of the direct-path signal 302
  • ⁇ ⁇ is an amplitude of the multipath signal 304
  • ⁇ is a time delay of the multipath signal 304 relative to the direct-path signal 302 ( ⁇ is also referred to as “multipath time delay”)
  • ⁇ is a phase delay of the multipath signal 304 relative to the direct-path signal 302 ( ⁇ is also referred to as “multipath phase delay”).
  • the time delay ⁇ is always larger than zero.
  • the multipath signal amplitude ⁇ ⁇ is typically smaller than the direct-path signal amplitude ⁇ ⁇ .
  • the multipath phase delay ⁇ can be any value ranging from -180° to 180°. When ⁇ is equal to 0°, the multipath signal is in- phase with the direct-path signal. When ⁇ is equal to 180°, the multipath signal is out-of- phase with the direct-path signal.
  • ⁇ ⁇ , ⁇ , and ⁇ are the three independent parameters that describe a multipath signal.
  • Equation 1 assumes that a direct-path signal always exists, i.e., ⁇ ⁇ is always larger than 0. If the direct-path signal is blocked (e.g., by trees or buildings), corresponding multipath signals, if any, may be detected and excluded by the GNSS receiver using multipath mitigation mechanisms different from those described in this document. Equation 1 also assumes that there is at most one multipath signal for each direct-path signal. If more than one multipath signal exists, it is assumed that the strongest multipath signal is a dominant signal that is much stronger than all the remaining multipath signals such that all the remaining multipath signals can be omitted or ignored for purposes of mitigating multipath-induced errors in GNSS receiver positioning.
  • Figure 4A illustrates a phasor diagram of the in-phase (I) and quadrature (Q) components of a direct-path signal in accordance with some embodiments.
  • Figure 4A shows that this phasor vector is aligned with the I axis, i.e., the value of Q component is equal to 0, and the magnitude of the vector, ⁇ ⁇ 402, is equal to the value of I component.
  • This pair of I and Q samples is just one of the n pairs of coherent integration I and Q samples as shown in Figure 2.
  • the carrier tracking loop utilizes this pair of I and Q samples to make the phase of the generated carrier signal phase-locked with that of the received signal by driving Q component to zero.
  • Figure 4A is a phasor diagram when the carrier tracking loop is in a steady state of phase lock.
  • the vector magnitude ⁇ ⁇ 402 may not be numerically equal to the direct-path signal magnitude ⁇ ⁇ in equation 1 due to receiver signal processing; however, they only differ by a scale factor since the receiver signal processing can be regarded as a linear system.
  • multipath errors can be introduced to range measurements, for example, pseudorange and carrier phase measurements.
  • Figure 4B illustrates a phasor diagram of I (401) and Q (403) components of a composite signal, i.e., a direct-path signal plus a multipath signal, where ⁇ ⁇ 404 is the magnitude of the direct-path signal phasor vector, ⁇ ⁇ 406 is the magnitude of the multipath signal phasor vector, and ⁇ ⁇ 408 is the magnitude of the composite signal phasor vector.
  • the resulted coherent integration I and Q samples are the linear sum of integration values from the direct-path signal and integration values from the multipath signal, corresponding to the sum of the two vectors in the phasor diagram.
  • the summation vector represented by the dashed line of the composite signal phasor vector, is aligned with (e.g., parallel to) the I axis.
  • the direct-path signal phasor vector would be aligned with the I axis instead, just as shown in Figure 4A. Therefore, the angle, ⁇ ⁇ 410 between the summation vector and the direct-path signal vector is the carrier phase multipath error introduced by the multipath signal.
  • An objective of methods and systems disclosed herein is to estimate the carrier phase multipath error ⁇ ⁇ 410, and correct and mitigate multipath errors in the corresponding carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ 410.
  • Figure 4B shows that the direct-path signal phasor vector, the multipath signal phasor vector, and the composite signal phasor vector form a triangle.
  • Phase difference angle ⁇ ⁇ 412 between the direct-path signal vector and the multipath signal vector is supplementary to the multipath phase delay ⁇ as defined in equation 1.
  • Equation 3 When I and Q samples are taken at code chip edge-up or edge-down transitions, equation 3 needs to be changed into the following form: where ⁇ ⁇ ( ⁇ ) is the magnitude step response of the last stage band pass filter in a receiver’s RF front end 174. Equation 6 contains both the step response ⁇ ⁇ ( ⁇ ) and its delayed version ⁇ ⁇ ( ⁇ ⁇ ⁇ ), and, as discussed below, can be used to estimate the multipath parameters (e.g., ⁇ ⁇ , ⁇ , and/or ⁇ ).
  • a GNSS signal is a sinusoidal carrier signal modulated at least by some code, where the code is a pseudo random sequence of ones and zeros, which are called chips.
  • the signal carrier phase will undergo a 180° phase reversal at a code chip edge transition, either from 0 to 1 or from 1 to 0.
  • GPS L1 Coarse Acquisition (CA) signal has a center carrier frequency of 1575.42 MHz, and its CA code has a chipping rate of 1.023 MHz.
  • the baseband DSP 176 carries on correlation and integration for a PDI period of 1 ms, a received GPS L1 CA signal would be expected to have about 1023 code chip edge transitions during each PDI period of 1 ms.
  • Figure 5 illustrates an example of the magnitude 502 of a filter step response function ⁇ ⁇ ( ⁇ ) versus time, in accordance with some embodiments.
  • the dashed line 508 represents an ideal code chip edge-up transition from 0 to 1 instantaneously at time 0. However, because of the RF front end filtering, or equivalently, because the received signal is filtered and thus bandwidth limited, it is impossible to have such an ideal transition which contains infinitely high frequency components in practice.
  • the solid line 506 represents an example of the actual magnitude of the filter step response function ⁇ ⁇ ( ⁇ ) versus time, where a code chip edge-up transition serves as a step input function. Compared to the corresponding ideal step response (represented by the dashed line 508), the actual step response ⁇ ⁇ ( ⁇ ) shown by the solid line 506 not only has a response time of ⁇ 510 but also has some overshoots before and after the transition.
  • step response characteristics are determined by the type and bandwidth of the final filter in the receiver’s RF front end module 174.
  • a 6-pole Butterworth filter with an IF equivalent bandwidth of 30 MHz has a step response time of about 50 ns and enters steady state after about 150 ns.
  • the code chip edge transition step response model ⁇ ⁇ ( ⁇ ) is determined and stored in memory of the receiver before the phase multipath mitigation algorithm (e.g., as described in Figures 9, 11, 12, and 13) is executed.
  • the baseband DSP module 176 of a GNSS receiver outputs I and Q samples along code chip edge transitions.
  • code phase management unit 272 can properly control and set the phase delays of n pairs of code replicas in such a way that the corresponding n pairs of correlation results, ⁇ ( ⁇ ⁇ ) 259-1I and ⁇ ( ⁇ ⁇ ) 259-1Q, ⁇ ( ⁇ ⁇ ) 259-2I and ⁇ ( ⁇ ⁇ ) 259-2Q, ..., and ⁇ ( ⁇ ⁇ ) 259-nI and ⁇ ( ⁇ ⁇ ) 259-nQ, are along code chip edge transitions.
  • the correlation measurements ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ) are the summed, or equivalently, averaged ones over all code chip edge transitions in one period of correlation duration.
  • only a subset of all these ⁇ pairs of correlators are arranged along code chip edge transitions for the purpose of phase multipath mitigation.
  • the effective correlators for the purpose of phase multipath mitigation are constructed in a different approach, for example, by utilizing multiple channel resources.
  • the number of pairs of correlation samples can be increased to cover a larger portion (e.g., even covering a tail portion) of the PN code sequence so as to accurately measure even long delayed multipath signals.
  • enough correlators are arranged to cover a entire code chip edge transition.
  • some aspects of the present disclosure estimate the unknown signal parameters (e.g., the carrier phase multipath error ⁇ ⁇ , the multipath time delay ⁇ , and the amplitude of the multipath signal ⁇ ⁇ ) in two procedures, also referred to herein as Procedure 1 and Procedure 2, which are described below (e.g., with respect to Figures 9 and 11).
  • the code chip edge transition of a direct-path signal can be distorted, and accordingly the code chip edge transition response time can be changed from ⁇ 510, as shown in Figure 5, to ⁇ + ⁇ ⁇ , where the time error ⁇ ⁇ is unknown but of interest since it largely reveals the code tracking error and code phase measurement error due to the influence of the multipath signal.
  • ⁇ ( ⁇ ) ⁇ ⁇ ⁇ ⁇ ( ⁇ ) sin ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ( ⁇ ⁇ ⁇ ) sin ⁇ ⁇ (11) [0089]
  • filter gains which are usually less than 1.
  • Equation 19 indicates that estimate ⁇ ⁇ ⁇ ⁇ is the sum of the previous estimate ⁇ ⁇ ⁇ and instantaneous solution ⁇ ⁇ , and then this essentially accumulated estimate ⁇ ⁇ ⁇ is used to convert ⁇ ( ⁇ ⁇ ) to ⁇ ( ⁇ ⁇ ) using equation 15.
  • the instantaneous solution ⁇ ⁇ will be very close to zero as the Least Squares fitting process iterates and converges. Once the Least Squares fitting process converges and ⁇ ⁇ is very close to zero, both equations 16 and 17 suggest that small errors in estimate ⁇ ⁇ ⁇ ⁇ should not be a concern, because the product of ⁇ ⁇ ⁇ ⁇ and ⁇ ⁇ will eventually approach to zero. [0095] Not all multipath scenarios can be handled by equation 17. For example, in some scenarios, the angle ⁇ ⁇ is equal or very close to either +90° or -90°, as illustrated in Figure 7A and Figure 7B, respectively.
  • tan ⁇ ⁇ will approach ⁇ infinity, and thus solving the corresponding matrix equation 17 will be numerically impossible.
  • the angle ⁇ ⁇ is equal or very close to either 0° or 180°, as illustrated Figure 8A and Figure 8B, respectively, and thus tan ⁇ has a value that is very close to zero. If tan ⁇ is equal or very close to zero, then all multiplications in the matrix equation 17 will produce values that are zero or close to zero, and the right-hand side ⁇ ( ⁇ ⁇ ) of equation 16 should also be zero in theory, and thus solving matrix equation 17 will be numerically unreliable or meaningless.
  • Equation 20 can be equivalently written into the following matrix form: [0098]
  • Figure 9 is a workflow 900 for Procedure 1 of the phase multipath mitigation algorithm, in accordance with some embodiments.
  • the workflow 900 is performed by a GNSS receiver, such as the GNSS receiver 170 in Figure 1B, or by a computing device (e.g., computer system 130, Figure 1C).
  • a GNSS receiver such as the GNSS receiver 170 in Figure 1B
  • a computing device e.g., computer system 130, Figure 1C.
  • ⁇ e.g., corresponding to ⁇ ( ⁇ ⁇ ) 259-1I and ⁇ ( ⁇ ⁇ ) 259-1Q, ⁇ ( ⁇ ⁇ ) 259-2I and ⁇ ( ⁇ ⁇ ) 259-2Q, ..., and ⁇ ( ⁇ ⁇ ) 259-nI and ⁇ ( ⁇ ⁇ ) 259-
  • the n pairs of ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ) are delivered by integrate and dump circuits 260 as illustrated in Figure 2.
  • the chip edge transition model ⁇ ⁇ ( ⁇ ) or ⁇ ⁇ ( ⁇ ⁇ ) is pre-determined and therefore available.
  • An example of a pre-determined chip edge transition model ⁇ ⁇ ( ⁇ ) or ⁇ ⁇ ( ⁇ ⁇ ) is been illustrated in Figures 5 and 6. As can be seen in these figures, the amplitude of ⁇ ⁇ ( ⁇ ) is normalized to 1.
  • the workflow 900 includes determining (904) a sign of carrier phase multipath error ⁇ ⁇ 410.
  • the sign of ⁇ ⁇ is positive if a triangle formed by the direct-path signal phasor vector, the multipath signal phasor vector, and the composite signal phasor vector is in the first quadrant (e.g., Quadrant 1416 in Figure 4B). Otherwise, the sign of ⁇ ⁇ is negative if the triangle is in the fourth quadrant (e.g., Quadrant 4418 in Figure 4B).
  • ⁇ ⁇ ⁇ 90°
  • ⁇ ⁇ +90° if the sign of ⁇ ⁇ has been determined as positive; otherwise, ⁇ ⁇ is ⁇ 90° if the sign of ⁇ ⁇ is negative.
  • the workflow 900 continues with the path indicated by step 914, to iteratively solve for ⁇ ⁇ and ⁇ ⁇ , using a first iterative process.
  • an iteration count e.g., k
  • an iteration count e.g., k
  • the values of parameters ⁇ ⁇ , ⁇ ⁇ and ⁇ ⁇ are estimated (916) (e.g., initial estimates of these values are set, determined or obtained).
  • ⁇ ⁇ (referred to as “x hat”) denotes an estimate of parameter ⁇ in iteration ⁇ .
  • ⁇ ⁇ is used as the initial estimate for a next iteration (i.e., iteration ⁇ + 1) of the first iterative process.
  • the initial estimate of ⁇ i.e., ⁇ ⁇ ⁇
  • the initial estimate ⁇ ⁇ ⁇ can be assigned using the magnitude of the first pair of I and Q samples, e.g., ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ) for a chip edge- down transition.
  • the absolute value of tan ⁇ ⁇ i.e.,
  • is smaller than the predefined threshold (e.g., Threshold2) (922), this is a case of ⁇ ⁇ 0° or 180° (924).
  • Threshold 2 another predefined threshold
  • the peak of their composited correlations will shift to right (i.e., be delayed) if a multipath signal is in phase with its direct-path signal.
  • the peak of their composited correlations will shift to left (i.e., be advanced) if a multipath signal is 180° out of phase with its direct-path signal.
  • the star markers represent ⁇ ( ⁇ ⁇ ) measurement samples
  • the dot markers represent the code chip edge transition model samples ⁇ ⁇ ( ⁇ ⁇ ).
  • the carrier phase multipath error ⁇ ⁇ is 0.
  • the uncertainty or error is tolerable and does not affect the accuracy of estimating the carrier phase multipath error ⁇ ⁇ .
  • the Least Squares fitting process updates (928) the estimate of tan ⁇ ⁇ and ⁇ ⁇ based on equations 18 and 19 to obtain new estimates tan ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ at this k th iteration, respectively.
  • the workflow 900 then proceeds to determine (930) whether the Least Squares fitting process has converged, diverged or not.
  • workflow 900 determines that the Least Squares fitting process has converged (932), outputs a solution for ⁇ ⁇ (934), and terminates the first iterative process.
  • the change of estimates of tan ⁇ in the last two iterations e.g., ⁇ tan ⁇ , ⁇ ⁇ tan ⁇ , ⁇ ⁇
  • workflow 900 determines that the Least Squares fitting process has converged (932), outputs a solution for ⁇ ⁇ (934), and terminates the iterative process.
  • the workflow 900 determines that the iterative process has diverged (or at least cannot converge anymore) and therefore terminates the iterative process without producing a solution for ⁇ ⁇ (934). Otherwise (e.g., k ⁇ Threshold3 or k ⁇ Threshold3), since the process has yet to converge but has not diverged, the workflow 900 repeats the iterative process after increasing the number of iterations ⁇ by 1 (938).
  • a threshold e.g., k > Threshold3
  • the sign of carrier phase error ⁇ ⁇ which was determined in step 904, is utilized again: if the sign of carrier phase error ⁇ ⁇ is positive while ⁇ ⁇ calculated from arctangent is negative, then ⁇ ⁇ + 180° is re-assigned to the final solution of ⁇ ⁇ ; similarly, if the sign of ⁇ ⁇ is negative while ⁇ ⁇ calculated from arctangent is positive, then ⁇ ⁇ ⁇ 180° is re-assigned to the final ⁇ ⁇ .
  • Procedure 1 has solved parameter ⁇ ⁇ and parameter ⁇ ⁇ , based on the results of the last iteration of Procedure 1 (e.g., using equations 18 and 19), while Procedure 2 (as described in Figures 11, 12, and 13) will continue to solve for the remaining signal parameters, e.g., ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , ⁇ , and the carrier phase error ⁇ ⁇ .
  • Procedure 1 may fail to converge and deliver a no-solution (e.g., no valid solution). If this happens, Procedure 2 will not be executed.
  • Equation 26 has four unknowns ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ . Since Equation 26 is non- linear, Taylor expansion and linearization techniques are applied here.
  • Equation 28 can be rewritten into the following form: where, ⁇ ⁇ ( ⁇ ⁇ ) and ⁇ ⁇ ( ⁇ ⁇ ) are the differences between I & Q measurements and their corresponding estimates, i.e., [00116]
  • equation 30 can be written into the following two matrix equations: [00117] As can be seen, the coefficient matrices in the matrix equations 32 and 33 happen to be same.
  • One approach is to solve equations 32 and 33 separately, while estimates of ⁇ ⁇ , ⁇ ⁇ and ⁇ for equation 32 are updated based on equations 34, 36 and 37, and estimates of ⁇ ⁇ , ⁇ ⁇ and ⁇ for equation 33 are updated based on equations 35, 36 and 37.
  • the final solution ⁇ ⁇ can be a weighted sum of the solution ⁇ ⁇ from equation 32 and the other ⁇ ⁇ from equation 33.
  • the final solution ⁇ can also be a weighted sum of the two ⁇ solutions from equations 32 and 33.
  • Equations 32 and 33 Another approach involves converting/combining the two ⁇ ⁇ 3 matrix equations (equations 32 and 33) into one ( 2 ⁇ ) ⁇ 4 matrix equation with a total of four unknowns (e.g., unknowns ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ ⁇ ), and solving the ( 2 ⁇ ) ⁇ 4 matrix equation (e.g., instead of combining solutions).
  • equation 32 is solved to obtain solutions ⁇ ⁇ , ⁇ ⁇ and ⁇ . Because the unknowns ⁇ ⁇ and ⁇ have been solved, their values are now known.
  • Equation 35 can still be used to update ⁇ ⁇ ⁇ , ⁇ iteratively.
  • the third approach solves matrix equation 41 (instead of matrix equation 33) to obtain ⁇ ⁇ .
  • the third approach offers two advantages.
  • the third approach solves a ⁇ ⁇ 3 matrix (equation 32) and a ⁇ ⁇ 1 matrix (equation 41), thus reducing the amount of computational power used to produce a solution.
  • an ⁇ ⁇ 1 matrix equation e.g., equation 41
  • an ⁇ ⁇ 1 matrix equation typically converges in fewer iterations (e.g., two or three iterations). Therefore, solutions for the unknown parameters can be obtained more quickly and efficiently.
  • Figure 11 illustrates a workflow 1100 for Procedure 2 of the phase multipath mitigation algorithm, in accordance with some embodiments.
  • the workflow 1100 is performed by a GNSS receiver (e.g., navigation signal receiver 120 or GNSS receiver 170) or by a computing device (e.g., computer system 130, Figure 1C).
  • a GNSS receiver e.g., navigation signal receiver 120 or GNSS receiver 170
  • a computing device e.g., computer system 130, Figure 1C
  • the steps in the workflow 1100 are executed after Procedure 1 (described in the workflow 900) outputs a valid estimation for the carrier phase multipath error ⁇ ⁇ .
  • the workflow 1100 includes determining (1102) initial estimates, ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ , of four unknown parameters, ⁇ , ⁇ , ⁇ , and ⁇ .
  • the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ can be assigned with the last pair of correlation measurements, i.e., ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ).
  • Procedure 1 e.g., the workflow 900
  • ⁇ ⁇ ⁇ ⁇ , ⁇ can assigned to the value of ⁇ ⁇ that has been solved for from equation 22.
  • the final solution of ⁇ ⁇ is also equal to this value.
  • ⁇ ⁇ ⁇ , ⁇ can assigned to 0.
  • the final solution of ⁇ ⁇ in this case is also 0.
  • the initial multipath signal amplitude estimate ⁇ ⁇ ⁇ ⁇ , ⁇ can be usually set at 0.5, e.g., half of the amplitude which is used to normalize the chip edge transition model ⁇ ⁇ ( ⁇ ) and correlation samples ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ).
  • the estimate ⁇ ⁇ ⁇ ⁇ , ⁇ can be scaled down accordingly.
  • Equation 26 shows that ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ) measurement models are linear with parameters ⁇ ⁇ , ⁇ ⁇ and ⁇ ⁇ , so these initial estimates typically work well. [00131] Equation 26 shows that both ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ ) measurement models are highly nonlinear functions in terms of the multipath time delay parameter ⁇ , and therefore it is critical to make its initial estimate ⁇ ⁇ ⁇ as accurate as possible.
  • the initial multipath time delay estimate ⁇ ⁇ (e.g., ⁇ ⁇ ⁇ 1010 in Figure 10A and ⁇ ⁇ 1020 in Figure 10B) can be set as [00132] If later the Least Squares fitting process cannot converge using this estimate ⁇ ⁇ ⁇ , it suggests this estimate ⁇ ⁇ ⁇ is probably too erroneous.
  • Figures 10A and 10B are examples of ⁇ ( ⁇ ⁇ ) measurement samples when angle ⁇ ⁇ is equal to 0° and 180°, respectively, this method of estimating ⁇ is valid for all multipath scenarios.
  • Procedure 2 Having determined the initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ , Procedure 2 selects an appropriate algorithm to iteratively solve for the unknown parameters, ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ , depending on whether ⁇ ⁇ is equal to ⁇ 90° or not (1104).
  • FIG. 12 illustrates a workflow 1200 for solving the unknown parameters, ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ when ⁇ ⁇ is not ⁇ 90°, in accordance with some embodiments.
  • an iterative process at iteration k first estimates (1204) I measurements, ⁇ ⁇ ⁇ ⁇ , ⁇ , based on initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ using equation 29, calculates its estimation error using equation 31, constructs I-matrix equation 32, and solves for unknowns ⁇ ⁇ , ⁇ ⁇ and ⁇ ⁇ using the Least-Squares method (1206).
  • the workflow 1200 updates (1208) ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ using equations 34, 36 and 37.
  • a post-fix ⁇ ⁇ ⁇ ⁇ , ⁇ at each iteration is also calculated as follows: [00135]
  • the workflow 1200 includes determining (1210) whether the iterative process has diverged or the number of iterations, ⁇ , is larger than a threshold (e.g., Threshold4).
  • a threshold e.g., Threshold4
  • the workflow 1200 determines that the iterative process has diverged (1224). This is because, since a multipath signal must travel through a longer path than its direct-path signal, the multipath time delay ⁇ should always be larger than 0. However, if ⁇ ⁇ ⁇ is larger than ⁇ , it means that the multipath signal effect is not observable even by the last pair of correlators.
  • ⁇ ⁇ ⁇ ⁇ , ⁇ is less than 0, it means multipath signal is either too weak or does not exist.
  • These thresholds may need some fine-tuning or margin, since, for example, ⁇ ⁇ ⁇ ⁇ may oscillate to a slightly negative value before converging to a small positive one.
  • the way of determining convergence here is similar to the one in Procedure 1: if changes in the estimates between two latest iterations are smaller than the corresponding thresholds, and/or if the change in ⁇ ⁇ ⁇ ⁇ , i.e., ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ , is smaller than a threshold, and/or if ⁇ ⁇ ⁇ ⁇ , ⁇ itself is smaller than a threshold, then the workflow 1200 determines that the iterative process has converged (1216); otherwise, the iterative process has not converged yet (1220).
  • Figure 12 also illustrates that in some instances, if the iterative process has diverged or if ⁇ is larger than a threshold (1224), then the current iteration is terminated and a no-solution (e.g., no valid solution) is output (1231). In some instances, divergence may be caused by a bad initial estimate ⁇ ⁇ .
  • a new, different initial estimate ⁇ ⁇ ⁇ is tried, as illustrated in steps 1226 and 1228.
  • a series of different initial estimates ⁇ ⁇ ⁇ ⁇ can be tried sequentially (e.g., in a “tree pattern”). This is illustrated in Figure 14. If the first estimate ⁇ ⁇ ⁇ 1402 results in a divergence, a second new estimate ⁇ ⁇ ′ 1404 can be adjusted to its right. If the second estimate ⁇ ⁇ ⁇ ′ 1404 also fails, then a third new estimate ⁇ ⁇ ⁇ ′′ 1406 with a value smaller than the first estimate ⁇ ⁇ ⁇ 1402 can be tried.
  • the adjustment step size can be a constant (e.g., equal to the correlator spacing).
  • the range for the adjustments say, from ⁇ ⁇ 1420 to ⁇ ⁇ 1422, should cover all correlation sample times, i.e., from 1424 to ⁇ ⁇ 1426, with some margin. Once an adjustment moves beyond this range, it means that all possible initial estimates in this corresponding direction have been tried.
  • an iteration count e.g., k, representing a number of iterations of the third iterative process is set or reset to an initial value, e.g., 1.
  • the third iterative process calculates (1240) Q measurements, ⁇ ⁇ ⁇ ⁇ , ⁇ , based on the initial estimate using equation 39, calculates (1240) the measurement estimation error ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ using equation 31, constructs Q-matrix equation 41, solves (1242) for unknowns ⁇ ⁇ using the Least-Squares method, and updates (1242) ⁇ ⁇ ⁇ , ⁇ using equation 35. Meanwhile, a post-fix ⁇ ⁇ ⁇ , ⁇ at each iteration is also calculated as follows: [00141] At the end of each iteration (e.g., iteration k), the workflow 1200 determines (1244) whether convergence of the third iterative process has occurred.
  • a way of checking convergence is similar to that described above with respect to step 1214: if the change of estimates of ⁇ in the last two iterations, e.g., is smaller than a threshold, and/or if the change in ⁇ ⁇ ⁇ ⁇ , i.e., ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ , is smaller than a threshold, and/or if ⁇ ⁇ ⁇ ⁇ , ⁇ itself is smaller than a threshold, then the third iterative process has successfully converged (1246), and thus ⁇ ⁇ is obtained (1248).
  • the workflow 1200 in Figure 12 describes the steps of solving for ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ when ⁇ ⁇ is not equal to ⁇ 90°.
  • the iterative process for determining parameters ⁇ ⁇ , ⁇ ⁇ and ⁇ , or the iterative process (e.g., the third iterative process) for determining ⁇ ⁇ can be stopped (e.g., when divergence occurs or when the numbers of iterations has exceeded a threshold number, such as Threshold4 or Threshold5) and information indicating that no valid solution has been determined is output.
  • Figure 13 illustrates a workflow 1300 for solving for the unknown parameters ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ when ⁇ ⁇ is equal to ⁇ 90° (step 1302, Figure 13), in accordance with some embodiments.
  • Procedure 1 (workflow 900) has determined that ⁇ ⁇ is equal to ⁇ 90° using ⁇ ( ⁇ ⁇ ) measurements based on equation 22, which is used to solve for ⁇ ⁇ , there is no need to repeat this process of solving for some unknowns from ⁇ ( ⁇ ⁇ ) measurements. Instead, the solution ⁇ ⁇ to equation 22 obtained in Procedure 1 can be assigned not only to the initial estimate but also to the final solution of ⁇ ⁇ .
  • the remaining steps in the workflow 1300 solve for ⁇ ⁇ , ⁇ ⁇ and ⁇ using matrix equation 33.
  • the steps for solving matrix equation 33 which are shown in the workflow 1300, are similar to those for solving matrix equation 32, which are shown on the left portion of the workflow 1200.
  • the workflow 1300 in Figure 13 is an iterative process (e.g., a fourth iterative process) that includes calculating (1304) (e.g., at iteration k) Q measurements, based on initial estimates ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ using equation 29, calculating (1304) the estimation error ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ using equation 31, constructing Q-matrix equation 33, and solving (1306) for unknowns ⁇ ⁇ , ⁇ ⁇ and ⁇ ⁇ using the Least-Squares method.
  • calculating (1304) e.g., at iteration k
  • Q measurements based on initial estimates ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ using equation 29, calculating (1304) the estimation error ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ using equation 31, constructing Q-matrix equation 33, and solving (1306) for unknowns ⁇ ⁇ , ⁇ ⁇ and ⁇ ⁇ using the Least
  • an iteration count e.g., k, representing a number of iterations of the four iterative process is set or reset to an initial value, e.g., 1.
  • the workflow 1300 includes updating (1308) ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ using equations 35, 36 and 37. Meanwhile, a post-fix ⁇ ⁇ ⁇ , ⁇ at each iteration is also calculated using equation 44.
  • the workflow 1300 determines (1310) whether the iterative process has diverged or the number of iterations, ⁇ , is larger than a threshold (e.g., Threshold6). If estimates are out of their valid range (e.g., if estimate ⁇ ⁇ ⁇ is larger than the allowably maximal time delay corresponding to the last pair of correlators, ⁇ , or if ⁇ ⁇ ⁇ is less than 0, or if estimate ⁇ ⁇ ⁇ , ⁇ is less than 0, then the iterative process is determined to have diverged (1320).
  • a threshold e.g., Threshold6
  • these thresholds may be fine-tuned to some margin, since ⁇ ⁇ ⁇ may oscillate to a slightly negative value before converging to a small positive one.
  • the workflow 1300 determines that the fourth iterative process has converged (1316). Otherwise, the workflow 1300 determines that the fourth iterative process has not converged. If the fourth iterative process has converged (1316), the Q-matrix equation 33 has been successfully solved, and the ⁇ ⁇ , ⁇ ⁇ and ⁇ are obtained (1318).
  • the workflow 1300 includes re-solving equation 33 using a new, different estimate ⁇ ⁇ ⁇ (1322). For example, in some embodiments, Procedure 2 selects different initial estimates ⁇ ⁇ ⁇ ⁇ using the approach that is explained above with respect to Figure 14.
  • Procedure 2 has solved for parameters ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ , regardless of whether ⁇ ⁇ is equal to ⁇ 90° or not. It is also possible that no valid solution is obtained in Procedure 1 and/or Procedure 2 (steps 110, 1112 and 1114). In some embodiments, if there is no valid solution, this means that no error correction to carrier or code phase measurements is available.
  • the phase multipath mitigation techniques output information indicating that there is not a valid solution, instead of a suspiciously bad solution, to protect GNSS receivers’ reliability performance.
  • angles ⁇ ⁇ , ⁇ ⁇ and ⁇ ⁇ are all positive if a triangle formed by three phasor vectors is in Quadrant 1 416, as illustrated in Figure 4B; otherwise, they are all negative if the triangle is in Quadrant 4418.
  • all of the unknown parameters including ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ ⁇ , and ⁇ , have been solved.
  • the carrier phase measurement error ⁇ ⁇ can be used to correct the corresponding carrier phase measurement, and thus the carrier phase multipath error is mitigated.
  • the time error ⁇ ⁇ if solved in Procedure 1, can be used as a reference to correct the corresponding code phase measurement and pseudorange measurement. In some embodiments, these parameters, whose values have been determined, can be further used, for example, to mitigate other errors using other signal tracking and multipath mitigation techniques. [00154] Some applications solve for the carrier phase multipath error ⁇ ⁇ and/or the time error ⁇ ⁇ . In these instances, the steps (described above) that are used to calculate other parameters can be skipped as soon as ⁇ ⁇ and ⁇ ⁇ are obtained.
  • Procedure 1 determines that ⁇ ⁇ is equal to 0° or 180°, then the entire Procedure 2 can be skipped and a correct solution to ⁇ ⁇ can be directly delivered ( ⁇ ⁇ is equal to 0° in this case).
  • ⁇ ⁇ is equal to 0° in this case.
  • computational power and time can be saved.
  • the geometry in Figure 15A is assumed to reflect a weak multipath scenario, where the multipath signal amplitude ⁇ ⁇ ⁇ ⁇ ⁇ , and the value of the angle ⁇ ⁇ (“true value”) is equal to 30°.
  • ⁇ ⁇ is detected to be 90° by Procedure 1, and then the steps described in Figure 13 are executed, where low SNR measurements ⁇ ( ⁇ ⁇ ) are utilized, and these nearly zero-valued measurements ⁇ ( ⁇ ⁇ ) result in small values of solutions ⁇ ⁇ and ⁇ ⁇ .
  • the comparison between Figure 15A and Figure 15B demonstrates that the solution for the carrier phase multipath error ⁇ ⁇ in Figure 15B is quite close to the ⁇ ⁇ carrier phase multipath error in Figure 15A, even though the solution ⁇ ⁇ could be very different from the “true value” of ⁇ ⁇ .
  • n data marks 602-i in Figure 6 (representing n pairs of correlation measurements) only cover a portion of the transition, this will result in an ambiguity between multipath time delay ⁇ and multipath signal amplitude ⁇ ⁇ . In other words, this phase multipath mitigation may converge to a wrong solution or even diverge. Therefore, it is recommended that n pairs of correlation measurements should be arranged along the whole code chip edge transition including its tail. [00160] In addition to a code chip edge transition model ⁇ ⁇ ( ⁇ ), Figure 6 also illustrates its time-delayed version ⁇ ⁇ ( ⁇ ⁇ ⁇ ).
  • Figure 6 shows that ⁇ ⁇ ( ⁇ ⁇ ⁇ ) and the corresponding multipath signal may not be observable even by the last pair of receiver correlators if multipath time delay ⁇ is too large, larger than the last correlator pair’s code phase ⁇ ⁇ .
  • more correlator resources are needed in receiver hardware design to cover and measure the very tail part of the code chip edge transition model ⁇ ⁇ ( ⁇ ), if correlator spacing remains unchanged.
  • phase multipath mitigation system needs to deal with very short- and very long-delayed multipath scenarios, it demands more correlator resources. This tradeoff between receiver hardware complexity and multipath mitigation performance should be considered in GNSS receiver design. In some instances. when the receiver hardware does not have enough resources to support multipath mitigation for some multipath scenarios, the phase multipath mitigation algorithm that is described in Procedure 1 and Procedure 2 should determine the MSE values for most if not all steps, as described above, and deliver a no- solution as needed, to ensure that the receiver’s reliability and integrity performance are protected.
  • FIGS 16A and 16B illustrate a flowchart diagram for a method 1600 of mitigating the effect of a multipath-induced error in a global navigation satellite system (GNSS), in accordance with some embodiments.
  • the method 1600 is performed at a GNSS receiver (e.g., navigation signal receiver 120 in Figure 1A or GNSS receiver 170 in Figure 1B), or at a computing device (e.g., computer system 130) that includes one or more processors (e.g., CPU(s) 202) and memory (e.g., memory 210).
  • a GNSS receiver e.g., navigation signal receiver 120 in Figure 1A or GNSS receiver 170 in Figure 1B
  • a computing device e.g., computer system 130
  • processors e.g., CPU(s) 202
  • memory e.g., memory 210
  • the GNSS receiver receives (1602) a band-limited composite signal corresponding to a respective satellite in the GNSS, including a band-limited direct-path signal and a band-limited multipath signal.
  • the direct-path signal and the multipath signal are modulated with a pseudorandom (PN) code
  • PN pseudorandom
  • the GNSS receiver obtains (1604), for a code chip edge transition of the PN code (e.g., a chip edge-up transition or a chip edge-down transition), n pairs of correlation samples (e.g., samples 602-i, Figure 6) of the composite signal during a corresponding n time instances, each having a different time offset from the code chip edge transition.
  • n is a positive integer such as 16, 32, 48, or the like.
  • the code chip edge transition has a predetermined filter step response function SR(t) (see, e.g., Figure 5).
  • Each respective pair of correlation samples of the n pairs of correlation samples consists of a respective in-phase component sample I(ti) and a respective quadrature component sample Q(ti), wherein i is a positive integer from one to n.
  • the n pairs of correlation samples consist of n in-phase component samples of the composite signal (I(t) samples) and n quadrature component samples of the composite signal (Q(t) samples).
  • the n pairs of correlation samples are uniformly arranged along the code chip edge transition.
  • the n pairs of correlation are arranged over one period of correlation duration.
  • a respective pair of correlation samples (I(ti), Q(ti)) is summed or averaged over multiple periods of correlation duration.
  • the number of the pairs of correlation samples can be increased to cover a larger portion (e.g., even covering a tail portion) of the PN code sequence so as to accurately measure even long delayed multipath.
  • the GNSS receiver obtains (1606), for each respective pair of correlation samples of the n pairs of correlation samples, a corresponding sample of the filter step response function SR(ti) corresponding in time with the respective pair of correlation samples, thereby obtaining n samples of the filter step response function (predetermined SR(t) samples) during the corresponding n time instances.
  • the GNSS receiver determines (1608) that a phase ⁇ ⁇ is ⁇ 90° within a first predefined margin (e.g., no greater than 1, 2, 3, 4 or 5 degrees).
  • the phase ⁇ ⁇ is equal to 180° minus a sum of (i) a carrier phase multipath error ⁇ ⁇ (e.g., which can be used to correct carrier phase measurement) and (ii) a phase difference ⁇ ⁇ between the direct-path signal and the multipath signal (see Equation 2).
  • ⁇ ⁇ is equal to 180 degrees minus the difference between the direct-path signal and the multipath signal.
  • the GNSS receiver solves (1610) a first set of matrix equations (e.g., Equation 17), thereby obtaining a solution for tan ⁇ ⁇ and for ⁇ ⁇ tan ⁇ ⁇ , wherein ⁇ ⁇ is a response time error of the code chip edge transition due to the multipath signal.
  • Equation 17 a first set of matrix equations
  • the GNSS receiver determines (1612) that the phase ⁇ ⁇ is 0° or 180° within a second predefined margin (e.g., no greater 1, 0.5, 0.1, 0.01, or 0.001 degrees).
  • the second threshold value e.g., Threshold2 in Figure 9, step 920
  • the second threshold value is a value such as 0.01 or less, with some examples being 0.0001, 0.001, 0.005, or 0.01.
  • the GNSS receiver determines (1614) ⁇ ⁇ in accordance with the solution for tan ⁇ ⁇ (see above discussion of Figure 9, steps subsequent to step 926). For example, in some embodiments,
  • the GNSS receiver adjusts (1616) a pseudorange measurement for the respective satellite in accordance with the determined ⁇ ⁇ ; and/or adjusts a carrier phase measurement for the respective satellite in accordance with a parameter (e.g., the carrier phase multipath error ⁇ ⁇ ) corresponding to the determined ⁇ ⁇ .
  • the GNSS receiver performs (1618) a navigation function (e.g., determines a location of the GNSS receiver) using the adjusted pseudorange and/or carrier phase measurements for the respective satellite.
  • the GNSS receiver determines, based on the Q(t) samples, a sign (e.g., positive or negative) of the carrier phase multipath error ⁇ ⁇ .
  • the phasor diagram is in the first quadrant (e.g., Quadrant 1416, Figure 4B) and the sign of the carrier phase multipath error ⁇ ⁇ is positive.
  • the phasor diagram is in the fourth quadrant (e.g., Quadrant 4418) and the sign of the carrier phase multipath error ⁇ ⁇ is negative.
  • some ⁇ ( ⁇ ⁇ ) are positive while some are negative.
  • the GNSS receiver finds the maximum of their magnitudes, say
  • the GNSS receiver determines that ⁇ ⁇ is +90° within the first predefined margin when the sign of the carrier phase multipath error ⁇ ⁇ is positive, and determines that ⁇ ⁇ is -90° within the first predefined margin when the sign of the carrier phase multipath error ⁇ ⁇ is negative.
  • the GNSS receiver determines a position of a peak of the I(t) samples relative to a position of a peak of the predetermined SR(t) samples. In accordance with a determination that the position of the peak of the I(t) samples (e.g., ⁇ ⁇ , ⁇ 1004 in Figure 10A) occurs later in time than the position of the peak of the SR(t) samples (e.g., ⁇ ⁇ , ⁇ 1008 in Figure 10A), the GNSS receiver determines that ⁇ ⁇ is 0° within the second predefined margin.
  • MSE mean squared error
  • ⁇ ⁇ ⁇ cos ⁇ ⁇ , where ⁇ ⁇ is a magnitude of the direct-path signal, and ⁇ ⁇ is the carrier phase multipath error.
  • solving the first set of matrix equations comprises: obtaining (i) an initial estimate ⁇ ⁇ ⁇ , ⁇ for the phase ⁇ , (ii) an initial estimate ⁇ ⁇ ⁇ for the response time error ⁇ ⁇ , and (iii) an initial estimate ⁇ ⁇ ⁇ for a magnitude ⁇ of the composite signal (See Figure 9, step 916).
  • the initial estimate ⁇ ⁇ ⁇ , ⁇ is set to zero.
  • the initial estimate ⁇ ⁇ ⁇ is set to zero.
  • the GNSS receiver uses the initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , and ⁇ ⁇ ⁇ to solve the first set of matrix equations (e.g., Equation 17) via a Least Squares fitting process to obtain the solution for tan ⁇ ⁇ and ⁇ ⁇ tan ⁇ ⁇ .
  • the GNSS receiver in accordance with the determination that
  • Each iteration of the first iterative process includes: obtaining an updated estimate tan ⁇ ⁇ ⁇ , ⁇ for the phase ⁇ ⁇ and an updated estimate ⁇ ⁇ ⁇ ⁇ for the response time error ⁇ ⁇ (e.g., according to Equations 18 and 19) (see, e.g., Figure 9, step 928).
  • a number of iterations e.g., the current number of iterations, k
  • a third threshold value e.g., k > Threshold3 or k ⁇ Threshold3 as illustrated in Figure 9, steps 930 and 932: the GNSS receiver outputs information that no valid solution for ⁇ ⁇ has been determined (see, e.g., Figure 9, step 934) and terminates the first iterative process without producing a solution for ⁇ ⁇ .
  • the first iterative process has diverged (or at least cannot converge anymore) and therefore should be cancelled with a no valid solution.
  • each iteration of the first iterative process further includes: in accordance with a determination that the number of iterations for the first iterative process does not satisfy the third threshold value (e.g., k ⁇ Threshold3 or k ⁇ Threshold3) (e.g., the first iterative process has not diverged, or the number of iterations, ⁇ , for the least squares fitting process does not equal or exceed a maximum value): in accordance with a determination that a magnitude based on the updated estimate ⁇ ⁇ ⁇ (e.g., ⁇ ⁇ ⁇ ) or based on the updated estimate tan ⁇ ⁇ ⁇ , ⁇ (e.g., a difference in magnitude between the updated estimate tan ⁇ , ⁇ and the estimate tan ⁇ , ⁇ from the previous iteration, ⁇ ) satisfies (e.g., is less than, or is less than or equal to) a fourth threshold value: determining that the Least Squares fitting process has
  • each iteration of the first iterative process further includes: in accordance with a determination that the magnitude based on the updated estimate ⁇ ⁇ ⁇ ⁇ or based on the updated estimate tan ⁇ ⁇ ⁇ , ⁇ does not satisfy the fourth threshold value, performing a next iteration of the first iterative process (see., e.g., Figure 9, step 938).
  • the magnitude based on the updated estimate ⁇ ⁇ ⁇ ⁇ or based on the updated estimate tan ⁇ ⁇ ⁇ , ⁇ does not satisfy the fourth threshold value when ⁇ ⁇ ⁇ is greater than, or is greater than or equal to, the fourth threshold value, or when is greater than, or is greater than or equal to, the fourth threshold value.
  • the GNSS receiver uses at least a subset of (i) the initial estimate ⁇ ⁇ ⁇ , ⁇ , (ii) the initial estimate ⁇ ⁇ ⁇ , ⁇ , (iii) the initial estimate ⁇ ⁇ ⁇ , ⁇ , and (iv) the initial estimate ⁇ ⁇ as starting values, performing one or more iterative processes to obtain solutions for at least a subset of ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ (e.g. using Procedure 2, as illustrated in Figures 11, 12, and 13).
  • the GNSS receiver computes the carrier phase multipath error ⁇ ⁇ in accordance with the obtained solutions for at least the subset of ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ , and corrects the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • the initial estimate ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a last in-phase component sample I(tn) of the I(t) samples; and the initial estimate ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a last quadrature component Q(tn) of the Q(t) samples.
  • the initial estimate ⁇ ⁇ ⁇ ⁇ , ⁇ corresponds to a magnitude of a first in-phase component sample I(t 1 ) of the I(t) samples; and the initial estimate ⁇ , ⁇ corresponds to a magnitude of a first quadrature component Q(t1) of the Q(t) samples.
  • the GNSS receiver in accordance with the determination that the phase ⁇ ⁇ corresponds to ⁇ 90° within the first predefined margin, the GNSS receiver obtains a solution for ⁇ ⁇ by solving a matrix equation (see Equation 22), and assigning the obtained solution as the initial estimate ⁇ ⁇ ⁇ , ⁇ . [00185] In some embodiments, in accordance with the determination that the phase ⁇ ⁇ is 0° or 180° within the second predefined margin, the GNSS receiver assigns the initial estimate ⁇ ⁇ ⁇ , ⁇ as zero.
  • the initial estimate ⁇ ⁇ ⁇ , ⁇ is 0.5 (e.g., half of the amplitude that is used to normalize the chip edge transition model ⁇ ⁇ ( ⁇ ) and correlation samples ⁇ ( ⁇ ⁇ ) and ⁇ ( ⁇ ⁇ )).
  • estimate ⁇ ⁇ ⁇ , ⁇ can be scaled down accordingly.
  • the GNSS receiver uses the initial estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as starting values (e.g., obtained using equation 29) to obtain a solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ via a second iterative process (see, e.g., workflow 1200 in Figure 12).
  • an iteration count e.g., k, representing a number of iterations of the second iterative process is set or reset to an initial value, e.g., 1.
  • the one or more predetermined first parameters include one or more of ⁇ ⁇ , ⁇ ⁇ , and ⁇ .
  • the determination that the second iterative process has converged is in accordance with a determination that (i) the change in an estimated value of ⁇ ⁇ between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e., ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ ⁇ ), and/or (ii) the change in an estimated value of ⁇ ⁇ between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e., ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ h ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ ⁇ ), and/or (ii) the change in an estimated value of
  • the one or more predetermined first parameters include a post-fix mean squared error (MSE) for the in-phase component estimation error, defined by the k th iteration (see Equation 43).
  • MSE post-fix mean squared error
  • the determining that the second iterative process has converged is in accordance with a determination that a change in the post-fix MSE for the in-phase component estimation error between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ , ⁇ )).
  • a k th iteration of the second iterative process further includes, in accordance with the determination that the change in value of the one or more predetermined parameters between the k th iteration and the (k-1) th iteration does not satisfy the respective corresponding threshold (e.g., the second iterative process has not converged), performing a next iteration of the second iterative process (see, e.g., Figure 12, steps 1220 and 1222).
  • the change in value of the one or more predetermined parameters between the k th iteration and the (k-1) th iteration does not satisfy the respective corresponding threshold when or ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ , ⁇ .
  • a k th iteration of the second iterative process further includes: in accordance with the determination that (i) the number of iterations ⁇ does not satisfy the fifth threshold value (e.g., k does not satisfy the fifth threshold value when k > Threshold4 or k ⁇ Threshold4, as illustrated in Figure 12, step 1210), terminating the second iterative process without producing a solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ (see, e.g., Figure 12, step 1231).
  • the fifth threshold value e.g., k does not satisfy the fifth threshold value when k > Threshold4 or k ⁇ Threshold4, as illustrated in Figure 12, step 1210
  • the GNSS receiver in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ , ⁇ is not within the corresponding first valid range or (ii) the updated estimate ⁇ ⁇ ⁇ is not within the corresponding second valid range (i.e., the second iterative process has diverged), the GNSS receiver adjusts the initial estimate for the time delay ⁇ from ⁇ ⁇ to ⁇ ⁇ ⁇ ′ (see, e.g., Figure 12, step 1226). The GNSS receiver executes the second iterative process using the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the adjusted initial estimate ⁇ ⁇ ⁇ ′ as the starting values.
  • ⁇ ⁇ can be tried using the pattern illustrated in Figure 14. For example, if the first estimate 1402 ( ⁇ ⁇ ⁇ ⁇ ) results in a divergence, a second new estimate 1404 ( ⁇ ⁇ ⁇ ′ ) can be adjusted to a value that is larger than (e.g., to the right of) that of the first estimate 1402 ( ⁇ ⁇ ⁇ ⁇ ). If the second new estimate 1404 also fails, then a third new estimate 1406 ( ⁇ ⁇ ⁇ ⁇ ) can be adjusted to a value that is smaller than (e.g., to left of) that of the first estimate 1402 ( ⁇ ⁇ ⁇ ⁇ ).
  • a fourth new estimate 1408 can be adjusted to a value that is larger than (e.g., to the right of) that of the second estimate 1404 ( ⁇ ⁇ ⁇ ′ ).
  • the adjustment step size can be a constant (e.g., equal to the correlator spacing).
  • executing the second iterative process using the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the adjusted initial estimate ⁇ ⁇ ′ as the starting values includes: in accordance with a determination that (i) the updated estimate ⁇ ⁇ ⁇ , ⁇ (e.g., produced during a k th iteration of the second iterative process) is not within the corresponding first valid range or (ii) the updated estimate ⁇ ⁇ ⁇ is not within the corresponding second valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, the GNSS receiver repeats the steps of adjusting the initial estimate for the time delay ⁇ and executing the second iterative process (see, e.g., the path corresponding to steps 1226, 1228, and 1204 in Figure 12), and otherwise (e.g., in accordance with a determination that a next adjusted value for the time delay would not fall
  • the adjustment range (e.g., ⁇ ⁇ 1420 to ⁇ ⁇ 1422 in Figure 14) for the time delay ⁇ should cover all correlation sample times (e.g., from ⁇ ⁇ to ⁇ ⁇ ) with some margin. Once an adjustment moves beyond this range, it means that all possible initial estimates in this corresponding direction have been tried. If adjustments have covered the whole range without producing an updated estimate ⁇ ⁇ ⁇ , ⁇ within the corresponding first valid range and an updated estimate ⁇ ⁇ ⁇ within the corresponding second valid range, the GNSS receiver stops the second iterative process and outputs information that there is a no valid solution (e.g., see, e.g., Figure 12, steps 1230, 1231).
  • the GNSS receiver corrects the pseudorandom phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • the GNSS receiver after outputting the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ , ⁇ , and ⁇ , respectively (see, e.g., Figure 12, step 1218), and in accordance with the determination that the phase ⁇ ⁇ is not 0° or 180° within the second predefined margin (see, e.g., Figure 12, step 1238): the GNSS receiver starts with the initial estimate for ⁇ ⁇ ⁇ , ⁇ (e.g., see above discussion of step 1102, Figure 11) and obtains a solution for ⁇ ⁇ via a third iterative process.
  • k is a positive integer (e.g., starting from one (e.g., reset to 1 at step 1240 of the flowchart in Figure 12) representing a number of iterations of the third iterative process.
  • the one or more predetermined second parameters include the parameter ⁇ ⁇ .
  • the determination that the third iterative process has converged is made in accordance with a determination that the change in an estimated value of ⁇ ⁇ between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ ⁇ ).
  • the one or more predetermined second parameters include a post-fix mean squared error (MSE) for the quadrature component estimation error, defined by ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ for the k th iteration (See Equation 44).
  • MSE mean squared error
  • the determining that the third iterative process has converged is in accordance with a determination that a change in the post-fix MSE for the quadrature component estimation error between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e., ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ , ⁇ ), and/or if ⁇ ⁇ ⁇ ⁇ , ⁇ itself is smaller than (or does not exceed) a corresponding threshold.
  • a corresponding threshold i.e., ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ , ⁇
  • the k th iteration of the third iterative process further includes, in accordance with a determination that the change in value of the one or more predetermined second parameters between the k th iteration and the (k-1) th iteration does not satisfy a respective corresponding threshold, performing a next iteration of the third iterative process (see, e.g., Figure 12, step 1252).
  • the change in value of the one or more predetermined second parameters between the k th iteration and the (k-1) th iteration does not satisfy a respective corresponding threshold when or when ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ), or when ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ , ⁇ ), or when [00200]
  • the k th iteration of the third iterative process further includes: in accordance with the determination that the number of iterations ⁇ does not satisfy the sixth threshold value (e.g., k > Threshold5 or k ⁇ Threshold5 in Figure 12), terminating the third iterative process without producing a solution for ⁇ ⁇ (see, e.g., Figure 12, steps
  • the GNSS receiver corrects the pseudorange measurement for the respective satellite in accordance with ⁇ ⁇ and/or corrects the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • the GNSS receiver uses the initial estimates as starting values and obtains a solution for ⁇ ⁇ , ⁇ ⁇ , and ⁇ via a fourth iterative process (see, using workflow 1300 in Figure 13).
  • k is a positive integer (e.g., starting from one when the fourth iterative process begins) representing a number of iterations of the fourth iterative process.
  • a k th iteration of the fourth iterative process includes: obtaining n estimated quadrature component values ( ⁇ ⁇ ⁇ ( ⁇ ) values) for the Q(t) samples based on initial estimates and ⁇ ⁇ ⁇ from a previous (k-1) th iteration of the fourth iterative process (e.g., using equation 29) (see, e.g., Figure 13, step 1304), or for a first iteration of the fourth iterative process, based on previously determine initial estimates (e.g., as explained above with reference to step 1102, Figure 11); determining, for each estimated quadrature component value ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ of the ⁇ ⁇ ⁇ ( ⁇ ) values, a respective quadrature component estimation error ⁇ ⁇ , ⁇ between (i) the respective estimated quadrature component value ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ and (ii) the respective quadrature component sample Q(t i ) (e.g., using equation
  • equations 32 and 33 are solved separately, while estimates of ⁇ ⁇ , ⁇ ⁇ and ⁇ for equation 33 are updated based on equations 35, 36 and 37.
  • the two ⁇ ⁇ 3 matrix equations (equations 32 and 33) are combined into one ( 2 ⁇ ) ⁇ 4 matrix equation with a total of four unknowns (e.g., unknowns ⁇ ⁇ , ⁇ ⁇ , ⁇ ⁇ , and ⁇ ⁇ ), and the ( 2 ⁇ ) ⁇ 4 matrix equation is solved (e.g., instead of combining solutions).
  • the one or more predetermined second parameters include one or more of ⁇ ⁇ , ⁇ ⁇ , and ⁇ .
  • the determination that the fourth iterative process has converged is made in accordance with a determination that (i) the change in an estimated value of ⁇ ⁇ between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., and/or (ii) the change in an estimated value of ⁇ between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ), and/or (iii) the change in an estimated value of ⁇ between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇
  • MSE mean squared error
  • determining that the fourth iterative process has converged is in accordance with a determination that a change in the post-fix MSE for the quadrature component estimation error between the k th iteration and the (k-1) th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ iteration of the fourth iterative process further includes: in accordance with the determination that the change in value of the one or more predetermined parameters between the k th iteration and the (k-1) th iteration does not satisfy the respective corresponding threshold (e.g., the fourth iterative process has not converged, see Figure 13, Step 1328), performing a next iteration of the fourth iterative process ( Figure 13, Step 1330).
  • the change in value of the one or more predetermined parameters between the k th iteration and the (k-1) th iteration does not satisfy the respective corresponding threshold when ⁇ ⁇ ⁇ ⁇ ⁇ , ⁇ ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ), or ⁇ h ⁇ ⁇ ⁇ h ⁇ ⁇ ⁇ ⁇ , ⁇ ).
  • the k th iteration of the fourth iterative process further includes: in accordance with the determination that the number of iterations ⁇ does not satisfy the seventh threshold value (e.g., k ⁇ Threshold6 or k > Threshold6 in Figure 13, step 1310) (e.g., the fourth iterative process has diverged or failed to converge), terminating the fourth iterative process without producing a solution for ⁇ ⁇ (see, e.g., Figure 13, step 1326).
  • the seventh threshold value e.g., k ⁇ Threshold6 or k > Threshold6 in Figure 13, step 1310
  • the GNSS receiver updates the initial estimate for the time delay ⁇ from ⁇ ⁇ ⁇ to ⁇ ⁇ ⁇ ′ .
  • the GNSS receiver executes the fourth iterative process (i.e., executes the fourth iterative process an additional time) using the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the updated initial estimate ⁇ ⁇ ⁇ ′ as the starting values.
  • the GNSS receiver uses the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the adjusted initial estimate ⁇ ⁇ ⁇ ′ as the starting values.
  • the GNSS receiver uses the initial estimates ⁇ ⁇ ⁇ , ⁇ and ⁇ ⁇ ⁇ , ⁇ and the adjusted initial estimate ⁇ ⁇ ⁇ ′ as the starting values.
  • the GNSS receiver repeats the steps of adjusting the initial estimate for the time delay ⁇ and executing the fourth iterative process with the updated initial estimate for the time delay ⁇ .
  • the GNSS receiver in conjunction with outputting the updated estimates ⁇ ⁇ ⁇ , ⁇ , ⁇ ⁇ ⁇ , ⁇ , and ⁇ ⁇ ⁇ as the solution for ⁇ , ⁇ , and ⁇ , respectively, the GNSS receiver obtains Because Procedure 1 (see workflow 900, Figure 9) determines that ⁇ ⁇ is equal to ⁇ 90° (see steps 906, 908, 910, 912, Figure 12, and step 1302, Figure 13) using ⁇ ( ⁇ ⁇ ) measurements based on the same equation, equation 22, use to solve for ⁇ ⁇ , in Procedure 2 (see, e.g., workflow 1100 in Figure 11 and workflow 1300 in Figure 13), there is no need to repeat this process of solving for some unknowns from ⁇ ( ⁇ ⁇ ) measurements.
  • the solution ⁇ ⁇ to equation 22 obtained in Procedure 1 can be assigned not only to the initial estimate ⁇ ⁇ ⁇ , ⁇ but also to the final solution of ⁇ ⁇ .
  • the GNSS receiver computes the carrier phase multipath error ⁇ ⁇ based on the solution for ⁇ ⁇ and the solution for ⁇ ⁇ (e.g., using Equation 46).
  • the GNSS receiver corrects the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ⁇ ⁇ .
  • ⁇ ⁇ and ⁇ ⁇ are known.
  • first means “first,” “second,” etc.
  • these elements should not be limited by these terms. These terms are only used to distinguish one element from another.
  • a first contact could be termed a second contact, and, similarly, a second contact could be termed a first contact, without changing the meaning of the description, so long as all occurrences of the “first contact” are renamed consistently and all occurrences of the second contact are renamed consistently.
  • the first contact and the second contact are both contacts, but they are not the same contact.
  • the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the claims.
  • the term “if” may be construed to mean “when” or “upon” or “in response to determining” or “in accordance with a determination” or “in response to detecting,” that a stated condition precedent is true, depending on the context.
  • the phrase “if it is determined [that a stated condition precedent is true]” or “if [a stated condition precedent is true]” or “when [a stated condition precedent is true]” may be construed to mean “upon determining” or “in response to determining” or “in accordance with a determination” or “upon detecting” or “in response to detecting” that the stated condition precedent is true, depending on the context.

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Digital Transmission Methods That Use Modulated Carrier Waves (AREA)

Abstract

A Global Navigation Satellite System (GNSS) receiver receives a band-limited composite signal corresponding to a respective satellite in the GNSS. The composite signal is composed of a direct-path signal and a multipath signal. The receiver obtains, for a code chip edge transition of a pseudorandom code, correlation samples of the composite signal (I(t) samples and Q(t) samples) and compares them with a pre-determined filter characteristic (SR(t) samples) corresponding to a filter used to band limit the composite signal. Based on the comparison, the receiver determines a phase ϕ3rd and/or a response time error of the code chip edge transition due to the multipath signal, Δt. The receiver adjusts a pseudorange measurement for the respective satellite and/or adjusts a carrier phase measurement for the respective satellite in accordance the determined ϕ3rd and/or Δt. The receiver performs a navigation function using the adjusted pseudorange and/or carrier phase measurement for the respective satellite.

Description

PHASE MULTIPATH MITIGATION IN A NAVIGATION SATELLITE SIGNAL RECEIVER RELATED APPLICATIONS [0001] This application is a continuation of U.S. Patent Application No.18/373,913, filed September 27, 2023, and claims priority to U.S. Provisional Patent Application No. 63/472,574, filed June 12, 2023, which is hereby incorporated by reference in its entirety. TECHNICAL FIELD [0002] The disclosed embodiments relate generally to a global navigation satellite systems (GNSS) and more specifically, to systems and methods for mitigating multipath- induced errors in GNSS receiver positioning. BACKGROUND [0003] Global Navigation Satellite Systems (GNSS), such as Global Positioning System (GPS), the Global Orbiting Navigation Satellite System (GLONASS), the European Galileo and EGNOS, Chinese BeiDou System (BDS), the Wide Area Augmentation System (WAAS), and Japanese Quasi-Zenith Satellite System (QZSS), are used in many applications. In a GNSS system, each satellite transmits a signal that identifies the satellite and allows a receiver to determine the time at which the signal was sent. To do so, a GNSS satellite transmits a pseudorandom code (PN code) (also called a pseudorandom noise code (PRN code)). The PN code is, for example, a series of ones and zeroes that looks random, but in fact uniquely identifies the satellite. The GNSS receivers acquire and track the GNSS satellite signals. From its code and carrier tracking loops, a GNSS receiver can then generate range or range rate measurements, such as pseudorange measurements, carrier phase measurements and doppler shift measurements. In some instances, using other data, such as ephemeris data and almanac data, that is encoded in the received signal, “on top of” the pseudorandom code, using various modulation schemes (e.g., BPSK, BOC), the receiver can determine the time at which the signal was sent by the satellite. Using signals received from four or more satellites, the receiver can determine its position (e.g., on Earth). [0004] Multipath is one of the major error sources in GNSS. When a signal from a GNSS satellite to a GNSS receiver is reflected by surroundings (e.g., buildings, mountains, lakes, other objects, and/or the ground), this reflected signal, called a multipath, can also be received by the GNSS receiver. Therefore, a GNSS receiver can receive a composite of a line-of-sight direct-path signal and zero, one or more multiple signals from each satellite. Just like direct-path signals, multipath signals are also internally processed by GNSS receivers, thus introducing errors to code and carrier phase measurements and reducing GNSS positioning accuracy. [0005] To improve GNSS range measurement accuracy and positioning accuracy, many techniques have been invented to mitigate or eliminate code multipath errors and carrier phase multipath errors. However, improving the performance of these multipath migration techniques, for example, their practicable implementation, reliability, effectiveness, and accuracy, is always of importance and of interest in the GNSS industry. Accordingly, there is a need for a technique to mitigate the effect of multipath-induced errors. SUMMARY [0006] In accordance with some embodiments, a method of mitigating the effect of a multipath-induced error in a GNSS is performed at a respective GNSS signal receiver. The method includes receiving a band-limited composite signal corresponding to a respective satellite in the GNSS. The band-limited composite signal includes a band-limited direct-path signal and a band-limited multipath signal. The direct-path signal and the multipath signal are modulated with a pseudorandom (PN) code. The method includes obtaining, for a code chip edge transition of the PN code, n pairs of correlation samples of the composite signal during a corresponding n time instances, each having a different time offset from the code chip edge transition. The code chip edge transition has a predetermined filter step response function SR(t). Each respective pair of correlation samples of the n pairs of correlation samples consists of a respective in-phase component sample I(ti) and a respective quadrature component sample Q(ti), wherein i is a positive integer from one to n. The n pairs of correlation samples consist of n in-phase component samples of the composite signal (I(t) samples) and n quadrature component samples of the composite signal (Q(t) samples). The method includes obtaining, for each respective pair of correlation samples of the n pairs of correlation samples, a corresponding sample of the filter step response function SR(ti) corresponding in time with the respective pair of correlation samples, thereby obtaining n samples of the filter step response function (predetermined SR(t) samples) during the corresponding n time instances. The method includes, in accordance with a determination that a computed similarity between the I(t) samples and the SR(t) samples satisfies a first threshold value, determining that a phase ^^ଷ^ௗ is ±90° within a first predefined margin, wherein the phase ^^ଷ^ௗ is equal to 180° minus a sum of (i) a carrier phase multipath error ^^ఌ and (ii) a phase difference ^^^ between the direct-path signal and the multipath signal. The method includes, in accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples does not satisfy the first threshold value: solving a first set of matrix equations, thereby obtaining a solution for tan ^^ଷ^ௗ and for ∆ ^^ tan ^^ଷ^ௗ, wherein ∆ ^^ is a response time error of the code chip edge transition due to the multipath signal. The method includes, in accordance with a determination that | ^^ ^^ ^^ ^^ଷ^ௗ| satisfies a second threshold value, determining that the phase ^^ଷ^ௗ is 0° or 180° within a second predefined margin. The method includes, in accordance with a determination that | ^^ ^^ ^^ ^^ଷ^ௗ| does not satisfy the second threshold value, determining ^^ଷ^ௗ in accordance with the solution for tan ^^ଷ^ௗ. The method includes adjusting a pseudorange measurement for the respective satellite in accordance with the determined ∆ ^^; and/or adjusting a carrier phase measurement for the respective satellite in accordance with a parameter corresponding to the determined ^^ଷ^ௗ. The method includes performing a navigation function using the adjusted pseudorange and/or the adjusted carrier phase measurement for the respective satellite. [0007] In some embodiments, the method further includes determining a sign of the carrier phase multipath error ^^ఌ. The method incudes, in accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples satisfies the first threshold value: determining that ^^ଷ^ௗ is +90° within the first predefined margin when the sign of the carrier phase multipath error ^^ఌ is positive; and determining that ^^ଷ^ௗ is -90° within the first predefined margin when the sign of the carrier phase multipath error ^^ఌ is negative. [0008] In some embodiments, the method further includes, in accordance with a determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: determining a position of a peak of the I(t) samples relative to a position of a peak of the predetermined SR(t) samples. The method includes, in accordance with a determination that the position of the peak of the I(t) samples occurs later in time than the position of the peak of the SR(t) samples, determining that ^^ଷ^ௗ is 0° within the second predefined margin; and in accordance with a determination that the position of the peak of the I(t) samples occurs earlier in time than the position of the peak of the predetermined SR(t) samples, determining that ^^ଷ^ௗ is 180° within the second predefined margin. [0009] In some embodiments, the computed similarity between the I(t) samples and the SR(t) samples is obtained using a mean squared error (MSE), defined by ^^ ^^ ^^^^ ൌ ^ ^ ∑^ ^ୀ^ ^ ^^ ( ^^^ ) − ^^ூ ^^ ^^( ^^^) ^ଶ , wherein ^^ூ ൌ ^^ௗ cos ^^ఌ, ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error. [0010] In some embodiments, solving the first set of matrix equations includes obtaining (i) an initial estimate ^^ ^ ଷ^ௗ,^ for the phase ^^ଷ^ௗ, (ii) an initial estimate ∆ ^ ^^^ for the response time error ∆ ^^, and (iii) an initial estimate ^ ^ ^ௗା^ for a magnitude ^^ௗା^ of the composite signal; and using the initial estimates ^^ ^ ଷ^ௗ,^ , ∆ ^ ^^^, and ^ ^ ^ௗା^ to solve the first set of matrix equations via a Least Squares fitting process to obtain the solution for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ. [0011] In some embodiments, the method further includes, in accordance with the determination that | ^^ ^^ ^^ ^^ଷ^ௗ | does not satisfy the second threshold value, obtaining a solution for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ via a first iterative process, each iteration of the first iterative process including: obtaining an updated estimate tan ^^ ^ ଷ^ௗ,^ for the phase ^^ଷ^ௗ and an updated estimate ∆^ ^^^ for the response time error ∆ ^^. The method includes, in accordance with a determination that a number of iterations for which the first iterative process has looped satisfies a third threshold value: outputting information that no valid solution for ^^ଷ^ௗ has been determined; and terminating the first iterative process without producing a solution for ^^ଷ^ௗ. The method further includes, in accordance with a determination that a number of iterations for which the first iterative process has looped does not satisfy the third threshold value: in accordance with a determination that a magnitude based on the updated estimate ∆ ^ ^^^ or based on the updated estimate tan ^^ ^ ଷ^ௗ,^ satisfies a fourth threshold value: determining that the Least Squares fitting process has converged; outputting the updated estimates tan ^^^ଷ^ௗ,^ and the updated estimate ∆^ ^^^ as the solutions for tan ^^ଷ^ௗ and ∆ ^^, respectively; calculating the phase ^^ଷ^ௗ from the updated estimate tan ^^^ ଷ^ௗ,^ ; and terminating the first iterative process. The method also includes, in accordance with a determination that (i) the magnitude based on the updated estimate ∆ ^ ^^^ or based on the updated estimate tan ^^ ^ ଷ^ௗ,^ does not satisfy the fourth threshold value, performing a next iteration of the first iterative process. [0012] In some embodiments, the method further includes, in accordance with a determination that a valid solution for the phase ^^ଷ^ௗ has been determined: determining (i) an initial estimate ^ ^ ^ூ,^ for ^^ூ, wherein ^^ூ ൌ ^^ௗ cos ^^ఌ, ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error; (ii) an initial estimate ^^^ ொ,^ for a magnitude ^^ொ, wherein ^^ொ ൌ ^^ௗ sin ^^ఌ; (iii) an initial estimate ^ ^ ^^,^ for a magnitude ^^^ of the multipath signal, and (iv) an initial estimate ^ ^ ^^ for a time delay ^^ of the multipath signal relative to the direct-path signal, or determine a subset of those initial estimates. The method includes using at least a subset of (i) the initial estimate ^ ^ ^ூ,^ , (ii) the initial estimate ^^ ^ ொ,^ , (iii) the initial estimate ^ ^ ^^,^, and (iv) the initial estimate ^ ^^^ as starting values, performing one or more iterative processes to obtain solutions for at least a subset of ^^ூ, ^^ொ, ^^^, and ^^; computing the carrier phase multipath error ^^ఌ in accordance with the obtained solutions for at least the subset of ^^ூ, ^^ொ, ^^^, and ^^; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [0013] In some embodiments, the method further includes, in accordance with a determination that a valid solution for the phase ^^ଷ^ௗ has been determined: determining (i) an initial estimate ^^^ூ,^ for ^^ூ, wherein ^^ூ ൌ ^^ௗ cos ^^ఌ, ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error; (ii) an initial estimate ^^^ ொ,^ for a magnitude ^^ொ, wherein ^^ொ ൌ ^^ௗ sin ^^ఌ; (iii) an initial estimate ^^^^,^ for a magnitude ^^^ of the multipath signal, and (iv) an initial estimate ^ ^ ^^ for a time delay ^^ of the multipath signal relative to the direct-path signal. [0014] In some embodiments, in accordance with a determination that the code chip edge transition is a code chip edge-up transition: the initial estimate ^ ^ ^ூ,^ corresponds to a magnitude of a last in-phase component sample I(tn) of the I(t) samples; and the initial estimate ^^ ^ ொ,^ corresponds to a magnitude of a last quadrature component Q(tn) of the Q(t) samples. [0015] In some embodiments, in accordance with the code chip edge transition is a code chip edge-down transition: the initial estimate ^^^ூ,^ corresponds to a magnitude of a first in-phase component sample I(t1) of the I(t) samples; and the initial estimate ^^^ ொ,^ corresponds to a magnitude of a first quadrature component Q(t1) of the Q(t) samples. [0016] In some embodiments, the method includes, in accordance with the determination that the phase ^^ଷ^ௗ corresponds to ±90° within the first predefined margin:
Figure imgf000007_0001
assigning the obtained solution as the initial estimate ^ ^ ^ ^^,0. [0017] In some embodiments, the method includes, in accordance with the determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: assigning the initial estimate ^^ ^ ொ,^ as zero. [0018] In some embodiments, the initial estimate ^ ^ ^^,^ is 0.5. [0019] In some embodiments, the method further includes, in accordance with the determination that the phase ^^ଷ^ௗ does not correspond to ±90° within the first predefined margin: using the initial estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^^^ as starting values to obtain a solution for ^^ூ, ^^^, and ^^ via a second iterative process, wherein k is a positive integer representing a number of iterations of the second iterative process. A kth iteration of the second iterative process includes: obtaining n estimated in-phase component values ( ^^^( ^^^ି^) values) for the I(t) samples based on estimates ^^^ , , an ^ th ூ,^ି^
Figure imgf000008_0001
d ^^^ି^ from a previous (k-1) iteration of the second iterative process; determining, for each estimated in-phase component value ^ ^ ^൫ ^^^,^ି^൯ of the ^ ^ ^( ^^^ି^) values, a corresponding in-phase component estimation error
Figure imgf000008_0002
between (i) the respective estimated in-phase component value ^^^൫ ^^^,^ି^൯ and (ii) the respective in-phase component sample I(ti) of the I(t) samples, thereby obtaining n in-phase component estimation errors ( ^^ூ ( ^^^ ) estimation errors); solving a second set of matrix equations based on the ^^ூ ( ^^^ ) estimation errors, thereby obtaining updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ for the k th iteration. The method includes, in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is within a corresponding first valid range and (ii) the updated estimate ^ ^ ^^ is within a corresponding second valid range: in accordance with a determination that the number of iterations ^^ satisfies a fifth threshold value: in accordance with a determination that a change in value of one or more predetermined first parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the second iterative process has converged; outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively; and terminating the second iterative process. The method includes, in accordance with the determination that the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold, performing a next iteration of the second iterative process. The method includes, in accordance with the determination that (i) the number of iterations ^^ does not satisfy the fifth threshold value, terminating the second iterative process without producing a solution for ^^ூ, ^^^, and ^^. [0020] In some embodiments, the method includes, in accordance with a determination that (i) the updated estimate ^^^^,^ is not within the corresponding first valid range or (ii) the updated estimate ^ ^ ^^ is not within the corresponding second valid range: adjusting the initial estimate for the time delay ^^ from ^ ^ ^^ to ^ ^ ^^′ ; and executing the second iterative process using the initial estimates ^ ^ ^ூ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^^^′ as the starting values. [0021] In some embodiments, the method includes using the initial estimates ^^^ூ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^ ^^′ as the starting values for the second iterative process (e.g., restarted from the beginning, with k=1). The method includes, in accordance with a determination that (i) the updated estimate ^^^^,^ is not within the corresponding first valid range or (ii) the updated estimate ^ ^ ^^ is not within the corresponding second valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, repeating the steps of adjusting the initial estimate for the time delay ^^ and executing the second iterative process; and otherwise, outputting information that no valid solution for ^^ூ, ^^^, and ^^ has been determined. [0022] In some embodiments, the method includes, after outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively, and in accordance with the determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: assigning a value of zero as the solution for ^^ொ. [0023] In some embodiments, the method includes determining the carrier phase multipath error ^^ఌ based on the solutions for ^^ூ and ^^ொ; and correcting a pseudorandom phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [0024] In some embodiments, the method includes, after outputting (e.g., at an end of the second iterative process) the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively, and in accordance with the determination that the phase ^^ଷ^ௗ is not 0° or 180° within the second predefined margin: starting with the initial estimate for ^^^ ொ,^ , obtaining a solution for ^^ொ via a third iterative process, wherein k is a positive integer representing a number of iterations of the third iterative process (e.g., k is set or reset to 1 at beginning of the third iterative process). A kth iteration of the third iterative process includes: obtaining n estimated quadrature component values ( ^ ^ ^ ( ^^^ି^ ) values) for the Q(t) samples based on an initial estimate
Figure imgf000010_0001
from a previous (k-1) th iteration of the third iterative process; determining, for each estimated quadrature component value ^^^൫
Figure imgf000010_0002
of the ^^^( ^^^ି^) values, a corresponding quadrature component estimation error ^^ொ൫ ^^^,^൯ between (i) the respective estimated quadrature component value ^ ^ ^൫ ^^^,^ି^൯ and (ii) the respective quadrature component sample Q(ti), thereby obtaining n quadrature component estimation errors estimation errors); solving a third set of matrix equations based on the
Figure imgf000010_0003
estimation errors, thereby obtaining an updated estimate ^^ ^ ொ,^ for the k th iteration. In accordance with a determination that the number of iterations ^^ satisfies a sixth threshold value, the method includes, in accordance with a determination that a change in value of one or more predetermined second parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the third iterative process has converged; outputting the updated estimate
Figure imgf000010_0004
the solution for ^^ொ; and terminating the third iterative process. The method includes, in accordance with a determination that the change in value of the one or more predetermined second parameters between the kth iteration and the (k-1)th iteration does not satisfy a respective corresponding threshold, performing a next iteration of the third iterative process. The method includes, in accordance with the determination that (i) the number of iterations ^^ does not satisfy the sixth threshold value, terminating the third iterative process without producing a solution for ^^ொ. [0025] In some embodiments, the method includes computing the carrier phase multipath error ^^ఌ based on the solutions for ^^ூ and ^^ொ; correcting the pseudorange measurement for the respective satellite in accordance with ∆ ^^; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [0026] In some embodiments, the method includes, in accordance with the determination that the phase ^^ଷ^ௗ corresponds to ±90° within the first predefined margin: using the initial estimates
Figure imgf000010_0005
as starting values, obtaining a solution for ^^ ^^, ^^ ^^, and ^^ via a fourth iterative process, wherein k is a positive integer representing a number of iterations of the fourth iterative process. A kth iteration of the fourth iterative process includes: obtaining n estimated quadrature component values ( ^ ^ ^ ( ^^^ି^ ) values) for the Q(t) samples based on initial estimates ^^ ^ ொ,^ି^ , ^^ ^ ^,^ି^ , and ^ ^ ^^ି^ from a previous (k-1) th iteration of the fourth iterative process; determining, for each estimated quadrature component value ^ ^ ^൫ ^^^,^ି^൯ of the ^ ^ ^ ( ^^^ି^ ) values, a respective quadrature component estimation error ^^ொ൫ ^^^,^൯ between (i) the respective estimated quadrature component value
Figure imgf000011_0001
and (ii) the respective quadrature component sample Q(ti), thereby obtaining n quadrature component estimation errors ( ^^ொ ( ^^^ ) estimation errors); and solving a fourth set of matrix equations based on the
Figure imgf000011_0002
estimation errors, thereby obtaining updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^^^^ for the kth iteration. The method includes, in accordance with a determination that (i) the updated estimate ^^^^,^ is within a corresponding third valid range and (ii) the updated estimate ^^^^ is within a corresponding fourth valid range: in accordance with a determination that the number of iterations ^^ satisfies a seventh threshold value: in accordance with a determination that a change in value of one or more predetermined second parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the fourth iterative process has converged; outputting the updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ ^^, ^^ ^^, and ^^, respectively; and terminating the fourth iterative process. The method includes, in accordance with the determination that the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold, performing a next iteration of the fourth iterative process. The method includes, in accordance with the determination that (i) the number of iterations ^^ does not satisfy the seventh threshold value, terminating the fourth iterative process without producing a solution for ^^ ^^. [0027] In some embodiments, the method includes, in accordance with a determination that (i) the updated estimate ^^^^,^ is not within a corresponding third valid range or (ii) the updated estimate ^^^^ is not within a corresponding fourth valid range: updating the initial estimate for the time delay ^^ from ^ ^ ^^ to ^ ^ ^^′ ; and executing the fourth iterative process using the initial estimates ^^ ^ ொ,^ and ^ ^ ^^,^ and the updated initial estimate ^ ^ ^^′ as the starting values. [0028] In some embodiments, the method includes using the initial estimates ^^^ ொ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^ ^^′ as the starting values. The method includes, in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is not within the corresponding third valid range or (ii) the updated estimate ^^^^ is not within the corresponding fourth valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, repeating the steps of updating the initial estimate for the time delay ^^ and executing the fourth iterative process; and otherwise, outputting information that no valid solution for ^^ ^^, ^^ ^^, and ^^ has been determined. [0029] In some embodiments, the method includes, in conjunction with outputting the updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ ^^, ^^ ^^, and ^^, respectively: obtaining a solution for ^^ூ by solving a matrix equation
Figure imgf000012_0001
[0030] In some embodiments, the method includes computing the carrier phase multipath error ^^ఌ based on the solution for ^^ூ and the solution for ^^^; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [0031] In accordance with some embodiments, a navigation module, or a global navigation satellite system (GNSS) receiver, or a device for mitigating the effect of a multipath-induced error in a global navigation satellite system (GNSS) includes one or more processors and memory. The memory stores instructions that, when executed by the one or more processors, cause the navigation module, or GNSS receiver, or device, to perform any of the methods disclosed herein. [0032] In accordance with some embodiments, a computer-readable storage medium is provided that stores instructions which, when executed by a navigation module, or a global navigation satellite system (GNSS) receiver, or a device, that includes one or more processors, cause the navigation module, or the GNSS receiver, or the device, to perform the any of the methods described herein. BRIEF DESCRIPTION OF THE DRAWINGS [0033] For a better understanding of the various described embodiments, reference should be made to the Description of Embodiments below, in conjunction with the following drawings in which like reference numerals refer to corresponding parts throughout the figures. [0034] Figure 1A is a block diagram illustrating a navigation system, according to some embodiments. [0035] Figure 1B illustrates a block diagram of a GNSS receiver in accordance with some embodiments. [0036] Figure 1C is a block diagram of a computer system, such as a computer system that is part of a moveable object’s navigation system, according to some embodiments. [0037] Figure 2 illustrates a block diagram of a baseband DSP channel in accordance with some embodiments. [0038] Figure 3 illustrates an example where both a direct-path signal and a multipath signal from a satellite are received by an antenna of a GNSS receiver, in accordance with some embodiments. [0039] Figure 4A illustrates a phasor diagram of the in-phase and quadrature components of a direct-path signal, in accordance with some embodiments. Figure 4B illustrates a phasor diagram of the in-phase and quadrature components of a composite signal, which is composed of a direct-path signal and a multipath signal, in accordance with some embodiments. [0040] Figure 5 illustrates the magnitude of a filter step response function ^^ ^^( ^^) versus time, in accordance with some embodiments. [0041] Figure 6 illustrates a correlator arrangement in which 16 pairs of samples are arranged along a code chip transition, and a delayed version of the filter step response function ^^ ^^( ^^ − ^^), in accordance with some embodiments. [0042] Figures 7A and 7B are phasor diagrams illustrating multipath scenarios where angle ^^ଷ^ௗ is equal to +90° and -90°, respectively, in accordance with some embodiments. [0043] Figures 8A and 8B are phasor diagrams illustrating multipath scenarios where angle ^^ଷ^ௗ is equal to 0° and 180°, respectively, in accordance with some embodiments. [0044] Figure 9 is a workflow for Procedure 1 of the phase multipath mitigation algorithm, in accordance with some embodiments [0045] Figures 10A and 10B show an example of ^^( ^^^) measurement samples when angle ^^ଷ^ௗ is equal to 0° and 180°, respectively, in accordance with some embodiments. [0046] Figure 11 illustrates a workflow 1100 for Procedure 2 of the phase multipath mitigation algorithm, in accordance with some embodiments. [0047] Figure 12 illustrates a workflow for solving the unknown parameters, ^^ூ, ^^ொ, ^^^, and ^^ when ^^ଷ^ௗ is not ±90°, in accordance with some embodiments. [0048] Figure 13 illustrates a workflow for solving for the unknown parameters ^^ூ, ^^ொ, ^^^, and ^^ when ^^ଷ^ௗ is equal to ±90°, in accordance with some embodiments. [0049] Figure 14 illustrates a series of different initial estimates ^^^^ that can be tried sequentially, in accordance with some embodiments. [0050] Figure 15A illustrates a phasor diagram for a very weak multipath scenario with ^^ଷ^ௗ ൌ 30°. Figure 15B illustrates a correct solution to angle ^^ఌ even when this multipath scenario is detected as a case of ^^ଷ^ௗ ൌ ^90°. [0051] Figures 16A and 16B illustrate a flowchart diagram for a method of mitigating the effect of a multipath-induced error in a GNSS, in accordance with some embodiments. DESCRIPTION OF EMBODIMENTS [0052] Figure 1A is a block diagram illustrating a navigation system 100, according to some embodiments. Navigation system 100 enables a moveable object 110 (e.g., a phone, specialized GNSS receiver, a boat, a truck or other vehicle, a farming appliance, mining appliance, drilling system, etc.) to determine, at any point of time, its current position 112 with respect to a global coordinate system (e.g., a coordinate system for Earth 114). Moveable object 110 is equipped with a satellite receiver (navigation signal receiver 120), typically including or coupled to one or more satellite antennas 140, to receive satellite navigation signals from at least four satellites 115 that are orbiting Earth. The satellite navigation signals received by moveable object 110 are typically global navigation satellite system (GNSS) signals. [0053] Moveable object 110 optionally receives satellite orbit correction information and satellite clock correction information (sometimes collectively called “correction information”) for the plurality of satellites. The correction information is typically broadcast by and received from one or more satellites 118 distinct from the GNSS satellites 115, using an antenna 142 and signal receiver 152 (see Figure 1C) distinct from antenna 140 and navigation signal receiver 120 used to receive the satellite navigation signals. However, in some embodiments, the same antenna and receiver are used to receive both satellite navigation signals and correction information [0054] Moveable object 110 determines a position of moveable object 110, using the received satellite navigation signals and, optionally, the received satellite orbit correction information and satellite clock correction information for the plurality of satellites. In some embodiments, received satellite navigation signals are processed by navigation signal receiver 120, including analog signal processing circuitry 122 and a digital signal processor 124, optionally taking into account the correction information, to determine code measurements and phase measurements for signals received from four or more satellites 115. Embedded computer system 130 determines the position of moveable object 110 based on those measurements. [0055] Figure 1B illustrates a block diagram of an example GNSS receiver 170 (e.g., navigation signal receiver 120) in accordance with some embodiments. The GNSS receiver 170 receives signals from one or more GNSS satellites (e.g., GNSS satellite 171), processes the signals, and delivers Position, Velocity and Time (PVT) results 180. In some embodiments, the GNSS receiver 170 includes an antenna 172, a radio frequency (RF) front end module 174, a baseband digital signal processing (DSP) module 176, and a positioning and navigation module 178. The antenna 172 is the first stage of the GNSS receiver 170 and receives RF signals from all visible GNSS satellites. In some embodiments, the antenna 172 is equipped with a built-in amplifier. In some embodiments, the GNSS receiver 170 is equipped with multiple antennas or an antenna array. The received RF signals are filtered, amplified, down-converted, and digitized by the RF front end module 174. The RF front end module 174 typically includes one or more stages of down-conversions and filters, to convert the RF signals to different intermediate frequency (IF) bands with designated bandwidths while excluding out-of-band interferences. The final filter in high-precision GNSS receivers typically has a wide bandwidth (e.g., 30 MHz or more). One important characteristic of the final filter of the RF front end module 174 is its step response, which will be described later. Analog-to-digital (A/D) converters in the RF front end module convert received analog signals, e.g., from the final filter, into digital Intermediate Frequency (IF) signals. In some embodiments, the RF front end electronics 174 are integrated into a RF Application-Specific Integrated Circuit (ASIC) chip. [0056] The RF front end module 174 outputs digital IF signals, which are processed by the baseband DSP module 176. By processing the digital IF signals, the baseband DSP module 176 acquires and tracks each individual visible satellite (e.g., each satellite for which there is a direct line of sight to the receiver 170), generates range measurements, including code phase measurements (and therefore pseudorange measurements), carrier phase measurements and doppler measurements, and decodes navigation data bits encoded in the signal(s) received from each satellite. According to some embodiments of the present disclosure, phase multipath mitigation is performed by the baseband DSP module 176. In some embodiments, the baseband DSP module 176 comprises both hardware and firmware, but the division of these two portions is often not clearly defined and depends on a variety of design considerations. In some embodiments, the baseband DSP electronics are integrated into a baseband ASIC chip. In some embodiments, the RF front-end electronics and the baseband DSP electronics are integrated into a single-die System-On-Chip (SOC) chip. [0057] The baseband DSP module 176 outputs range measurements and navigation data that are used by the positioning and navigation module 178 to calculate PVT results 180. In some embodiments, the PVT results 180 are used by user applications. In some embodiments, the positioning and navigation module 178 is implemented in software executed by a microprocessor, while the firmware portion of the baseband DSP module 176 is implemented in either the same microprocessor as the microprocessor of the positioning and navigation module 178, or in a separate microprocessor. In some embodiments, the GNSS receiver 170 includes more than one microprocessor. [0058] Figure 1C is a block diagram of computer system 130 in, and used by, a moveable object (e.g., moveable object 110, Figure 1A) to determine the position of the moveable object, according to some embodiments. [0059] Computer system 130 typically includes one or more processors (sometimes called CPUs, or hardware processors) 202 for executing programs or instructions; memory 210; one or more communications interfaces 206; and one or more communication buses 205 for interconnecting these components. Computer system 130 optionally includes a user interface 209 comprising a display device 211 and one or more input devices 213 (e.g., one or more of a keyboard, mouse, touch screen, keypad, etc.) coupled to other components of computer system 130 by the one or more communication buses 205. Navigation signal receiver 120 (or GNSS receiver 170), and supplemental receiver(s) 152, if provided, are also coupled to other components of computer system 130 by the one or more communication buses 205. The one or more communication buses 205 may include circuitry (sometimes called a chipset) that interconnects and controls communications between system components. [0060] Communication interface 206 (e.g., a receiver or transceiver) is used by computer system 130, and more generally moveable object 110, to convey information to external systems (e.g., external system(s) 160, Figure 1A), and to receive communications from external systems. [0061] Memory 210 includes high-speed random access memory, such as DRAM, SRAM, DDR RAM or other random access solid state memory devices; and may include non-volatile memory, such as one or more magnetic disk storage devices, optical disk storage devices, flash memory devices, or other non-volatile solid state storage devices. Memory 210 optionally includes one or more storage devices remotely located from the CPU(s) 202. Memory 210, or alternately the non-volatile memory device(s) within memory 210, comprises a computer readable storage medium (e.g., a non-transitory computer readable storage medium). In some embodiments, memory 210 or the computer readable storage medium of memory 210 stores the following programs, modules and data structures, or a subset thereof: • an operating system 212 that includes procedures for handling various basic system services and for performing hardware dependent tasks; • a communications module 214 that operates in conjunction with communication interface 206 (e.g., a receiver and/or transceiver) to handle communications between moveable object 110 and external systems 160 (Figure 1); the connection between computer system 130 and external systems 160 may include a communication network 162, such as the internet or a public or proprietary wireless network; • optionally, a user interface module 216 for receiving information from one or more input device 213 of user interface 209, and to convey information to a user of moveable object 110 via one or more display or output devices 211; • a navigation module 218 for determining a position of the moveable object and performing one or more navigation functions (e.g., routing of the mobile object; displaying on a map one or more suggested routes for moving the mobile object from a current location to a specified or target location; and/or providing information regarding locations nearby the mobile object or along a proposed route of the mobile object); • an acquisition engine 220 for acquiring satellite signals from GNSS satellites, including determining a tracking frequency and a code shift (sometimes called a code offset) for each of several respective GNSS satellites, or for managing operation of an acquisition engine that is implemented, at least in part, in digital signal processor 124; and • a tracking module 222, sometimes called the satellite signal tracking module, which tracks GNSS satellite signals using acquisition information (e.g., the tracking frequency and code shift) handed over from the acquisition module 220, or for managing operation of a tracking module that is implemented, at least in part, in digital signal processor 124. For example, in some embodiments, tracking module 222 samples a GNSS satellite signal using a tracking frequency and code shift initially determined by the acquisition module 220, and then adjusted by the tracking module itself as the tracking module determines code phase corrections for respective channels of the received satellite navigation signals. [0062] Operating system 212 and each of the above identified modules and applications correspond to a set of instructions for performing a function described above. The set of instructions can be executed by the one or more processors 202 of computer system 130. The above identified modules, applications or programs (i.e., sets of instructions) need not be implemented as separate software programs, procedures or modules, and thus various subsets of these modules may be combined or otherwise re-arranged in various embodiments. In some embodiments, memory 210 stores a subset of the modules and data structures identified above. Furthermore, memory 210 optionally stores additional modules and data structures not described above. [0063] Figure 1C is intended more as a functional description of the various features which may be present in a computer system 130 of a moveable object 110 than as a structural schematic of the embodiments described herein. In practice, and as recognized by those of ordinary skill in the art, items shown separately could be combined and some items could be separated. For example, some items shown separately in Figure 1C could be combined into a single module or component, and single items could be implemented using two or more modules or components. The actual number of modules and components, and how features are allocated among them will vary from one implementation to another. [0064] In addition, in some embodiments, some or all of the above-described functions may be implemented with hardware circuits (e.g., which may comprise graphics processors for efficiently performing discrete Fourier transforms (DFTs), field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), a “system on a chip” that includes processors and memory, or the like). To that end, in some embodiments, CPUs 202 include specialized hardware for performing these and other tasks. In some embodiments, some of the operations described above with respect to computer system 130 are performed by navigation signal receiver 120 and/or GNSS receiver 170 rather than computer system 130. [0065] In some embodiments, the baseband DSP module 176 processes digital IF signals channel by channel (e.g., one channel per satellite being tracked). Figure 2 illustrates a block diagram of a baseband DSP channel 250 in accordance with some embodiments. Each channel may include hardware and software components configured to acquire and track one code signal of one satellite. In some multi-constellation, multi-frequency receivers, there can be tens or hundreds of nearly identical channels running and sharing the baseband DSP hardware resources, for example using time multiplexing. [0066] The baseband DSP module 176 receives (e.g., fetches) digital in-phase (I) IF samples 252 and quadrature (Q) IF samples 254 from the RF front-end module 174 of the GNSS receiver 170. The digital I IF samples 252 and Q IF samples 254 are first mixed with the replicas of the carrier signal in carrier removal mixers 256, so the carrier frequency of the received signals is stripped off, and the digital IF signals are down-converted to baseband I and Q samples, respectively. To completely strip off the carrier frequency, including doppler shift of the carrier frequency, the receiver generated carrier signal is synchronized with the frequency of the received IF signals. Maintaining carrier frequency synchronization is a primary goal of the carrier tracking loop in each channel of the receiver. Furthermore, when the carrier tracking loop of a DSP channel 250 is configured to align the phase of the receiver generated carrier signal with the carrier phase of the received IF signals, the carrier tracking loop is called a phase lock loop or phase-locked loop (PLL). [0067] The I and Q signals output by the carrier removal mixers 256 are then correlated with n replicas of code (e.g., the pseudorandom code corresponding to the satellite for which IF signals are being processed, provided by code phase management unit 272) in n pairs of correlators 258 (e.g., correlator pair 258-1, correlator pair 258-2, … and correlator pair 258-n). Controlled by code phase management unit 272, each of the n replicas of the code has a different phase delay. Typically, the code phase of one of the n replicas is better aligned with the code phase of the received signals than any of the other replicas, and the correlation of that replica with the received signals has a peak value relative to the other correlations. The strongest correlation is identified and used by the code tracking loop’s discriminator and filter 264, and the carrier tracking loop’s discriminator and filter 266, to adjust the code tracking loop and carrier tracking loop, and to remove the code (sometimes called a spreading spectrum code) from the received I and Q samples. After the received signals are de-spread by the correlators 258, the resulting signals represent the navigation data bits or symbols that are encoded (e.g., embedded) in the received satellite signal. Controlled by code phase management unit 272, each of the n replicas of the code has a different phase delay. Typically, the code phase of one of the n replicas is aligned with the code phase of the received signals, therefore their correlation can achieve a peak value. The goal of Code Tracking Loop in the receiver is to maintain this code phase synchronization or code phase lock. [0068] All 2n correlation results are coherently integrated for a pre-detection integration time in integrate and dump circuits 260, and then dumped for further signal detection tracking and/or multipath mitigation (e.g., via phase multipath mitigation unit 262), etc. The length of the pre-detection integration (PDI) time depends on the received signal’s navigation data bit rate or symbol rate and other design considerations. A typical value of the PDI is 1 ms. The n pairs of coherent integration I and Q samples are denoted as ^^( ^^^) 259-1I and ^^( ^^^) 259-1Q, ^^( ^^ଶ) 259-2I and ^^( ^^ଶ) 259-2Q, …, and ^^( ^^^) 259-nI and ^^( ^^^) 259-nQ. The coherent integration samples can be further coherently or non-coherently integrated for a longer duration in a receiver to achieve a higher signal to noise ratio (SNR). [0069] In some embodiments, a conventional design generates three replicas of code (n = 3), called Early, Prompt and Late, indicating their code phases relative to the code phase of the received signals, and thus this design obtains three samples of the signal’s auto- correlation function (ACF). Based on the Early, Prompt and Late correlation results, the code tracking loop’s discriminator and filter 264 can identify the code phase difference between the received signals and the Prompt replica of the code, filter it, and use it to adjust and drive code numerically controlled oscillator (NCO) 268. The output of the code NCO 268 is used by code generator or loader 270 to control the frequency and phase of the Prompt code replica such that it is kept aligned with those of the received signals. Based on the phase of the Prompt code, the code tracking loop’s discriminator and filter 264 can generate code phase measurements of the received signals. Using time information obtained by combining information obtained from the satellite signals of four or more satellites, the code phase measurements are used to generate pseudorange measurements, which are one type of key satellite range measurement for GNSS positioning (e.g., determining the position of a mobile object, such as mobile object 110, Figure 1A). [0070] Meanwhile, the carrier tracking loop’s discriminator and filter 266 typically utilizes only the pair of Prompt correlations, discriminates the frequency difference and/or phase difference between the received signals and the receiver generated carrier signal, filters the difference(s), and drives carrier NCO 274 so as to maintain synchronization between the received signals (e.g., from the satellite being tracked by the baseband DSP channel 250) and the carrier signal replica produced by carrier NCO 274 and carrier phase summation unit 276. In addition to carrier phase measurements, which is another type of key satellite range measurement for GNSS positioning, the carrier tracking loop’s discriminator and filter 266 can also output raw navigation data bits or symbols embedded in the satellite signal being processed by the baseband DSP channel 250. [0071] Figure 3 illustrates an example where both a direct-path signal 302 and a multipath signal 304 from a satellite 170 are received by an antenna 172 of a GNSS receiver. GNSS signals from satellites – whether line-of-sight (LOS) direct-path signals or multipath signals – that enter a receiver’s antenna will all be processed by the receiver. When a direct- path signal is reflected by surroundings, such as buildings 306, or mountains, or lakes, and/or other objects, it is called a multipath (or multipath signal). [0072] In some embodiments, a composite signal ^^( ^^) of a direct-path signal and its multipath signal can be mathematically modeled as: ^^( ^^) = ^^ௗ sin( ^^ ^^) + ^^^ sin( ^^( ^^ − ^^) + ^^) (1) where t represents the time that the composite signal is sampled or measured, ^^ௗ is an amplitude of the direct-path signal 302, ^^^ is an amplitude of the multipath signal 304, ^^ is a time delay of the multipath signal 304 relative to the direct-path signal 302 ( ^^ is also referred to as “multipath time delay”), and ^^ is a phase delay of the multipath signal 304 relative to the direct-path signal 302 ( ^^ is also referred to as “multipath phase delay”). Since the path- length of a multipath signal 304 is always longer than that of the corresponding direct-path signal 302, the time delay ^^ is always larger than zero. The multipath signal amplitude ^^^ is typically smaller than the direct-path signal amplitude ^^ௗ. The multipath phase delay ^^ can be any value ranging from -180° to 180°. When ^^ is equal to 0°, the multipath signal is in- phase with the direct-path signal. When ^^ is equal to 180°, the multipath signal is out-of- phase with the direct-path signal. In some embodiments, ^^^, ^^, and ^^ are the three independent parameters that describe a multipath signal. [0073] The composite signal model in equation 1 assumes that a direct-path signal always exists, i.e., ^^ௗ is always larger than 0. If the direct-path signal is blocked (e.g., by trees or buildings), corresponding multipath signals, if any, may be detected and excluded by the GNSS receiver using multipath mitigation mechanisms different from those described in this document. Equation 1 also assumes that there is at most one multipath signal for each direct-path signal. If more than one multipath signal exists, it is assumed that the strongest multipath signal is a dominant signal that is much stronger than all the remaining multipath signals such that all the remaining multipath signals can be omitted or ignored for purposes of mitigating multipath-induced errors in GNSS receiver positioning. [0074] Figure 4A illustrates a phasor diagram of the in-phase (I) and quadrature (Q) components of a direct-path signal in accordance with some embodiments. Figure 4A shows that this phasor vector is aligned with the I axis, i.e., the value of Q component is equal to 0, and the magnitude of the vector, ^^ௗ 402, is equal to the value of I component. This pair of I and Q samples is just one of the n pairs of coherent integration I and Q samples as shown in Figure 2. As described above, the carrier tracking loop utilizes this pair of I and Q samples to make the phase of the generated carrier signal phase-locked with that of the received signal by driving Q component to zero. In other words, Figure 4A is a phasor diagram when the carrier tracking loop is in a steady state of phase lock. The vector magnitude ^^ௗ 402 may not be numerically equal to the direct-path signal magnitude ^^ௗ in equation 1 due to receiver signal processing; however, they only differ by a scale factor since the receiver signal processing can be regarded as a linear system. [0075] When a direct-path signal 302 and a multipath signal 304 are both received and processed by a receiver, multipath errors can be introduced to range measurements, for example, pseudorange and carrier phase measurements. Figure 4B illustrates a phasor diagram of I (401) and Q (403) components of a composite signal, i.e., a direct-path signal plus a multipath signal, where ^^ௗ 404 is the magnitude of the direct-path signal phasor vector, ^^^ 406 is the magnitude of the multipath signal phasor vector, and ^^ௗା^ 408 is the magnitude of the composite signal phasor vector. When the composite signal is processed by a receiver, the resulted coherent integration I and Q samples are the linear sum of integration values from the direct-path signal and integration values from the multipath signal, corresponding to the sum of the two vectors in the phasor diagram. Thus, at the steady state of a PLL, the summation vector, represented by the dashed line of the composite signal phasor vector, is aligned with (e.g., parallel to) the I axis. However, if the multipath signal does not exist, the direct-path signal phasor vector would be aligned with the I axis instead, just as shown in Figure 4A. Therefore, the angle, ^^ఌ 410 between the summation vector and the direct-path signal vector is the carrier phase multipath error introduced by the multipath signal. An objective of methods and systems disclosed herein is to estimate the carrier phase multipath error ^^ఌ 410, and correct and mitigate multipath errors in the corresponding carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ 410. [0076] Figure 4B shows that the direct-path signal phasor vector, the multipath signal phasor vector, and the composite signal phasor vector form a triangle. Phase difference angle ^^^ 412 between the direct-path signal vector and the multipath signal vector is supplementary to the multipath phase delay ^^ as defined in equation 1. Angle ^^ଷ^ௗ 414 is the third angle of the triangle, so the sum of carrier phase multipath error ^^ఌ 410, the phase difference angle ^^^ 412, and the third angle
Figure imgf000023_0001
414 is equal to 180°, and therefore phase difference angle ^^^ 412 is equal to 180° minus the sum of the other two angles: ^^^ = 180° − ^^ଷ^ௗ − ^^ఌ (2) [0077] If the triangle lies in the first quadrant (i.e., in Quadrant 1416) of this I-Q phasor coordinate system, then all the three angles ^^ఌ 410, ^^^412, and ^^ଷ^ௗ 414, are defined to be positive; otherwise, if the triangle lies in the fourth quadrant (i.e., in Quadrant 4418), then all three angles are defined to be negative. [0078] Based on the geometry in Figure 4B, the following relationships among all these parameters are obtained:
Figure imgf000023_0002
[0079] As explained above, coherent integration ^^ is zero ideally at the steady state of a phase lock loop if noise does not exist, in which case: ^^ௗ sin ^^ఌ = ^^^ sin ^^ଷ^ௗ (4) [0080] Meanwhile, from the geometry in Figure 4B, the following relationship can also be obtained: ^^ௗା^ = ^^ௗ cos ^^ఌ + ^^^ cos ^^ଷ^ௗ (5) [0081] Equation 3 does not contain the multipath time delay ^^, since it assumes that I and Q integration samples are taken when the values of the direct-path signal and the multipath signal both have reached their amplitudes ^^ௗ and ^^^, respectively. When I and Q samples are taken at code chip edge-up or edge-down transitions, equation 3 needs to be changed into the following form:
Figure imgf000023_0003
where ^^ ^^( ^^) is the magnitude step response of the last stage band pass filter in a receiver’s RF front end 174. Equation 6 contains both the step response ^^ ^^( ^^) and its delayed version ^^ ^^( ^^ − ^^), and, as discussed below, can be used to estimate the multipath parameters (e.g., ^^^, ^^, and/or ^^). [0082] In general, a GNSS signal is a sinusoidal carrier signal modulated at least by some code, where the code is a pseudo random sequence of ones and zeros, which are called chips. The signal carrier phase will undergo a 180° phase reversal at a code chip edge transition, either from 0 to 1 or from 1 to 0. For example, GPS L1 Coarse Acquisition (CA) signal has a center carrier frequency of 1575.42 MHz, and its CA code has a chipping rate of 1.023 MHz. Thus, if the baseband DSP 176 carries on correlation and integration for a PDI period of 1 ms, a received GPS L1 CA signal would be expected to have about 1023 code chip edge transitions during each PDI period of 1 ms. [0083] Figure 5 illustrates an example of the magnitude 502 of a filter step response function ^^ ^^( ^^) versus time, in accordance with some embodiments. The dashed line 508 represents an ideal code chip edge-up transition from 0 to 1 instantaneously at time 0. However, because of the RF front end filtering, or equivalently, because the received signal is filtered and thus bandwidth limited, it is impossible to have such an ideal transition which contains infinitely high frequency components in practice. The solid line 506 represents an example of the actual magnitude of the filter step response function ^^ ^^( ^^) versus time, where a code chip edge-up transition serves as a step input function. Compared to the corresponding ideal step response (represented by the dashed line 508), the actual step response ^^ ^^( ^^) shown by the solid line 506 not only has a response time of ^^ 510 but also has some overshoots before and after the transition. These step response characteristics are determined by the type and bandwidth of the final filter in the receiver’s RF front end module 174. For example, a 6-pole Butterworth filter with an IF equivalent bandwidth of 30 MHz has a step response time of about 50 ns and enters steady state after about 150 ns. The code chip edge transition step response model ^^ ^^( ^^) is determined and stored in memory of the receiver before the phase multipath mitigation algorithm (e.g., as described in Figures 9, 11, 12, and 13) is executed. [0084] In some embodiments, the baseband DSP module 176 of a GNSS receiver outputs I and Q samples along code chip edge transitions. When a satellite signal is tracked and locked, code phase management unit 272 can properly control and set the phase delays of n pairs of code replicas in such a way that the corresponding n pairs of correlation results, ^^( ^^^) 259-1I and ^^( ^^^) 259-1Q, ^^( ^^ଶ) 259-2I and ^^( ^^ଶ) 259-2Q, …, and ^^( ^^^) 259-nI and ^^( ^^^) 259-nQ, are along code chip edge transitions. The ith pair of I and Q samples can be rewritten as
Figure imgf000024_0001
where ^^ = 1,2, … , ^^. [0085] Figure 6 illustrates a correlator arrangement in which 16 pairs of samples 602- i (i = 1 to 16 in this example) are arranged (e.g., uniformly) along a code chip transition, and a delayed version of the filter step response function ^^ ^^( ^^ − ^^) 604. In some embodiments, the correlation measurements ^^( ^^^) and ^^( ^^^) are the summed, or equivalently, averaged ones over all code chip edge transitions in one period of correlation duration. In some embodiments, only a subset of all these ^^ pairs of correlators are arranged along code chip edge transitions for the purpose of phase multipath mitigation. In some embodiments, the effective correlators for the purpose of phase multipath mitigation are constructed in a different approach, for example, by utilizing multiple channel resources. In some embodiments, the number of pairs of correlation samples can be increased to cover a larger portion (e.g., even covering a tail portion) of the PN code sequence so as to accurately measure even long delayed multipath signals. Optionally, enough correlators are arranged to cover a entire code chip edge transition. [0086] Given the step response model ^^ ^^( ^^) and ^^ pairs of correlation measurements ^^( ^^^) and ^^( ^^^) along code chip edge transitions, some aspects of the present disclosure estimate the unknown signal parameters (e.g., the carrier phase multipath error ^^ఌ, the multipath time delay ^^, and the amplitude of the multipath signal ^^^) in two procedures, also referred to herein as Procedure 1 and Procedure 2, which are described below (e.g., with respect to Figures 9 and 11). [0087] In some embodiments, due to the influence of a multipath signal, the code chip edge transition of a direct-path signal can be distorted, and accordingly the code chip edge transition response time can be changed from ^^ 510, as shown in Figure 5, to ^^ + ∆ ^^, where the time error ∆ ^^ is unknown but of interest since it largely reveals the code tracking error and code phase measurement error due to the influence of the multipath signal. To identify the unknown time error ∆ ^^, a pair of new measurements ^^ and ^^ are defined as the following:
Figure imgf000025_0001
where t represents time, and ^^′ = ^^ + ∆ ^^ (9) [0088] Before Procedure 1 (as will be described with reference to Figure 9) starts to iteratively solve for the unknown time error ∆ ^^, the relationship between ^^( ^^) and ^^( ^^) needs to be derived first as follows. Based on equations 5 and 6, equation 8 can be rewritten as follows: ^^( ^^) = ( ^^ௗ cos ^^ఌ + ^^^ cos ^^ଷ^ௗ) ^^ ^^( ^^′) − ^ ^^ௗ ^^ ^^( ^^) cos ^^ఌ + ^^^ ^^ ^^( ^^ − ^^) cos ^^ଷ^ௗ^ = ^^ௗ^ ^^ ^^( ^^′) − ^^ ^^( ^^)^ cos ^^ఌ + ^^^^ ^^ ^^( ^^′) − ^^ ^^( ^^ − ^^)^ cos ^^ଷ^ௗ (10) ^^( ^^) = ^^ௗ ^^ ^^( ^^) sin ^^ఌ − ^^^ ^^ ^^( ^^ − ^^) sin ^^ଷ^ௗ (11) [0089] Multiplying equation 10 by
Figure imgf000026_0001
yields ^^( ^^) tan ^^ଷ^ௗ = ^^ௗ ^ ^^ ^^ ( ^^′ ) − ^^ ^^( ^^) ^ cos ^^ఌ tan ^^ଷ^ௗ + ^^^ ^ ^^ ^^ ( ^^′ ) − ^^ ^^( ^^ − ^^) ^ sin ^^ଷ^ௗ (13) [0090] Subtracting equation 13 from equation 11 results in ^^( ^^) = ^^( ^^) tan ^^ଷ^ௗ + ^^ௗ ^^ ^^( ^^) sin ^^ఌ − ^^ௗ^ ^^ ^^( ^^′) − ^^ ^^( ^^)^ cos ^^ఌ tan ^^ଷ^ௗ − ^^^ ^^ ^^( ^^′) sin ^^ଷ^ௗ = ^^( ^^) tan ^^ଷ^ௗ − ^^ௗ^ ^^ ^^( ^^′) − ^^ ^^( ^^)^ cos ^^ఌ tan ^^ଷ^ௗ + ^^^^ ^^ ^^( ^^) − ^^ ^^( ^^′)^ sin ^^ଷ^ௗ = ^^( ^^) tan ^^ଷ^ௗ − ^ ^^ ^^( ^^′) − ^^ ^^( ^^)^( ^^ௗ cos ^^ఌ tan ^^ଷ^ௗ + ^^^ sin ^^ଷ^ௗ) = ^^( ^^) tan ^^ଷ^ௗ − ^ ^^ ^^( ^^′) − ^^ ^^( ^^)^( ^^ௗ cos ^^ఌ + ^^^ cos ^^ଷ^ௗ) tan ^^ଷ^ௗ = ^^( ^^) tan ^^ଷ^ௗ − ^ ^^ ^^( ^^′) − ^^ ^^( ^^)^ ^^ௗା^ tan ^^ଷ^ௗ
Figure imgf000026_0002
where ௗ௧ stands for the value of
Figure imgf000026_0003
at ^^′. On one hand, based on directly given measurement samples of ^^( ^^) and ^^( ^^) at ^^^ with ^^ = 1,2, … , ^^, indirect measurements ^^( ^^) and ^^( ^^) at ^^^ can be calculated as follows:
Figure imgf000026_0004
[0091] On the other hand, equation 14 can also be discretized into the following form:
Figure imgf000026_0005
[0092] The above equation can be equivalently written into the following matrix equation form:
Figure imgf000026_0006
[0093] The goal of Procedure 1 is to solve for ^^ଷ^ௗ and ∆ ^^ from matrix equation 17. Matrix equation 17 can be solved in several different ways. One way is a Least Squares fitting process. Starting with a set of initial estimates for ^^ଷ^ௗ, ∆ ^^, and ^^ௗା^ at iteration k, say, ^^ଷ ^ ^ௗ,^ି^ , ∆ ^ ^^^ି^ and ^ ^ ^ௗା^ , the Least Squares fitting process converts n pairs of ^^( ^^^) and ^^( ^^^) samples into ^^( ^^^) and ^^( ^^^) values based on equation 15, where ^^ = 1,2, … , ^^, and forms the matrix equation 17. Then, the Least Squares fitting process solves for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ from the matrix equation, and generates the k-th estimates of tan ^^ଷ^ௗ and ∆ ^^, i.e., tan ^^^ଷ^ௗ,^ and ∆^ ^^^, respectively, according to tan^ ^^ଷ^ௗ,^ = tan ^^^ଷ^ௗ,^ି^ + ^^థ൫tan ^^ଷ^ௗ − tan ^^^ଷ^ௗ,^ି^ ൯ (18)
Figure imgf000027_0001
where,
Figure imgf000027_0002
are filter gains, which are usually less than 1. Like equations 18 and 19, many other smoothing or filtering methods can also be validly used to filter instantaneous solutions from equation 17. Starting with newly updated estimates tan ^ ^^ଷ^ௗ,^ and ∆ ^ ^^^, the Least Squares fitting process performs a new iteration (e.g., iteration ^^ + 1), if needed (e.g., if the Least Squares fitting process has not yet converged), as discussed below with reference to Figure 9. [0094] Equation 19 indicates that estimate ∆^ ^^^ is the sum of the previous estimate ∆ ^ ^^^ି^ and instantaneous solution ∆ ^^, and then this essentially accumulated estimate ∆ ^ ^^^ is used to convert ^^( ^^^) to ^^( ^^^) using equation 15. In some embodiments, the instantaneous solution ∆ ^^ will be very close to zero as the Least Squares fitting process iterates and converges. Once the Least Squares fitting process converges and ∆ ^^ is very close to zero, both equations 16 and 17 suggest that small errors in estimate ^^^ௗା^ should not be a concern, because the product of ^^^ௗା^ and ∆ ^^ will eventually approach to zero. [0095] Not all multipath scenarios can be handled by equation 17. For example, in some scenarios, the angle ^^ଷ^ௗ is equal or very close to either +90° or -90°, as illustrated in Figure 7A and Figure 7B, respectively. Accordingly, tan ^^ଷ^ௗ will approach ± infinity, and thus solving the corresponding matrix equation 17 will be numerically impossible. In another scenario, the angle ^^ଷ^ௗ is equal or very close to either 0° or 180°, as illustrated Figure 8A and Figure 8B, respectively, and thus tan ^^ଷ^ௗ has a value that is very close to zero. If tan ^^ଷ^ௗ is equal or very close to zero, then all multiplications in the matrix equation 17 will produce values that are zero or close to zero, and the right-hand side ^^( ^^^) of equation 16 should also be zero in theory, and thus solving matrix equation 17 will be numerically unreliable or meaningless. Accordingly, in some embodiments, as discussed below, different equations and algorithms are used to solve for ^^ଷ^ௗ and ∆ ^^ in these scenarios, which cannot be covered by equation 17. [0096] Further, when angle ^^ଷ^ௗ is equal or very close to either +90° or -90°, i.e., cos ^^ଷ^ௗ = 0, the I component of Equation 7 becomes ^^( ^^^) = ^^ௗ ^^ ^^( ^^^) cos ^^ఌ = ^^ூ ^^ ^^( ^^^) (20) where ^^ூ is defined as ^^ூ = ^^ௗ cos ^^ఌ (21) [0097] Equation 20 indicates a mechanism of detecting whether a received satellite signal includes a multipath signal for which ^^ଷ^ௗ = ±90°: if function ^^( ^^^) has the same shape as the chip edge transition model ^^ ^^( ^^^), then this (e.g., the received satellite signal, including both direct and multipath signals) is a multipath scenario of ^^ଷ^ௗ = ±90°. Equation 20 can be equivalently written into the following matrix form:
Figure imgf000028_0001
[0098] Unknown ^^ூ can be easily solved from the above linear matrix equation by, for example, a Least Squares method. With ^^ூ being solved, the Mean Squared Error (MSE) of this estimation can be calculated as
Figure imgf000028_0002
[0099] The smaller the value of ^^ ^^ ^^^^, the higher the similarity between measurement samples ^^( ^^^) and step response model samples ^^ ^^( ^^^). If ^^ ^^ ^^^^ is small enough, then the multipath scenario for the received satellite signal can be treated as a case of ^^ଷ^ௗ = ±90°. [00100] Figure 9 is a workflow 900 for Procedure 1 of the phase multipath mitigation algorithm, in accordance with some embodiments. In some embodiments, the workflow 900 is performed by a GNSS receiver, such as the GNSS receiver 170 in Figure 1B, or by a computing device (e.g., computer system 130, Figure 1C). [00101] The workflow 900 includes obtaining (902) inputs corresponding to n pairs of ^^( ^^^) and ^^( ^^^) along chip edge transitions with ^^ = 1,2, … , ^^ (e.g., corresponding to ^^( ^^^) 259-1I and ^^( ^^^) 259-1Q, ^^( ^^ଶ) 259-2I and ^^( ^^ଶ) 259-2Q, …, and ^^( ^^^) 259-nI and ^^( ^^^) 259-nQ). In some embodiments, the n pairs of ^^( ^^^) and ^^( ^^^) are delivered by integrate and dump circuits 260 as illustrated in Figure 2. Meanwhile, the chip edge transition model ^^ ^^( ^^) or ^^ ^^( ^^^) is pre-determined and therefore available. An example of a pre-determined chip edge transition model ^^ ^^( ^^) or ^^ ^^( ^^^) is been illustrated in Figures 5 and 6. As can be seen in these figures, the amplitude of ^^ ^^( ^^) is normalized to 1. Accordingly, in some embodiments, all correlation samples ^^( ^^^) and ^^( ^^^) are normalized by the same amount before they are input to the phase multipath mitigation algorithm. [00102] The workflow 900 includes determining (904) a sign of carrier phase multipath error ^^ఌ 410. As described with reference to Figure 4B, the sign of ^^ఌ is positive if a triangle formed by the direct-path signal phasor vector, the multipath signal phasor vector, and the composite signal phasor vector is in the first quadrant (e.g., Quadrant 1416 in Figure 4B). Otherwise, the sign of ^^ఌ is negative if the triangle is in the fourth quadrant (e.g., Quadrant 4418 in Figure 4B). Given ^^ samples of ^^( ^^^), one embodiment calculates their magnitudes | ^^( ^^^)|, searches for the maximum, say, ห ^^( ^^^)ห is the maximum, checks the sign of this sample ^^( ^^^), and then the sign of that sample ^^( ^^^) is assigned to the sign of ^^ఌ. The importance of determining the sign of ^^ఌ is described below. [00103] The workflow 900 includes checking (906) the similarity between ^^( ^^^) and ^^ ^^( ^^^) to determine whether the current multipath scenario is a case of ^^ଷ^ௗ = ±90°. This involves solving for ^^ூ from equation 22, calculating the Mean Squared Error, ^^ ^^ ^^^^ using equation 23, and comparing (908) the value of ^^ ^^ ^^^^ to a pre-determined threshold value (e.g., Threshold1). If ^^ ^^ ^^^^ is smaller than (or smaller than or equal to) the threshold (910), it is a multipath case of ^^ଷ^ௗ = ±90° (912). Otherwise, if ^^ ^^ ^^^^ is larger than (or larger than or equal to) the pre-determined threshold (908-No), then ^^ଷ^ௗ is not = ±90°. If the current multipath scenario is detected as a case of ^^ଷ^ௗ = ±90°, then ^^ଷ^ௗ is equal to +90° if the sign of ^^ఌ has been determined as positive; otherwise, ^^ଷ^ௗ is −90° if the sign of ^^ఌ is negative. [00104] If the current multipath scenario is not a case of ^^ଷ^ௗ = ±90°, the workflow 900 continues with the path indicated by step 914, to iteratively solve for ^^ଷ^ௗ and ∆ ^^, using a first iterative process. At the beginning of the first iterative process, an iteration count, e.g., k, representing a number of iterations of the first iterative process is set or reset to an initial value, e.g., 1. In addition, before the first iterative process begins, the values of parameters ^^ଷ^ௗ, ∆ ^^ and ^^ௗା^ are estimated (916) (e.g., initial estimates of these values are set, determined or obtained). In some embodiments, the initial estimates of both ^^ଷ^ௗ and ∆ ^^ are set to 0, i.e., ^^ ^ ଷ^ௗ,^ = 0 and ∆ ^ ^^^ = 0 for the first iteration, i.e., ^^ = 1. As used herein, the term ^^^^ (referred to as “x hat”) denotes an estimate of parameter ^^ in iteration ^^. In some embodiments, ^^^^ is used as the initial estimate for a next iteration (i.e., iteration ^^ + 1) of the first iterative process. [00105] In some embodiments, the initial estimate of ^^ௗା^, i.e., ^ ^ ^ௗା^ , can be set as the magnitude of the last pair of I and Q samples for a code chip edge-up transition, e.g., ^^ ^ௗା^ = ^ ^ ^^ ( ^^^ )^ଶ + ^ ^^ ( ^^^ )^ଶ (24) [00106] Otherwise, in some embodiments, the initial estimate ^ ^ ^ௗା^ can be assigned using the magnitude of the first pair of I and Q samples, e.g., ^^( ^^^) and ^^( ^^^) for a chip edge- down transition. The estimate ^ ^ ^ௗା^ is not updated in the Least Squares fitting process; however, the estimation of ^^ௗା^ using equation 24 is quite accurate and reliable, and as explained above, any error in estimate ^ ^ ^ௗା^ should not be a concern as the Least Squares fitting process converges. [00107] Starting with initial estimates ^^ଷ ^ ^ௗ,^ି^ , ∆ ^ ^^^ି^ and ^ ^ ^ௗା^ at iteration k, the Least Squares fitting process converts n pairs of ^^( ^^^) and ^^( ^^^) samples into ^^( ^^^) and ^^( ^^^) (in step 918) based on equation 15 with ^^ = 1,2, … , ^^, constructs the matrix equation 17, and solves for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ using the Least Squares method. To avoid any potential “divided by zero” error, and more importantly, to check whether the current multipath scenario is a case of ^^ଷ^ௗ = 0° or 180°, the absolute value of tan ^^ଷ^ௗ, i.e., |tan ^^ଷ^ௗ|, is compared to (920) another predefined threshold (e.g., Threshold 2) which is typically a very small positive value. If |tan ^^ଷ^ௗ| is smaller than the predefined threshold (e.g., Threshold2) (922), this is a case of ^^ଷ^ௗ = 0° or 180° (924). Otherwise, if |tan ^^ଷ^ௗ| is larger than the predefined threshold (e.g., Threshold2) (926), then this is not a case of ^^ଷ^ௗ = 0° or 180°. It is noted that this document references multiple thresholds, each of which may be different than the other threshold, and each of which may be set so as to make the multipath mitigation process both efficient and reliable. [00108] If the multipath scenario for a received satellite signal is detected as a case of ^^ଷ^ௗ = 0° or 180°, further determination on whether ^^ଷ^ௗ is equal to 0° or 180° is needed. Recall that Figure 8A shows a multipath scenario of ^^ଷ^ௗ = 0° where a multipath signal is in phase with a direct-path signal and Figure 8B shows a multipath scenario of ^^ଷ^ௗ = 180° where a multipath signal is 180° out of phase with a direct-path signal. As basic knowledge on GPS/GNSS code correlation and multipath, the peak of their composited correlations will shift to right (i.e., be delayed) if a multipath signal is in phase with its direct-path signal. In contrast, the peak of their composited correlations will shift to left (i.e., be advanced) if a multipath signal is 180° out of phase with its direct-path signal. Using this GPS/GNSS knowledge, a mechanism of determining whether ^^ଷ^ௗ is 0° or 180° is devised. As illustrated in Figure 10A, if the peak ^^ூ,^^^^ 1004 of ^^( ^^^) samples 1002 is located to the right-hand side of the peak ^^^ௗ^,^^^^ 1008 of ^^ ^^( ^^^) samples 1006, then ^^ଷ^ௗ = 0°; otherwise, as illustrated in Figure 10B, if the peak ^^ூ,^^^^ 1014 of ^^( ^^^) samples 1012 is located to the left-hand side of the peak ^^^ௗ^,^^^^ 1018 of ^^ ^^( ^^^) samples 1016, then ^^ଷ^ௗ = 180°. In Figure 10A and 10B, the star markers represent ^^( ^^^) measurement samples, and the dot markers represent the code chip edge transition model samples ^^ ^^( ^^^). In either case, the carrier phase multipath error ^^ఌ is 0. In other words, if there is some uncertainty or error (e.g., within predefined margins) in determining whether ^^ଷ^ௗ is either 0° or 180°, the uncertainty or error is tolerable and does not affect the accuracy of estimating the carrier phase multipath error ^^ఌ. [00109] With continued reference to Figure 9, if a multipath scenario is not a case of ^^ଷ^ௗ = 0° or 180° (920-No, 926) and not a case of ^^ଷ^ௗ = ±90° (908-No, 914), the Least Squares fitting process updates (928) the estimate of tan ^^ଷ^ௗ and ∆ ^^ based on equations 18 and 19 to obtain new estimates tan ^^ ^ ଷ^ௗ,^ and ∆ ^ ^^^ at this k th iteration, respectively. The workflow 900 then proceeds to determine (930) whether the Least Squares fitting process has converged, diverged or not. Several different aspects can be considered. In some embodiments, if the magnitude of estimate ∆ ^ ^^^ (i.e., ห∆ ^ ^^^ห) is smaller than a threshold, the workflow 900 determines that the Least Squares fitting process has converged (932), outputs a solution for ^^ଷ^ௗ (934), and terminates the first iterative process. In some embodiments, if the change of estimates of tan ^^ଷ^ௗ in the last two iterations (e.g., หtan ^^^ଷ^ௗ,^ − tan ^^ଷ^^ௗ,^ି^ ห) is smaller than a predetermined threshold, then workflow 900 determines that the Least Squares fitting process has converged (932), outputs a solution for ^^ଷ^ௗ (934), and terminates the iterative process. [00110] In some embodiments, if the number of iterations, ^^, for which the Least Squares fitting process has looped is larger than a threshold (930) (e.g., k > Threshold3), the workflow 900 determines that the iterative process has diverged (or at least cannot converge anymore) and therefore terminates the iterative process without producing a solution for ^^ଷ^ௗ (934). Otherwise (e.g., k < Threshold3 or k ≤ Threshold3), since the process has yet to converge but has not diverged, the workflow 900 repeats the iterative process after increasing the number of iterations ^^ by 1 (938). [00111] If the Least Squares fitting process has successfully converged, the updated estimates of tan ^^ଷ^ௗ and ∆ ^^ obtained from the last iteration are regarded as the final solutions. The value of parameter ^^ଷ^ௗ can be calculated from the solution tan ^^ଷ^ௗ using arctangent function, i.e., ^^ଷ^ௗ = tanି^(tan ^^ଷ^ௗ) (25) [00112] However, arctangent function tanି^(∙) returns an angle in a range of −90° to +90°, while ^^ଷ^ௗ can be in a range of −180° to +180°. To resolve this ambiguity, the sign of carrier phase error ^^ఌ, which was determined in step 904, is utilized again: if the sign of carrier phase error ^^ఌ is positive while ^^ଷ^ௗ calculated from arctangent is negative, then ^^ଷ^ௗ + 180° is re-assigned to the final solution of ^^ଷ^ௗ; similarly, if the sign of ^^ఌ is negative while ^^ଷ^ௗ calculated from arctangent is positive, then ^^ଷ^ௗ − 180° is re-assigned to the final ^^ଷ^ௗ. [00113] Procedure 1 has solved parameter ^^ଷ^ௗ and parameter ∆ ^^, based on the results of the last iteration of Procedure 1 (e.g., using equations 18 and 19), while Procedure 2 (as described in Figures 11, 12, and 13) will continue to solve for the remaining signal parameters, e.g., ^^ௗ, ^^^, ^^^, ^^, and the carrier phase error ^^ఌ. In some instances, Procedure 1 may fail to converge and deliver a no-solution (e.g., no valid solution). If this happens, Procedure 2 will not be executed. [00114] The derivation of equations for Procedure 2 also starts with I and Q measurement model equation 7, which is now rewritten into the following form:
Figure imgf000032_0001
where, ^^ூ is defined in Equation 21, and ^^ொ = ^^ௗ sin ^^ఌ (27) [00115] Equation 26 has four unknowns ^^ூ, ^^ொ, ^^^, and ^^. Since Equation 26 is non- linear, Taylor expansion and linearization techniques are applied here. Given initial estimates of these four unknowns at iteration k, noted as
Figure imgf000032_0002
, ^^ ^ ொ,^ି^ , ^^ ^ ^,^ି^ , and ^ ^ ^^ି^, the Taylor expansion of ^^( ^^^) and ^^( ^^^) at these initial estimates can be expressed as: where, the 2 nd and higher order terms are omitted, ^ ^ ^൫ ^^^,^ି^൯ and ^ ^ ^൫ ^^^,^ି^൯ are the estimated values
Figure imgf000033_0001
Figure imgf000033_0002
^ ^ ^^ି^. Equation 28 can be rewritten into the following form:
Figure imgf000033_0003
where, ^^ூ( ^^^) and ^^ொ( ^^^) are the differences between I & Q measurements and their corresponding estimates, i.e.,
Figure imgf000033_0004
[00116] Equivalently, equation 30 can be written into the following two matrix equations:
Figure imgf000033_0005
[00117] As can be seen, the coefficient matrices in the matrix equations 32 and 33 happen to be same. Once ∆ , ∆^ೂ, and ∆ are solved from equations 32 and estimates of ^^ூ, ^^ொ, ^^^, and ^^ at iteration k can be updated as follows: [00118] Given the updated estimates ^ ^ ^ூ,^ , ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ from iteration k, Procedure 2 can start a new iteration (e.g., iteration k+1) and so on, until the process converges. [00119] Since ∆^^ and ∆ఋ are common unknowns for equations 32 and 33, there are several different approaches for solving these two equations and combining their solutions. One approach is to solve equations 32 and 33 separately, while estimates of ^^ூ, ^^^ and ^^ for equation 32 are updated based on equations 34, 36 and 37, and estimates of ^^ொ, ^^^ and ^^ for equation 33 are updated based on equations 35, 36 and 37. Once the solutions to these two equations both have converged separately, the final solution ^^^ can be a weighted sum of the solution ^^^ from equation 32 and the other ^^^ from equation 33. Similarly, the final solution ^^ can also be a weighted sum of the two ^^ solutions from equations 32 and 33. [00120] Another approach involves converting/combining the two ^^ × 3 matrix equations (equations 32 and 33) into one (2 ^^) × 4 matrix equation with a total of four unknowns (e.g., unknowns ∆^^, ∆^ೂ, ∆^^, and ∆ఋ), and solving the (2 ^^) × 4 matrix equation (e.g., instead of combining solutions). [00121] A third approach for solving equations 32 and 33 and combining their solutions is described here. First, equation 32 is solved to obtain solutions ^^ூ, ^^^ and ^^. Because the unknowns ^^^ and ^^ have been solved, their values are now known. Thus, when ^^( ^^^) in equation 26 is Taylor expanded and linearized, it becomes
Figure imgf000034_0001
where ^ ^ ^൫ ^^^,^ି^൯ is the estimated value of ^^
Figure imgf000034_0002
at the beginning of iteration ^^, i.e.,
Figure imgf000034_0003
[00122] Accordingly, the equation regarding ^^ொ( ^^^,^), which is defined in equation 31, now becomes the following:
Figure imgf000034_0004
or equivalently in a matrix form: [00123] Equation 35 can still be used to update ^^ ^ ொ,^ iteratively. In other words, the third approach solves matrix equation 41 (instead of matrix equation 33) to obtain ^^ொ. The third approach offers two advantages. First, it solves for all possible unknowns (i.e., ^^ூ, ^^^ and ^^) based solely on higher signal-to-noise ratio (SNR) measurements ^^( ^^^) and the corresponding ^^ூ ( ^^^ ) , while measurements ^^( ^^^) and the corresponding ^^ொ ( ^^^ ) are often close to zero in value and have a much lower SNR. Therefore, this approach can deliver more accurate solutions. Second, instead of solving two ^^ × 3 matrix equations (questions 32 and 33) (using the first approach) or a solving (2 ^^) × 4 matrix equation (using the second approach), the third approach solves a ^^ × 3 matrix (equation 32) and a ^^ × 1 matrix (equation 41), thus reducing the amount of computational power used to produce a solution. Furthermore, compared to other larger-sized matrices, an ^^ × 1 matrix equation (e.g., equation 41) typically converges in fewer iterations (e.g., two or three iterations). Therefore, solutions for the unknown parameters can be obtained more quickly and efficiently. [00124] Figure 11 illustrates a workflow 1100 for Procedure 2 of the phase multipath mitigation algorithm, in accordance with some embodiments. In some embodiments, the workflow 1100 is performed by a GNSS receiver (e.g., navigation signal receiver 120 or GNSS receiver 170) or by a computing device (e.g., computer system 130, Figure 1C). Typically, the steps in the workflow 1100 are executed after Procedure 1 (described in the workflow 900) outputs a valid estimation for the carrier phase multipath error ^^ଷ^ௗ. [00125] The workflow 1100 includes determining (1102) initial estimates, ^ ^ ^ூ,^ , ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^, of four unknown parameters, ^^ூ, ^^ொ, ^^^, and ^^. [00126] In some embodiments, the initial estimates ^ ^ ^ூ,^ and ^^ ^ ொ,^ can be assigned with the last pair of correlation measurements, i.e., ^^( ^^^) and ^^( ^^^). [00127] In some embodiments, if Procedure 1 (e.g., the workflow 900) determines that this is a multipath scenario with ^^ଷ^ௗ = ±90°, then ^^^ூ,^ can assigned to the value of ^^ூ that has been solved for from equation 22. As will be shown later, the final solution of ^^ூ is also equal to this value. [00128] In some embodiments, if Procedure 1 determines that this is a multipath scenario with ^^ଷ^ௗ = 0° or 180°, then ^^^ ொ,^ can assigned to 0. As will be shown later, the final solution of ^^ொ in this case is also 0. [00129] In some embodiments, the initial multipath signal amplitude estimate ^^^^,^ can be usually set at 0.5, e.g., half of the amplitude which is used to normalize the chip edge transition model ^^ ^^( ^^) and correlation samples ^^( ^^^) and ^^( ^^^). [00130] In some embodiments, for a very weak direct-path signal partially indicated by a small value of ^^ூ,^, the estimate ^^^^,^ can be scaled down accordingly. Equation 26 shows that ^^( ^^^) and ^^( ^^^) measurement models are linear with parameters ^^ூ, ^^ொ and ^^^, so these initial estimates typically work well. [00131] Equation 26 shows that both ^^( ^^^) and ^^( ^^^) measurement models are highly nonlinear functions in terms of the multipath time delay parameter ^^, and therefore it is critical to make its initial estimate^ ^^^ as accurate as possible. As illustrated in Figures 10A and 10B, assuming the peak sample of the chip edge transition model ^^ ^^( ^^^) is located at time ^^^ௗ^,^^^^ (e.g., ^^^ௗ^,^^^^ 1008 in Figure 10A and ^^^ௗ^,^^^^ 1018 in Figure 10B), and the peak value of ^^( ^^^) measurement samples is found to be at time ^^ூ,^^^^ (e.g., ^^ூ,^^^^ 1004 in Figure 10A and ^^ூ,^^^^ 1014 in Figure 10B), the initial multipath time delay estimate ^ ^^^ (e.g., ^ ^ ^^ 1010 in Figure 10A and ^ ^^^ 1020 in Figure 10B) can be set as
Figure imgf000036_0001
[00132] If later the Least Squares fitting process cannot converge using this estimate ^ ^ ^^, it suggests this estimate ^ ^ ^^ is probably too erroneous. As will be explained later, in Figure 14, there exists a mechanism of trying other values of ^^^^ if needed. Although Figures 10A and 10B are examples of ^^( ^^^) measurement samples when angle ^^ଷ^ௗ is equal to 0° and 180°, respectively, this method of estimating ^^ is valid for all multipath scenarios. [00133] Having determined the initial estimates ^ ^ ^ூ,^ , ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^, Procedure 2 selects an appropriate algorithm to iteratively solve for the unknown parameters, ^^ூ, ^^ொ, ^^^, and ^^, depending on whether ^^ଷ^ௗ is equal to ±90° or not (1104). If Procedure 1 determines that ^^ଷ^ௗ is not equal to ±90° (1104-No, 1106), then the steps in workflow 1200 (Figure 12) are executed. Otherwise, if Procedure 1 determines that ^^ଷ^ௗ is equal to ±90° (1104-Yes, 1108), then the steps in workflow 1300 (Figure 13) are executed. [00134] Figure 12 illustrates a workflow 1200 for solving the unknown parameters, ^^ூ, ^^ொ, ^^^, and ^^ when ^^ଷ^ௗ is not ±90°, in accordance with some embodiments. When ^^ଷ^ௗ is not equal to ±90° (1202), an iterative process at iteration k first estimates (1204) I measurements, ^ ^ ^൫ ^^^,^ି^൯, based on initial estimates ^^ ^ ூ,^ି^ , ^^ ^ ^,^ି^ and ^ ^ ^^ି^ using equation 29, calculates its estimation error
Figure imgf000037_0001
using equation 31, constructs I-matrix equation 32, and solves for unknowns ∆^^, ∆^^ and ∆ఋ using the Least-Squares method (1206). The workflow 1200 updates (1208) ^ ^ ^ூ,^ , ^ ^ ^^,^ and ^ ^ ^^ using equations 34, 36 and 37. Meanwhile, a post-fix ^^ ^^ ^^ூ,^ at each iteration is also calculated as follows:
Figure imgf000037_0002
[00135] At the end of each iteration (e.g., iteration k), the workflow 1200 includes determining (1210) whether the iterative process has diverged or the number of iterations, ^^, is larger than a threshold (e.g., Threshold4). If estimates are out of their valid ranges (e.g., estimate ^ ^ ^^ is larger than the allowably maximal time delay corresponding to the last pair of correlators (i.e., ^^^), or ^ ^ ^^ is less than 0, or estimate ^ ^ ^^,^ is less than 0), then the workflow 1200 determines that the iterative process has diverged (1224). This is because, since a multipath signal must travel through a longer path than its direct-path signal, the multipath time delay ^^ should always be larger than 0. However, if ^ ^ ^^ is larger than ^^^, it means that the multipath signal effect is not observable even by the last pair of correlators. If ^^^^,^ is less than 0, it means multipath signal is either too weak or does not exist. These thresholds may need some fine-tuning or margin, since, for example, ^^^^ may oscillate to a slightly negative value before converging to a small positive one. [00136] With continued reference to Figure 12, if the iterative process has not diverged and ^^ is smaller than a threshold (1212), the workflow 1200 includes determining whether the iterative process has converged (1214). The way of determining convergence here is similar to the one in Procedure 1: if changes in the estimates between two latest iterations are smaller than the corresponding thresholds, and/or if the change in ^^ ^^ ^^ூ, i.e., ห ^^ ^^ ^^ூ,^ − ^^ ^^ ^^ூ,^ି^ห, is smaller than a threshold, and/or if ^^ ^^ ^^ூ,^ itself is smaller than a threshold, then the workflow 1200 determines that the iterative process has converged (1216); otherwise, the iterative process has not converged yet (1220). If the iterative process has converged, the I-matrix equation 32 has been successfully solved, and thus parameters ^^ூ, ^^^ and ^^ are obtained (1218). Otherwise, if the iterative process has not converged yet, a new iteration with ^^ = ^^ + 1 is performed (1222). [00137] Figure 12 also illustrates that in some instances, if the iterative process has diverged or if ^^ is larger than a threshold (1224), then the current iteration is terminated and a no-solution (e.g., no valid solution) is output (1231). In some instances, divergence may be caused by a bad initial estimate ^ ^^^. In these situations, a new, different initial estimate ^ ^ ^^ is tried, as illustrated in steps 1226 and 1228. [00138] In some embodiments, a series of different initial estimates ^^^^ can be tried sequentially (e.g., in a “tree pattern”). This is illustrated in Figure 14. If the first estimate^ ^^^ 1402 results in a divergence, a second new estimate ^ ^^^′ 1404 can be adjusted to its right. If the second estimate ^ ^ ^^′ 1404 also fails, then a third new estimate ^ ^ ^^′′ 1406 with a value smaller than the first estimate ^ ^ ^^ 1402 can be tried. If the third estimate ^ ^^^′′ 1406 also fails, then a fourth estimate ^^^^′′′ 1408 having a value that is larger (e.g., to the right of) than the second estimate ^ ^^^′ 1404 can be tried. In some embodiments, the adjustment step size can be a constant (e.g., equal to the correlator spacing). The range for the adjustments, say, from ^^^^^ 1420 to ^^^^௫1422, should cover all correlation sample times, i.e., from
Figure imgf000038_0001
1424 to ^^^ 1426, with some margin. Once an adjustment moves beyond this range, it means that all possible initial estimates in this corresponding direction have been tried. If the adjustments have covered the whole range, the iterative process terminates with an output that there is no valid solution. This is shown as steps 1230 and 1231 in Figure 12. [00139] Once the workflow 1200 has successfully solved for parameters ^^ூ, ^^^ and ^^ from I-matrix equation 32 (1218), it will continue to solve for the remaining parameter ^^ொ from Q-matrix equation 41. However, if Procedure 1 has determined that ^^ଷ^ௗ is equal to 0° or 180°, then ^^ொ can be directly assigned to 0. See steps 1232, 1234, and 1236. [00140] When ^^ଷ^ௗ is not equal to 0° or 180° (1238), the workflow 1200 proceeds to iteratively solve for ^^ொ using a third iterative process. Note that the iterations for solving equation 41 are independent of iterations for solving equation 32. Before the third iterative process, or at the beginning of the third iterative process (e.g., starting from step 1240), an iteration count, e.g., k, representing a number of iterations of the third iterative process is set or reset to an initial value, e.g., 1. At iteration k, the third iterative process calculates (1240) Q measurements, ^ ^ ^൫ ^^^,^ି^൯, based on the initial estimate
Figure imgf000038_0002
using equation 39, calculates (1240) the measurement estimation error ^^ொ൫ ^^^,^൯ using equation 31, constructs Q-matrix equation 41, solves (1242) for unknowns ∆^ೂ using the Least-Squares method, and updates (1242) ^^ ^ ொ,^ using equation 35. Meanwhile, a post-fix ^^ ^^ ^^ொ,^ at each iteration is also calculated as follows:
Figure imgf000039_0001
[00141] At the end of each iteration (e.g., iteration k), the workflow 1200 determines (1244) whether convergence of the third iterative process has occurred. A way of checking convergence is similar to that described above with respect to step 1214: if the change of estimates of ^^ொ in the last two iterations, e.g.,
Figure imgf000039_0002
is smaller than a threshold, and/or if the change in ^^ ^^ ^^ொ, i.e., ห ^^ ^^ ^^ொ,^ − ^^ ^^ ^^ொ,^ି^ห, is smaller than a threshold, and/or if ^^ ^^ ^^ொ,^ itself is smaller than a threshold, then the third iterative process has successfully converged (1246), and thus ^^ொ is obtained (1248). Otherwise, if the third iterative process has not converged while the number of iterations, ^^, is larger than a threshold (e.g., Threshold5), this means that this process cannot converge and thus should be terminated with an output indicating that there is no valid solution (1248). A third possibility is that the third iterative process has not converged (1250) and that ^^ is smaller than a threshold, then a new iteration with ^^ = ^^ + 1 commences (1252). [00142] To recap, the workflow 1200 in Figure 12 describes the steps of solving for ^^ூ, ^^ொ, ^^^, and ^^ when ^^ଷ^ௗ is not equal to ±90°. In some instances, the iterative process (e.g., the second iterative process) for determining parameters ^^ூ, ^^^ and ^^, or the iterative process (e.g., the third iterative process) for determining ^^ொ can be stopped (e.g., when divergence occurs or when the numbers of iterations has exceeded a threshold number, such as Threshold4 or Threshold5) and information indicating that no valid solution has been determined is output. [00143] Figure 13 illustrates a workflow 1300 for solving for the unknown parameters ^^ூ, ^^ொ, ^^^, and ^^ when ^^ଷ^ௗ is equal to ±90° (step 1302, Figure 13), in accordance with some embodiments. As discussed above, because Procedure 1 (workflow 900) has determined that ^^ଷ^ௗ is equal to ±90° using ^^( ^^^) measurements based on equation 22, which is used to solve for ^^ூ, there is no need to repeat this process of solving for some unknowns from ^^( ^^^) measurements. Instead, the solution ^^ூ to equation 22 obtained in Procedure 1 can be assigned not only to the initial estimate
Figure imgf000039_0003
but also to the final solution of ^^ூ. The remaining steps in the workflow 1300 solve for ^^ொ, ^^^ and ^^ using matrix equation 33. In some embodiments, the steps for solving matrix equation 33, which are shown in the workflow 1300, are similar to those for solving matrix equation 32, which are shown on the left portion of the workflow 1200. [00144] The workflow 1300 in Figure 13 is an iterative process (e.g., a fourth iterative process) that includes calculating (1304) (e.g., at iteration k) Q measurements,
Figure imgf000040_0001
based on initial estimates ^^ ^ ொ,^ି^ ,
Figure imgf000040_0002
and ^ ^ ^^ି^ using equation 29, calculating (1304) the estimation error ^^ொ൫ ^^^,^൯ using equation 31, constructing Q-matrix equation 33, and solving (1306) for unknowns ∆^ೂ, ∆^^ and ∆ఋ using the Least-Squares method. Before the fourth iterative process, or at the beginning of the fourth iterative process (e.g., at step 1302), an iteration count, e.g., k, representing a number of iterations of the four iterative process is set or reset to an initial value, e.g., 1. The workflow 1300 includes updating (1308) ^^ ^ ொ,^ , ^ ^ ^^,^ and ^ ^ ^^ using equations 35, 36 and 37. Meanwhile, a post-fix ^^ ^^ ^^ொ,^ at each iteration is also calculated using equation 44. [00145] At the end of each iteration (e.g., iteration k) of the fourth iterative process, the workflow 1300 determines (1310) whether the iterative process has diverged or the number of iterations, ^^, is larger than a threshold (e.g., Threshold6). If estimates are out of their valid range (e.g., if estimate ^ ^ ^^ is larger than the allowably maximal time delay corresponding to the last pair of correlators, ^^^, or if ^ ^ ^^ is less than 0, or if estimate ^ ^ ^^,^ is less than 0, then the iterative process is determined to have diverged (1320). In some embodiments, these thresholds may be fine-tuned to some margin, since ^ ^ ^^ may oscillate to a slightly negative value before converging to a small positive one. [00146] If the fourth iterative process has not diverged and the iteration count, ^^, for the fourth iterative process is smaller than a threshold (1312), the workflow 1300 determines whether the fourth iterative process has converged (1314). For example, if changes in the estimates between two latest iterations are smaller than the corresponding thresholds, and/or if the change in ^^ ^^ ^^ொ, i.e.,
Figure imgf000040_0003
is smaller than a threshold, and/or if ^^ ^^ ^^ொ,^ itself is smaller than a threshold, then the workflow 1300 determines that the fourth iterative process has converged (1316). Otherwise, the workflow 1300 determines that the fourth iterative process has not converged. If the fourth iterative process has converged (1316), the Q-matrix equation 33 has been successfully solved, and the ^^ொ, ^^^ and ^^ are obtained (1318). Otherwise, the fourth iterative process has not converged (1328) and a new iteration with ^^ = ^^ + 1 (1330) is executed (1330). [00147] In some embodiments, if the fourth iterative process has diverged or if ^^ is larger than a threshold (1320), then the current iteration is terminated, as illustrated in steps 1324 and 1326. In some embodiments, the workflow 1300 includes re-solving equation 33 using a new, different estimate^ ^^^ (1322). For example, in some embodiments, Procedure 2 selects different initial estimates ^^^^ using the approach that is explained above with respect to Figure 14. In some embodiments, if all candidates of initial estimates ^ ^ ^^ have been tried, the fourth iterative process is terminated and the workflow 1300 outputs information indicating that no valid solution has been obtained (1326). [00148] Referring again to Figure 11, at this point, Procedure 2 has solved for parameters ^^ூ, ^^ொ, ^^^, and ^^, regardless of whether ^^ଷ^ௗ is equal to ±90° or not. It is also possible that no valid solution is obtained in Procedure 1 and/or Procedure 2 (steps 110, 1112 and 1114). In some embodiments, if there is no valid solution, this means that no error correction to carrier or code phase measurements is available. The multipath error (if any) cannot be mitigated, and the receiver’s performance is the same as that without phase multipath mitigation. According to some embodiments of the present disclosure, the phase multipath mitigation techniques output information indicating that there is not a valid solution, instead of a suspiciously bad solution, to protect GNSS receivers’ reliability performance. [00149] With continued reference to Figure 11, if a valid solution for parameters ^^ூ, ^^ொ, ^^^, and ^^ is delivered (1115), all the remaining parameters can be calculated as well. Steps 1116, 1118 and 1120 of Figure 11 show that when ^^ଷ^ௗ is not equal to ±90°, the carrier phase measurement error ^^ఌ can be calculated as follows:
Figure imgf000041_0001
[00150] Steps 1116, 1124, and 1126 of Figure 11 show that when ^^ଷ^ௗ is equal to ±90°, error ^^ఌ can be calculated as follows instead:
Figure imgf000041_0002
[00151] Regardless of whether ^^ଷ^ௗ is equal to ±90° or not, the amplitude ^^ௗ of the direct-path signal is: ^^ ௗ = ^ ^^ଶ ூ + ^^ଶ ொ (47) [00152] And ^^^ can be calculated (1122) using equation 2. As a reminder, angles ^^ଷ^ௗ, ^^ఌ and ^^^ are all positive if a triangle formed by three phasor vectors is in Quadrant 1 416, as illustrated in Figure 4B; otherwise, they are all negative if the triangle is in Quadrant 4418. [00153] So far, all of the unknown parameters, including ^^ଷ^ௗ, ^^ఌ, ^^^, ^^ூ, ^^ொ, ^^^, ^^ௗ, ∆ ^^, and ^^, have been solved. The carrier phase measurement error ^^ఌ can be used to correct the corresponding carrier phase measurement, and thus the carrier phase multipath error is mitigated. The time error ∆ ^^, if solved in Procedure 1, can be used as a reference to correct the corresponding code phase measurement and pseudorange measurement. In some embodiments, these parameters, whose values have been determined, can be further used, for example, to mitigate other errors using other signal tracking and multipath mitigation techniques. [00154] Some applications solve for the carrier phase multipath error ^^ఌ and/or the time error ∆ ^^. In these instances, the steps (described above) that are used to calculate other parameters can be skipped as soon as ^^ఌ and ∆ ^^ are obtained. For example, if Procedure 1 determines that ^^ଷ^ௗ is equal to 0° or 180°, then the entire Procedure 2 can be skipped and a correct solution to ^^ఌ can be directly delivered ( ^^ఌ is equal to 0° in this case). Thus, depending which of the unknown parameters are of interest, computational power and time can be saved. [00155] The above procedures can also work in instances where a multipath signal is weak or does not exist. Considering weak multipath or no-multipath scenarios, for example when ^^^ is close to 0, the I component of equation 7 will become equation 20. In other words, weak multipath and no-multipath scenarios will be detected as cases of ^^ଷ^ௗ = ±90° in Procedure 1. [00156] Figure 15A illustrates a phasor diagram for a very weak multipath scenario with ^^ଷ^ௗ = 30°. The geometry in Figure 15A is assumed to reflect a weak multipath scenario, where the multipath signal amplitude ^^^ ≪ ^^ௗ, and the value of the angle ^^ଷ^ௗ (“true value”) is equal to 30°. [00157] Figure 15B illustrates a correct solution to the carrier phase multipath error ^^ఌ even when the multipath scenario is detected as a case of ^^ଷ^ௗ = +90°. In this example, ^^ଷ^ௗ is detected to be 90° by Procedure 1, and then the steps described in Figure 13 are executed, where low SNR measurements ^^( ^^^) are utilized, and these nearly zero-valued measurements ^^( ^^^) result in small values of solutions ^^^ and ^^ఌ. The comparison between Figure 15A and Figure 15B demonstrates that the solution for the carrier phase multipath error ^^ఌ in Figure 15B is quite close to the ^^ఌ carrier phase multipath error in Figure 15A, even though the solution ^^ଷ^ௗ could be very different from the “true value” of ^^ଷ^ௗ. [00158] In some instances, weak multipath and no-multipath scenarios can also possibly be detected as cases of ^^ଷ^ௗ = 0° or 180°, if they are not detected as cases of ^^ଷ^ௗ = ±90° in Procedure 1. If this happens, Figure 12 shows that solution ^^ொ will be equal to 0 for cases of ^^ଷ^ௗ = 0° or 180°, and then equation 45 will result in a correct solution of ^^ఌ = 0°. [00159] Now that the phase multipath mitigation algorithm has been fully described, the code chip edge transition model ^^ ^^( ^^) needs to be further discussed. Recall that Figure 6 illustrates a code chip edge transition model ^^ ^^( ^^). If the n data marks 602-i in Figure 6 (representing n pairs of correlation measurements) only cover a portion of the transition, this will result in an ambiguity between multipath time delay ^^ and multipath signal amplitude ^^^. In other words, this phase multipath mitigation may converge to a wrong solution or even diverge. Therefore, it is recommended that n pairs of correlation measurements should be arranged along the whole code chip edge transition including its tail. [00160] In addition to a code chip edge transition model ^^ ^^( ^^), Figure 6 also illustrates its time-delayed version ^^ ^^( ^^ − ^^). Figure 6 shows that ^^ ^^( ^^ − ^^) and the corresponding multipath signal may not be observable even by the last pair of receiver correlators if multipath time delay ^^ is too large, larger than the last correlator pair’s code phase ^^^. In order to mitigate a long-delayed multipath, more correlator resources are needed in receiver hardware design to cover and measure the very tail part of the code chip edge transition model ^^ ^^( ^^), if correlator spacing remains unchanged. [00161] On the other hand, when multipath time delay ^^ is very small, the function ^^ ^^( ^^) and ^^ ^^( ^^ − ^^) will almost overlap in time, ^^( ^^^) and ^^( ^^^) measurements of a direct- path signal and of a multipath can have a very similar pattern. Therefore, it can be very challenging for this phase multipath mitigation technique to estimate multipath parameters based on ^^( ^^^) and ^^( ^^^) measurements of the corresponding composite signal. More correlators need to be densely arranged along code chip edge transitions if the system needs to mitigate a short-delayed multipath. [00162] When the phase multipath mitigation system needs to deal with very short- and very long-delayed multipath scenarios, it demands more correlator resources. This tradeoff between receiver hardware complexity and multipath mitigation performance should be considered in GNSS receiver design. In some instances. when the receiver hardware does not have enough resources to support multipath mitigation for some multipath scenarios, the phase multipath mitigation algorithm that is described in Procedure 1 and Procedure 2 should determine the MSE values for most if not all steps, as described above, and deliver a no- solution as needed, to ensure that the receiver’s reliability and integrity performance are protected. [00163] Figures 16A and 16B illustrate a flowchart diagram for a method 1600 of mitigating the effect of a multipath-induced error in a global navigation satellite system (GNSS), in accordance with some embodiments. In some embodiments, the method 1600 is performed at a GNSS receiver (e.g., navigation signal receiver 120 in Figure 1A or GNSS receiver 170 in Figure 1B), or at a computing device (e.g., computer system 130) that includes one or more processors (e.g., CPU(s) 202) and memory (e.g., memory 210). For ease of explanation, method 1600 will be explained as being performed by a GNSS receiver, but it will be understood that at least some portions of method 1600 could be performed by an external system. [00164] The GNSS receiver receives (1602) a band-limited composite signal corresponding to a respective satellite in the GNSS, including a band-limited direct-path signal and a band-limited multipath signal. The direct-path signal and the multipath signal are modulated with a pseudorandom (PN) code [00165] The GNSS receiver obtains (1604), for a code chip edge transition of the PN code (e.g., a chip edge-up transition or a chip edge-down transition), n pairs of correlation samples (e.g., samples 602-i, Figure 6) of the composite signal during a corresponding n time instances, each having a different time offset from the code chip edge transition. In some embodiments, n is a positive integer such as 16, 32, 48, or the like. The code chip edge transition has a predetermined filter step response function SR(t) (see, e.g., Figure 5). Each respective pair of correlation samples of the n pairs of correlation samples consists of a respective in-phase component sample I(ti) and a respective quadrature component sample Q(ti), wherein i is a positive integer from one to n. The n pairs of correlation samples consist of n in-phase component samples of the composite signal (I(t) samples) and n quadrature component samples of the composite signal (Q(t) samples). [00166] In some embodiments, the n pairs of correlation samples are uniformly arranged along the code chip edge transition. In some embodiments, the n pairs of correlation are arranged over one period of correlation duration. In some embodiments, a respective pair of correlation samples (I(ti), Q(ti)) is summed or averaged over multiple periods of correlation duration. In some embodiments, the number of the pairs of correlation samples can be increased to cover a larger portion (e.g., even covering a tail portion) of the PN code sequence so as to accurately measure even long delayed multipath. [00167] The GNSS receiver obtains (1606), for each respective pair of correlation samples of the n pairs of correlation samples, a corresponding sample of the filter step response function SR(ti) corresponding in time with the respective pair of correlation samples, thereby obtaining n samples of the filter step response function (predetermined SR(t) samples) during the corresponding n time instances. [00168] In accordance with a determination that a computed similarity between the I(t) samples and the SR(t) samples satisfies a first threshold value, the GNSS receiver determines (1608) that a phase ^^ଷ^ௗ is ±90° within a first predefined margin (e.g., no greater than 1, 2, 3, 4 or 5 degrees). The phase ^^ଷ^ௗ is equal to 180° minus a sum of (i) a carrier phase multipath error ^^ఌ (e.g., which can be used to correct carrier phase measurement) and (ii) a phase difference ^^^ between the direct-path signal and the multipath signal (see Equation 2). In some embodiments, ^^^ is equal to 180 degrees minus the difference between the direct-path signal and the multipath signal. [00169] In accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples does not satisfy the first threshold value, the GNSS receiver solves (1610) a first set of matrix equations (e.g., Equation 17), thereby obtaining a solution for tan ^^ଷ^ௗ and for ∆ ^^ tan ^^ଷ^ௗ, wherein ∆ ^^ is a response time error of the code chip edge transition due to the multipath signal. [00170] In accordance with a determination that | ^^ ^^ ^^ ^^ଷ^ௗ| satisfies a second threshold value (e.g., | ^^ ^^ ^^ ^^ଷ^ௗ | < Threshold2 or | ^^ ^^ ^^ ^^ଷ^ௗ | ≤ Threshold2, as illustrated in Figure 9, step 920), the GNSS receiver determines (1612) that the phase ^^ଷ^ௗ is 0° or 180° within a second predefined margin (e.g., no greater 1, 0.5, 0.1, 0.01, or 0.001 degrees). In some embodiments, the second threshold value (e.g., Threshold2 in Figure 9, step 920) is a small positive value to avoid any potential “divided by zero” error. In some embodiments, the second threshold value is a value such as 0.01 or less, with some examples being 0.0001, 0.001, 0.005, or 0.01. [00171] In accordance with a determination that | ^^ ^^ ^^ ^^ଷ^ௗ | does not satisfy the second threshold value (see Figure 9, step 926), the GNSS receiver determines (1614) ^^ଷ^ௗ in accordance with the solution for tan ^^ଷ^ௗ (see above discussion of Figure 9, steps subsequent to step 926). For example, in some embodiments, | ^^ ^^ ^^ ^^ଷ^ௗ| does not satisfy the threshold value when | ^^ ^^ ^^ ^^ଷ^ௗ | > Threshold2 or | ^^ ^^ ^^ ^^ଷ^ௗ | ≥ Threshold2. [00172] The GNSS receiver adjusts (1616) a pseudorange measurement for the respective satellite in accordance with the determined ∆ ^^; and/or adjusts a carrier phase measurement for the respective satellite in accordance with a parameter (e.g., the carrier phase multipath error ^^ఌ) corresponding to the determined ^^ଷ^ௗ. [00173] The GNSS receiver performs (1618) a navigation function (e.g., determines a location of the GNSS receiver) using the adjusted pseudorange and/or carrier phase measurements for the respective satellite. [00174] In some embodiments, the GNSS receiver determines, based on the Q(t) samples, a sign (e.g., positive or negative) of the carrier phase multipath error ^^ఌ. For example, in some embodiments, if all ^^( ^^^) samples are positive, then the phasor diagram is in the first quadrant (e.g., Quadrant 1416, Figure 4B) and the sign of the carrier phase multipath error ^^ఌ is positive. In another example, if all ^^( ^^^) samples are negative, the phasor diagram is in the fourth quadrant (e.g., Quadrant 4418) and the sign of the carrier phase multipath error ^^ఌ is negative. In some embodiments, some ^^( ^^^) are positive while some are negative. In this instance, the GNSS receiver finds the maximum of their magnitudes, say |Q( ^^^)| is maximal, then the phasor diagram is in the first quadrant if Q(tj) is positive and the phasor diagram is in the fourth quadrature if Q(tj) is negative. [00175] In some embodiments, the GNSS receiver determines a sign of the carrier phase multipath error ^^ఌ (e.g., positive or negative). In accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples satisfies (e.g., is less than, or is less than or equal to) the first threshold value (see Figure 9, steps 906 and 908), the GNSS receiver determines that ^^ଷ^ௗ is +90° within the first predefined margin when the sign of the carrier phase multipath error ^^ఌ is positive, and determines that ^^ଷ^ௗ is -90° within the first predefined margin when the sign of the carrier phase multipath error ^^ఌ is negative. [00176] In some embodiments, in accordance with a determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin, the GNSS receiver determines a position of a peak of the I(t) samples relative to a position of a peak of the predetermined SR(t) samples. In accordance with a determination that the position of the peak of the I(t) samples (e.g., ^^ூ,^^^^ 1004 in Figure 10A) occurs later in time than the position of the peak of the SR(t) samples (e.g., ^^^ௗ^,^^^^ 1008 in Figure 10A), the GNSS receiver determines that ^^ଷ^ௗ is 0° within the second predefined margin. In accordance with a determination that the position of the peak of the I(t) samples (e.g., ^^ூ,^^^^ 1014) occurs earlier in time than the position of the peak of the predetermined SR(t) samples (e.g., ^^^ௗ^,^^^^ 1018), the GNSS receiver determines that ^^ଷ^ௗ is 180° within the second predefined margin. As described above with respect to Figures 10A and 10B, when ^^ଷ^ௗ = 0°, the multipath signal is in phase with the direct-path signal and the peak of their composited correlations will shift to the right (i.e., get delayed). When ^^ଷ^ௗ = 180°, the multipath signal is out of phase with the direct- path signal, and the peak of their composited correlations will shift to the left (e.g., be advanced). [00177] In some embodiments, the computed similarity between the I(t) samples and the SR(t) samples is obtained using a mean squared error (MSE), defined by: ^^ ^^ ^^^^ = ^ ^ ∑^ ^ୀ^ ^ ^^ ( ^^^ ) − ^^ூ ^^ ^^( ^^^) ^ଶ (see Equation 23). Furthermore, in such embodiments, ^^ூ = ^^ௗ cos ^^ఌ, where ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error. [00178] In some embodiments, solving the first set of matrix equations comprises: obtaining (i) an initial estimate ^^ ^ ଷ^ௗ,^ for the phase ^^ଷ^ௗ, (ii) an initial estimate ∆ ^ ^^^ for the response time error ∆ ^^, and (iii) an initial estimate ^ ^ ^ௗା^ for a magnitude ^^ௗା^ of the composite signal (See Figure 9, step 916). For example, in some embodiments, the initial estimate ^^ ^ ଷ^ௗ,^ is set to zero. In some embodiments, the initial estimate ∆ ^ ^^^ is set to zero. In some embodiments, the initial estimate ^^^ௗା^ is set as the magnitude of the last pair of correlation samples of the n pairs of correlation samples for code chip edge-up transitions, i.e., ^ ^ ^ௗା^ = ^ ^ ^^ ( ^^^ )^ଶ + ^ ^^ ( ^^^ )^ଶ (see Equation 24). In some embodiments, the initial estimate ^^^ௗା^ is set as the magnitude of the first pair of correlation samples of the n pairs of correlation samples for code chip edge-down transitions, i.e., ^^^ௗା^ = ^^ ^^( ^^^)^ଶ + ^ ^^( ^^^)^ଶ). In some embodiments, the GNSS receiver uses the initial estimates ^^ ^ ଷ^ௗ,^ , ∆ ^ ^^^, and ^ ^ ^ௗା^ to solve the first set of matrix equations (e.g., Equation 17) via a Least Squares fitting process to obtain the solution for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ. [00179] In some embodiments, in accordance with the determination that | ^^ ^^ ^^ ^^ଷ^ௗ | does not satisfy the second threshold value (see, e.g., Figure 9, step 926), the GNSS receiver obtains a solution for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ via a first iterative process. Each iteration of the first iterative process includes: obtaining an updated estimate tan ^^ ^ ଷ^ௗ,^ for the phase ^^ଷ^ௗ and an updated estimate ∆^ ^^^ for the response time error ∆ ^^ (e.g., according to Equations 18 and 19) (see, e.g., Figure 9, step 928). In accordance with a determination that a number of iterations (e.g., the current number of iterations, k) for the first iterative process (e.g., iterations completed so far) satisfies a third threshold value (e.g., k > Threshold3 or k ≥ Threshold3 as illustrated in Figure 9, steps 930 and 932): the GNSS receiver outputs information that no valid solution for ^^ଷ^ௗ has been determined (see, e.g., Figure 9, step 934) and terminates the first iterative process without producing a solution for ^^ଷ^ௗ. For example, in some embodiments, if the number of iterations, ^^, for which the least squares fitting process satisfies the third threshold, the first iterative process has diverged (or at least cannot converge anymore) and therefore should be cancelled with a no valid solution. In some embodiments, each iteration of the first iterative process further includes: in accordance with a determination that the number of iterations for the first iterative process does not satisfy the third threshold value (e.g., k < Threshold3 or k ≤ Threshold3) (e.g., the first iterative process has not diverged, or the number of iterations, ^^, for the least squares fitting process does not equal or exceed a maximum value): in accordance with a determination that a magnitude based on the updated estimate ∆ ^ ^^^ (e.g., ห∆ ^ ^^^ห) or based on the updated estimate tan ^^ ^ ଷ^ௗ,^ (e.g., a difference in magnitude between the updated estimate tan ^^^ଷ^ௗ,^ and the estimate tan ^^ଷ^^ௗ,^ି^ from the previous iteration,
Figure imgf000048_0001
ห) satisfies (e.g., is less than, or is less than or equal to) a fourth threshold value: determining that the Least Squares fitting process has converged (Figure 9, step 932); outputting the updated estimate tan ^^ ^ ଷ^ௗ,^ and the updated estimate ∆^ ^^^ as the solutions for tan ^^ଷ^ௗ and ∆ ^^, respectively (e.g., Figure 9, step 934); calculating the phase ^^ଷ^ௗ from the updated estimate tan ^^ ^ ଷ^ௗ,^ (i.e., ^^ଷ^ௗ = ^^ ^^ ^^ି^൫tan
Figure imgf000048_0002
(Figure 9, Step 934); and terminating the first iterative process. In some embodiments, each iteration of the first iterative process further includes: in accordance with a determination that the magnitude based on the updated estimate ∆^ ^^^ or based on the updated estimate tan ^^^ ଷ^ௗ,^ does not satisfy the fourth threshold value, performing a next iteration of the first iterative process (see., e.g., Figure 9, step 938). For example, in some embodiments, the magnitude based on the updated estimate ∆^ ^^^ or based on the updated estimate tan ^^ ^ ଷ^ௗ,^ does not satisfy the fourth threshold value when ห∆ ^ ^^^ห is greater than, or is greater than or equal to, the fourth threshold value, or when
Figure imgf000048_0003
is greater than, or is greater than or equal to, the fourth threshold value. [00180] In some embodiments, in accordance with a determination that a valid solution for the phase ^^ଷ^ௗ has been determined (e.g., when ^^ଷ^ௗ is one of: ±90° within the first predefined margin, 0° or 180° within the second predefined margin, or the updated estimate tan ^^^ଷ^ௗ,^ ), the GNSS receiver determines (i) an initial estimate (e.g., a value) ^^^ூ,^ for ^^ூ, wherein ^^ூ = ^^ௗ cos ^^ఌ (see Equation 21), ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error; (ii) an initial estimate ^^^ ொ,^ for a magnitude ^^ொ, wherein ^^ொ = ^^ௗ sin ^^ఌ (see Equation 27); (iii) an initial estimate ^^^^,^ for a magnitude ^^^ of the multipath signal, and (iv) an initial estimate ^ ^ ^^ for a time delay ^^ of the multipath signal relative to the direct-path signal (see, e.g., Figure 11, step 1102). The GNSS receiver uses at least a subset of (i) the initial estimate ^ ^ ^ூ,^ , (ii) the initial estimate ^^ ^ ொ,^ , (iii) the initial estimate ^ ^ ^^,^, and (iv) the initial estimate ^ ^^^ as starting values, performing one or more iterative processes to obtain solutions for at least a subset of ^^ூ, ^^ொ, ^^^, and ^^ (e.g. using Procedure 2, as illustrated in Figures 11, 12, and 13). The GNSS receiver computes the carrier phase multipath error ^^ఌ in accordance with the obtained solutions for at least the subset of ^^ூ, ^^ொ, ^^^, and ^^, and corrects the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [00181] In some embodiments, in accordance with a determination that a valid solution for the phase ^^ଷ^ௗ has been determined (e.g., when ^^ଷ^ௗ is one of: ±90° within the first predefined margin, 0° or 180° within the second predefined margin, or the updated estimate tan ^^^ଷ^ௗ,^ ), the GNSS receiver determines (i) an initial estimate (e.g., a value) ^^^ூ,^ for ^^ூ, wherein ^^ூ = ^^ௗ cos ^^ఌ (see Equation 21), ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error; (ii) an initial estimate ^^ ^ ொ,^ for a magnitude ^^ொ, wherein ^^ொ = ^^ௗ sin ^^ఌ (see Equation 27); (iii) an initial estimate ^^^^,^ for a magnitude ^^^ of the multipath signal, and (iv) an initial estimate ^ ^ ^^ for a time delay ^^ of the multipath signal relative to the direct-path signal. This is also illustrated in Figure 11, step 1102. [00182] In some embodiments, in accordance with a determination that the code chip edge transition is a code chip edge-up transition (and in accordance with the determination that the phase ^^ଷ^ௗ does not correspond to ±90° within the first predefined margin): the initial estimate ^ ^ ^ூ,^ corresponds to a magnitude of a last in-phase component sample I(tn) of the I(t) samples; and the initial estimate ^^^ ொ,^ corresponds to a magnitude of a last quadrature component Q(tn) of the Q(t) samples. [00183] In some embodiments, in accordance with the code chip edge transition is a code chip edge-down transition (and in accordance with the determination that the phase ^^ଷ^ௗ does not correspond to ±90° within the first predefined margin): the initial estimate ^^^ூ,^ corresponds to a magnitude of a first in-phase component sample I(t1) of the I(t) samples; and the initial estimate ^^^ொ,^ corresponds to a magnitude of a first quadrature component Q(t1) of the Q(t) samples. [00184] In some embodiments, in accordance with the determination that the phase ^^ଷ^ௗ corresponds to ±90° within the first predefined margin, the GNSS receiver obtains a solution for ^^ூ by solving a matrix equation (see Equation 22), and
Figure imgf000050_0001
assigning the obtained solution as the initial estimate ^ ^ ^ூ,^ . [00185] In some embodiments, in accordance with the determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin, the GNSS receiver assigns the initial estimate ^^ ^ ொ,^ as zero. [00186] In some embodiments, the initial estimate ^ ^ ^^,^ is 0.5 (e.g., half of the amplitude that is used to normalize the chip edge transition model ^^ ^^( ^^) and correlation samples ^^( ^^^) and ^^( ^^^)). In some embodiments, for a very weak direct-path signal partially indicated by a small value of ^^ூ,^, estimate ^ ^ ^^,^ can be scaled down accordingly. [00187] In some embodiments, in accordance with the determination that the phase ^^ଷ^ௗ does not correspond to ±90° within the first predefined margin (see, e.g., Figure 11, step 1106), the GNSS receiver uses the initial estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as starting values (e.g., obtained using equation 29) to obtain a solution for ^^ூ, ^^^, and ^^ via a second iterative process (see, e.g., workflow 1200 in Figure 12). At the beginning of the second iterative process, an iteration count, e.g., k, representing a number of iterations of the second iterative process is set or reset to an initial value, e.g., 1. A kth iteration of the second iterative process includes: obtaining n estimated in-phase component values ( ^ ^ ^ ( ^^^ି^ ) values) for the I(t) samples based on estimates ^^ ^ ூ,^ି^ , ^^ ^ ^,^ି^ , and ^ ^ ^^ି^ from a previous (k-1) th iteration of the second iterative process (e.g., using equation 29: ^ ^ ^( ^^^,^ି^) = ^^ ^ ூ,^ି^ ^^ ^^( ^^^) + ^^ ^ ^,^ି^ ^^ ^^( ^^^ − ^ ^ ^^ି^) cos ^^ଷ^ௗ) (see, e.g., Figure 12, step 1204), determining, for each estimated in-phase component value ^ ^ ^൫ ^^^,^ି^൯ of the ^ ^ ^( ^^^ି^) values, a corresponding in-phase component estimation error ^^ூ൫ ^^^,^൯ between (i) the respective estimated in-phase component value
Figure imgf000051_0001
and (ii) the respective in-phase component sample I
Figure imgf000051_0002
the I(t) samples (e.g., using equation 31:
Figure imgf000051_0003
), thereby obtaining n in-phase component estimation errors ( ^^ூ ( ^^^ ) estimation errors) (see, e.g., Figure 12, step 1204); solving a second set of matrix equations (e.g., Equation 32) based on the ^^ூ ( ^^^ ) estimation errors, thereby obtaining (e.g., using the results produced by solving the second set of matrix equations) updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ for the k th iteration (e.g.., using equations 34, 36 and 37) (see, e.g., Figure 12, step 1208). In accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is within a corresponding first valid range (e.g., the updated estimate ^ ^ ^^,^ is greater than 0) and (ii) the updated estimate ^^^^ is within a corresponding second valid range (e.g., ^ ^ ^^ is within the allowably maximal time delay corresponding to the last pair of correlators, ^ ^ ^^ is greater than 0 but less than ^^^) (i.e., the second iterative process has not diverged): in accordance with a determination that the number of iterations ^^ satisfies a fifth threshold value (e.g., k satisfies the fifth threshold value when k < Threshold4 or k ≤ Threshold4, as illustrated in Figure 12, steps 1210 and 1212): in accordance with a determination that a change in value of one or more predetermined first parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold, determining that the second iterative process has converged (see, e.g., Figure 12, step 1216); outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively (see, e.g., Figure 12, step 1218); and terminating the second iterative process. [00188] In some embodiments, the one or more predetermined first parameters include one or more of ^^ூ, ^^^, and ^^. In some embodiments, the determination that the second iterative process has converged is in accordance with a determination that (i) the change in an estimated value of ^^ூ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e.,
Figure imgf000051_0004
^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^^ூ), and/or (ii) the change in an estimated value of ^^^ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e., ห ^ ^ ^ெ,^ − ^^ ^ ெ,^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^^ெ), and/or (iii) the change in an estimated value of ^^ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e., ห ^ ^ ^^ − ^ ^ ^^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ఋ). In some embodiments, the one or more predetermined first parameters include a post-fix mean squared error (MSE) for the in-phase component estimation error, defined by the kth iteration (see Equation 43). In some embodiments, the determining that the second iterative process has converged is in accordance with a determination that a change in the post-fix MSE for the in-phase component estimation error between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold
Figure imgf000052_0001
^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ெௌா,ூ)). [00189] In some embodiments, a kth iteration of the second iterative process further includes, in accordance with the determination that the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold (e.g., the second iterative process has not converged), performing a next iteration of the second iterative process (see, e.g., Figure 12, steps 1220 and 1222). For example, in some embodiments, the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold when
Figure imgf000052_0002
or
Figure imgf000052_0003
^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ெௌா,ூ . [00190] In some embodiments, a kth iteration of the second iterative process further includes: in accordance with the determination that (i) the number of iterations ^^ does not satisfy the fifth threshold value (e.g., k does not satisfy the fifth threshold value when k > Threshold4 or k ≥ Threshold4, as illustrated in Figure 12, step 1210), terminating the second iterative process without producing a solution for ^^ூ, ^^^, and ^^ (see, e.g., Figure 12, step 1231). [00191] In some embodiments, in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is not within the corresponding first valid range or (ii) the updated estimate ^ ^ ^^ is not within the corresponding second valid range (i.e., the second iterative process has diverged), the GNSS receiver adjusts the initial estimate for the time delay ^^ from ^ ^^^ to ^ ^ ^^′ (see, e.g., Figure 12, step 1226). The GNSS receiver executes the second iterative process using the initial estimates ^ ^ ^ூ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^ ^^′ as the starting values. [00192] In some embodiments, different initial estimates ^ ^^^ can be tried using the pattern illustrated in Figure 14. For example, if the first estimate 1402 ( ^ ^ ^^) results in a divergence, a second new estimate 1404 ( ^ ^ ^^′ ) can be adjusted to a value that is larger than (e.g., to the right of) that of the first estimate 1402 ( ^^^^). If the second new estimate 1404 also fails, then a third new estimate 1406 ( ^ ^ ^^ ᇱᇱ ) can be adjusted to a value that is smaller than (e.g., to left of) that of the first estimate 1402 ( ^^^^). If the third new estimate 1406 also fails, then a fourth new estimate 1408
Figure imgf000053_0001
can be adjusted to a value that is larger than (e.g., to the right of) that of the second estimate 1404 ( ^ ^ ^^′ ). In some embodiments, the adjustment step size can be a constant (e.g., equal to the correlator spacing). [00193] In some embodiments, executing the second iterative process using the initial estimates ^ ^ ^ூ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^^^′ as the starting values includes: in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ (e.g., produced during a kth iteration of the second iterative process) is not within the corresponding first valid range or (ii) the updated estimate ^ ^ ^^ is not within the corresponding second valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, the GNSS receiver repeats the steps of adjusting the initial estimate for the time delay ^^ and executing the second iterative process (see, e.g., the path corresponding to steps 1226, 1228, and 1204 in Figure 12), and otherwise (e.g., in accordance with a determination that a next adjusted value for the time delay would not fall within a range of time values for the n pairs of correlation samples), the GNSS receiver outputs information that no valid solution for ^^ூ, ^^^, and ^^ has been determined (see, e.g., Figure 12, steps 1226, 1230, and 1231). In some embodiments, the adjustment range (e.g., ^^^^^ 1420 to ^^^^௫1422 in Figure 14) for the time delay ^^ should cover all correlation sample times (e.g., from ^^^ to ^^^) with some margin. Once an adjustment moves beyond this range, it means that all possible initial estimates in this corresponding direction have been tried. If adjustments have covered the whole range without producing an updated estimate ^ ^ ^^,^ within the corresponding first valid range and an updated estimate ^ ^ ^^ within the corresponding second valid range, the GNSS receiver stops the second iterative process and outputs information that there is a no valid solution (e.g., see, e.g., Figure 12, steps 1230, 1231). [00194] In some embodiments, after outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively, and in accordance with the determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: the GNSS receiver assigns a value of zero as the solution for ^^ொ (see, e.g., Figure 12, steps 1218, 1232, 1234, and 1236). [00195] In some embodiments, the GNSS receiver determines the carrier phase multipath error ^^ఌ based on the solutions for ^^ூ and ^^ொ (e.g., using Equation 45: ^^ఌ = tanି^ ^ೂ ^^ ). The GNSS receiver corrects the pseudorandom phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [00196] In some embodiments, after outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively (see, e.g., Figure 12, step 1218), and in accordance with the determination that the phase ^^ଷ^ௗ is not 0° or 180° within the second predefined margin (see, e.g., Figure 12, step 1238): the GNSS receiver starts with the initial estimate for ^^^ ொ,^ (e.g., see above discussion of step 1102, Figure 11) and obtains a solution for ^^ொ via a third iterative process. In the third iterative process, k is a positive integer (e.g., starting from one (e.g., reset to 1 at step 1240 of the flowchart in Figure 12) representing a number of iterations of the third iterative process. A kth iteration of the third iterative process includes: obtaining n estimated quadrature component values ( ^^^( ^^^ି^) values) for the Q(t) samples based on an initial estimate ^^^ ொ,^ି^ (using Equation 39) from a previous (k-1)th iteration of the third iterative process (see, e.g., Figure 12, step 1240), or, in the case of the first iteration (k=1), using the aforementioned initial estimate
Figure imgf000054_0001
; determining, for each estimated quadrature component value
Figure imgf000054_0002
of the ^ ^ ^ ( ^^^ି^ ) values, a corresponding quadrature component estimation error
Figure imgf000054_0003
between (i) the respective estimated quadrature component value ^^^൫ ^^^,^ି^൯ and (ii) the respective quadrature component sample Q(ti) (e.g., using equation 31), thereby obtaining n quadrature component estimation errors ( ^^ொ ( ^^^ ) estimation errors) (see, e.g., Figure 12, step 1240); solving a third set of matrix equations (Q-matrix equation 41) based on the ^^ொ( ^^^) estimation errors, thereby obtaining an updated estimate ^^^ ொ,^ for the kth iteration (using equation 35); in accordance with a determination that the number of iterations ^^ satisfies a sixth threshold value (e.g., k < Threshold5 or k ≤ Threshold5, as illustrated in Figure 12, step 1244): in accordance with a determination that a change in value of one or more predetermined second parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold, determining that the third iterative process has converged (see. e.g., Figure 12, step 1246); outputting the updated estimate ^^ ^ ொ,^ as the solution for ^^ொ (see, e.g., Figure 12 step 1248); and terminating the third iterative process. [00197] For example, in some embodiments, the one or more predetermined second parameters include the parameter ^^ொ. In some embodiments, the determination that the third iterative process has converged is made in accordance with a determination that the change in an estimated value of ^^ொ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g.,
Figure imgf000055_0001
− ^^ ^ ொ,^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^^ொ). ). In some embodiments, the one or more predetermined second parameters include a post-fix mean squared error (MSE) for the quadrature component estimation error, defined by ^^ ^^ ^^ொ,^
Figure imgf000055_0002
^^^,^ା^൯൧ଶ for the kth iteration (See Equation 44). In some embodiments, the determining that the third iterative process has converged is in accordance with a determination that a change in the post-fix MSE for the quadrature component estimation error between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (i.e.,
Figure imgf000055_0003
− ^^ ^^ ^^ொ,^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ெௌா,ொ), and/or if ^^ ^^ ^^ொ,^ itself is smaller than (or does not exceed) a corresponding threshold. [00198] In some embodiments, the kth iteration of the third iterative process further includes, in accordance with a determination that the change in value of the one or more predetermined second parameters between the kth iteration and the (k-1)th iteration does not satisfy a respective corresponding threshold, performing a next iteration of the third iterative process (see, e.g., Figure 12, step 1252). [00199] For example, in some embodiments, the change in value of the one or more predetermined second parameters between the kth iteration and the (k-1)th iteration does not satisfy a respective corresponding threshold when
Figure imgf000055_0004
or when ห ^^ ^ ொ,ఫ − ^^ ^ ொ,^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^^ொ), or when ห ^^ ^^ ^^ொ,^ − ^^ ^^ ^^ொ,^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ெௌா,ொ), or when
Figure imgf000055_0005
[00200] In some embodiments, the kth iteration of the third iterative process further includes: in accordance with the determination that the number of iterations ^^ does not satisfy the sixth threshold value (e.g., k > Threshold5 or k ≥ Threshold5 in Figure 12), terminating the third iterative process without producing a solution for ^^ொ (see, e.g., Figure 12, steps 1244, 1246, and 1248). [00201] In some embodiments, the GNSS receiver computes the carrier phase multipath error ^^ఌ based on the solutions for ^^ூ and ^^ொ (using Equation 45: ^^ఌ = tanି^ ^ೂ ^^ ). The GNSS receiver corrects the pseudorange measurement for the respective satellite in accordance with ∆ ^^ and/or corrects the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. [00202] In some embodiments, in accordance with the determination that the phase ^^ଷ^ௗ corresponds to ±90° within the first predefined margin (see, e.g., Figure 11, step 1108): the GNSS receiver uses the initial estimates
Figure imgf000056_0001
as starting values and obtains a solution for ^^ொ, ^^^, and ^^ via a fourth iterative process (see, using workflow 1300 in Figure 13). k is a positive integer (e.g., starting from one when the fourth iterative process begins) representing a number of iterations of the fourth iterative process. A kth iteration of the fourth iterative process includes: obtaining n estimated quadrature component values ( ^ ^ ^ ( ^^^ି^ ) values) for the Q(t) samples based on initial estimates
Figure imgf000056_0002
and ^ ^ ^^ି^ from a previous (k-1)th iteration of the fourth iterative process (e.g., using equation 29) (see, e.g., Figure 13, step 1304), or for a first iteration of the fourth iterative process, based on previously determine initial estimates (e.g., as explained above with reference to step 1102, Figure 11); determining, for each estimated quadrature component value ^^^൫ ^^^,^ି^൯ of the ^ ^ ^ ( ^^^ି^ ) values, a respective quadrature component estimation error ^^ொ൫ ^^^,^൯ between (i) the respective estimated quadrature component value ^^^൫ ^^^,^ି^൯ and (ii) the respective quadrature component sample Q(ti) (e.g., using equation 31) (see, e.g., Figure 13, step 1304), thereby obtaining n quadrature component estimation errors ( ^^ொ( ^^^) estimation errors); solving a fourth set of matrix equations (e.g., using Least Squares method) based on the ^^ொ ( ^^^ ) estimation errors, thereby obtaining updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ for the k th iteration (see, e.g., Figure 13, step 1308). In accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is within a corresponding third valid range (e.g., the updated estimate ^ ^ ^^,^ is greater than 0) and (ii) the updated estimate ^^^^ is within a corresponding fourth valid range (e.g., ^ ^ ^^ is within the allowably maximal time delay corresponding to the last pair of correlators (e.g., ^ ^ ^^ is greater than 0 but less than ^^^)): in accordance with a determination that the number of iterations ^^ satisfies a seventh threshold value (e.g., k < Threshold6 or k ≤ Threshold6, as illustrated in Figure 13, steps 1310, 1312): in accordance with a determination that a change in value of one or more predetermined second parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the fourth iterative process has converged (see, e.g., Figure 13, step 1316); outputting the updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ொ, ^^^, and ^^, respectively (see, e.g., Figure 13, step 1318); and terminating the fourth iterative process. [00203] For example, in some embodiments there are a few approaches for solving the fourth set of matrix equations. In one approach, equations 32 and 33 are solved separately, while estimates of ^^ொ, ^^^ and ^^ for equation 33 are updated based on equations 35, 36 and 37. In another approach, the two ^^ × 3 matrix equations (equations 32 and 33) are combined into one (2 ^^) × 4 matrix equation with a total of four unknowns (e.g., unknowns ∆^^, ∆^ೂ, ∆^^, and ∆ఋ), and the (2 ^^) × 4 matrix equation is solved (e.g., instead of combining solutions). In yet another approach, equation 32 is solved first, and then equation 41 is solved. In the latter approach, there is no need to solve equation 33 because it has been reduced to equation 41. [00204] For example, in some embodiments, the one or more predetermined second parameters include one or more of ^^ொ, ^^^, and ^^. In some embodiments, the determination that the fourth iterative process has converged is made in accordance with a determination that (i) the change in an estimated value of ^^ொ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g.,
Figure imgf000057_0001
and/or (ii) the change in an estimated value of ^^^ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g.,
Figure imgf000057_0002
ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^^ெ), and/or (iii) the change in an estimated value of ^^ between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., ห ^ ^ ^^ − ^ ^ ^^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ఋ). In some embodiments, the one or more predetermined second parameters include a post-fix mean squared error (MSE) for the quadrature component estimation error, defined by ^^ ^^ ^^ ^ ொ,^ = th ^
Figure imgf000057_0003
the k iteration (See Equation 4a). In some embodiments, determining that the fourth iterative process has converged is in accordance with a determination that a change in the post-fix MSE for the quadrature component estimation error between the kth iteration and the (k-1)th iteration has a magnitude that is smaller than (or does not exceed) a corresponding threshold (e.g., ห ^^ ^^ ^^ொ,^ −
Figure imgf000057_0004
iteration of the fourth iterative process further includes: in accordance with the determination that the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold (e.g., the fourth iterative process has not converged, see Figure 13, Step 1328), performing a next iteration of the fourth iterative process (Figure 13, Step 1330). [00206] For example, in some embodiments, the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold when ห
Figure imgf000058_0001
− ^^ ^ ொ,^ି^ ห ^ ^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^^ூ), or
Figure imgf000058_0002
^^ℎ ^^ ^^ ^^ℎ ^^ ^^ ^^ெௌா,ொ). [00207] In some embodiments, the kth iteration of the fourth iterative process further includes: in accordance with the determination that the number of iterations ^^ does not satisfy the seventh threshold value (e.g., k ≥ Threshold6 or k > Threshold6 in Figure 13, step 1310) (e.g., the fourth iterative process has diverged or failed to converge), terminating the fourth iterative process without producing a solution for ^^ொ (see, e.g., Figure 13, step 1326). [00208] In some embodiments, in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is not within a corresponding third valid range or (ii) the updated estimate ^ ^ ^^ is not within a corresponding fourth valid range: the GNSS receiver updates the initial estimate for the time delay ^^ from ^ ^ ^^ to ^ ^ ^^′ . The GNSS receiver executes the fourth iterative process (i.e., executes the fourth iterative process an additional time) using the initial estimates ^^^ ொ,^ and ^ ^ ^^,^ and the updated initial estimate ^ ^ ^^′ as the starting values. [00209] In some embodiments, the GNSS receiver uses the initial estimates ^^ ^ ொ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^ ^^′ as the starting values. In accordance with a determination that (i) the updated estimate ^^^^,^ is not within the corresponding third valid range or (ii) the updated estimate ^^^^ is not within the corresponding fourth valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, the GNSS receiver repeats the steps of adjusting the initial estimate for the time delay ^^ and executing the fourth iterative process with the updated initial estimate for the time delay ^^. Otherwise, the GNSS receiver outputs information that no valid solution for ^^ொ, ^^^, and ^^ has been determined. [00210] In some embodiments, in conjunction with outputting the updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ொ, ^^^, and ^^, respectively, the GNSS receiver obtains
Figure imgf000059_0001
Because Procedure 1 (see workflow 900, Figure 9) determines that ^^ଷ^ௗ is equal to ±90° (see steps 906, 908, 910, 912, Figure 12, and step 1302, Figure 13) using ^^( ^^^) measurements based on the same equation, equation 22, use to solve for ^^ூ, in Procedure 2 (see, e.g., workflow 1100 in Figure 11 and workflow 1300 in Figure 13), there is no need to repeat this process of solving for some unknowns from ^^( ^^^) measurements. Instead, the solution ^^ூ to equation 22 obtained in Procedure 1 can be assigned not only to the initial estimate ^ ^ ^ூ,^ but also to the final solution of ^^ூ. [00211] In some embodiments, the GNSS receiver computes the carrier phase multipath error ^^ఌ based on the solution for ^^ூ and the solution for ^^^ (e.g., using Equation 46). The GNSS receiver corrects the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. At this point, ^^ఌ and ^^ଷ^ௗ are known. The phase difference ^^^ can be computing by: ^^^ = 180° - ( ^^ఌ + ^^ଷ^ௗ). [00212] It will be understood that although the terms “first,” “second,” etc. may be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another. For example, a first contact could be termed a second contact, and, similarly, a second contact could be termed a first contact, without changing the meaning of the description, so long as all occurrences of the “first contact” are renamed consistently and all occurrences of the second contact are renamed consistently. The first contact and the second contact are both contacts, but they are not the same contact. [00213] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the claims. As used in the description of the embodiments and the appended claims, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will also be understood that the term "and/or" as used herein refers to and encompasses any and all possible combinations of one or more of the associated listed items. It will be further understood that the terms "comprises" and/or "comprising," when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. [00214] As used herein, the term “if” may be construed to mean “when” or “upon” or “in response to determining” or “in accordance with a determination” or “in response to detecting,” that a stated condition precedent is true, depending on the context. Similarly, the phrase “if it is determined [that a stated condition precedent is true]” or “if [a stated condition precedent is true]” or “when [a stated condition precedent is true]” may be construed to mean “upon determining” or “in response to determining” or “in accordance with a determination” or “upon detecting” or “in response to detecting” that the stated condition precedent is true, depending on the context. [00215] The foregoing description, for purpose of explanation, has been described with reference to specific embodiments. However, the illustrative discussions above are not intended to be exhaustive or to limit the invention to the precise forms disclosed. Many modifications and variations are possible in view of the above teachings. The embodiments were chosen and described in order to best explain the principles of the invention and its practical applications, to thereby enable others skilled in the art to best utilize the invention and various embodiments with various modifications as are suited to the particular use contemplated.

Claims

What is claimed is: 1. A method of mitigating the effect of a multipath-induced error in a global navigation satellite system (GNSS), performed at a respective GNSS signal receiver, comprising: receiving a band-limited composite signal corresponding to a respective satellite in the GNSS, including a band-limited direct-path signal and a band-limited multipath signal, wherein the direct-path signal and the multipath signal are modulated with a pseudorandom (PN) code; obtaining, for a code chip edge transition of the PN code, n pairs of correlation samples of the composite signal during a corresponding n time instances, each having a different time offset from the code chip edge transition, wherein: the code chip edge transition has a predetermined filter step response function SR(t); each respective pair of correlation samples of the n pairs of correlation samples consists of a respective in-phase component sample I(ti) and a respective quadrature component sample Q(ti), wherein i is a positive integer from one to n; and the n pairs of correlation samples consist of n in-phase component samples of the composite signal (I(t) samples) and n quadrature component samples of the composite signal (Q(t) samples); obtaining, for each respective pair of correlation samples of the n pairs of correlation samples, a corresponding sample of the filter step response function SR(ti) corresponding in time with the respective pair of correlation samples, thereby obtaining n samples of the filter step response function (predetermined SR(t) samples) during the corresponding n time instances; in accordance with a determination that a computed similarity between the I(t) samples and the SR(t) samples satisfies a first threshold value, determining that a phase ^^ଷ^ௗ is ±90° within a first predefined margin, wherein the phase ^^ଷ^ௗ is equal to 180° minus a sum of (i) a carrier phase multipath error ^^ఌ and (ii) a phase difference ^^^ between the direct- path signal and the multipath signal; in accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples does not satisfy the first threshold value: solving a first set of matrix equations, thereby obtaining a solution for tan ^^ଷ^ௗ and for ∆ ^^ tan ^^ଷ^ௗ, wherein ∆ ^^ is a response time error of the code chip edge transition due to the multipath signal; in accordance with a determination that | ^^ ^^ ^^ ^^ଷ^ௗ | satisfies a second threshold value, determining that the phase ^^ଷ^ௗ is 0° or 180° within a second predefined margin; and in accordance with a determination that | ^^ ^^ ^^ ^^ଷ^ௗ| does not satisfy the second threshold value, determining ^^ଷ^ௗ in accordance with the solution for tan ^^ଷ^ௗ; adjusting a pseudorange measurement for the respective satellite in accordance with the determined ∆ ^^; and/or adjusting a carrier phase measurement for the respective satellite in accordance with a parameter corresponding to the determined ^^ଷ^ௗ; and performing a navigation function using the adjusted pseudorange and/or the adjusted carrier phase measurement for the respective satellite. 2. The method of claim 1, further comprising: determining a sign of the carrier phase multipath error ^^ఌ; in accordance with a determination that the computed similarity between the I(t) samples and the SR(t) samples satisfies the first threshold value: determining that ^^ଷ^ௗ is +90° within the first predefined margin when the sign of the carrier phase multipath error ^^ఌ is positive; and determining that ^^ଷ^ௗ is -90° within the first predefined margin when the sign of the carrier phase multipath error ^^ఌ is negative. 3. The method of claim 1 or claim 2, further comprising: in accordance with a determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: determining a position of a peak of the I(t) samples relative to a position of a peak of the predetermined SR(t) samples; in accordance with a determination that the position of the peak of the I(t) samples occurs later in time than the position of the peak of the SR(t) samples, determining that ^^ଷ^ௗ is 0° within the second predefined margin; and in accordance with a determination that the position of the peak of the I(t) samples occurs earlier in time than the position of the peak of the predetermined SR(t) samples, determining that ^^ଷ^ௗ is 180° within the second predefined margin. 4. The method of any of claims 1-3, wherein the computed similarity between the I(t) samples and the SR(t) samples is obtained using a mean squared error (MSE), defined by
Figure imgf000062_0001
wherein ^^ூ = ^^ௗ cos ^^ఌ, ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error. 5. The method of any of claims 1-4, wherein solving the first set of matrix equations comprises: obtaining (i) an initial estimate ^^ ^ ଷ^ௗ,^ for the phase ^^ଷ^ௗ, (ii) an initial estimate ∆ ^ ^^^ for the response time error ∆ ^^, and (iii) an initial estimate ^ ^ ^ௗା^ for a magnitude ^^ௗା^ of the composite signal; and using the initial estimates ^^ ^ ଷ^ௗ,^ , ∆ ^ ^^^, and ^ ^ ^ௗା^ to solve the first set of matrix equations via a Least Squares fitting process to obtain the solution for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ. 6. The method of claim 5, further comprising: in accordance with the determination that | ^^ ^^ ^^ ^^ଷ^ௗ | does not satisfy the second threshold value, obtaining a solution for tan ^^ଷ^ௗ and ∆ ^^ tan ^^ଷ^ௗ via a first iterative process, each iteration of the first iterative process including: obtaining an updated estimate tan ^^^ ଷ^ௗ,^ for the phase ^^ଷ^ௗ and an updated estimate ∆ ^ ^^^ for the response time error ∆ ^^; in accordance with a determination that a number of iterations for the first iterative process satisfies a third threshold value: outputting information that no valid solution for ^^ଷ^ௗ has been determined; and terminating the first iterative process without producing a solution for ^^ଷ^ௗ; in accordance with a determination that a number of iterations for first iterative process does not satisfy the third threshold value: in accordance with a determination that a magnitude based on the updated estimate ∆ ^ ^^^ or based on the updated estimate tan ^^ ^ ଷ^ௗ,^ satisfies a fourth threshold value: determining that the Least Squares fitting process has converged; outputting the updated estimate tan ^^ ^ ଷ^ௗ,^ and the updated estimate ∆ ^ ^^^ as the solutions for tan ^^ଷ^ௗ and ∆ ^^, respectively; calculating the phase ^^ଷ^ௗ from the updated estimate tan
Figure imgf000063_0001
terminating the first iterative process; and in accordance with a determination that (i) the magnitude based on the updated estimate ∆ ^ ^^^ or based on the updated estimate tan ^^ ^ ଷ^ௗ,^ does not satisfy the fourth threshold value, performing a next iteration of the first iterative process. 7. The method of any of claims 1-6, further comprising: in accordance with a determination that a valid solution for the phase ^^ଷ^ௗ has been determined: determining (i) an initial estimate ^^^ூ,^ for ^^ூ, wherein ^^ூ = ^^ௗ cos ^^ఌ, ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error; (ii) an initial estimate ^^ ^ ொ,^ for a magnitude ^^ொ, wherein ^^ொ = ^^ௗ sin ^^ఌ; (iii) an initial estimate ^ ^ ^^,^ for a magnitude ^^^ of the multipath signal, and (iv) an initial estimate ^^^^ for a time delay ^^ of the multipath signal relative to the direct-path signal; using at least a subset of (i) the initial estimate ^ ^ ^ூ,^ , (ii) the initial estimate ^^ ^ ொ,^ , (iii) the initial estimate ^ ^ ^^,^, and (iv) the initial estimate ^ ^^^ as starting values, performing one or more iterative processes to obtain solutions for at least a subset of ^^ூ, ^^ொ, ^^^, and ^^; computing the carrier phase multipath error ^^ఌ in accordance with the obtained solutions for at least the subset of ^^ூ, ^^ொ, ^^^, and ^^; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. 8. The method of any of claims 1-6, further comprising: in accordance with a determination that a valid solution for the phase ^^ଷ^ௗ has been determined: determining (i) an initial estimate ^ ^ ^ூ,^ for ^^ூ, wherein ^^ூ = ^^ௗ cos ^^ఌ, ^^ௗ is a magnitude of the direct-path signal, and ^^ఌ is the carrier phase multipath error; (ii) an initial estimate ^^ ^ ொ,^ for a magnitude ^^ொ, wherein ^^ொ = ^^ௗ sin ^^ఌ; (iii) an initial estimate ^ ^ ^^,^ for a magnitude ^^^ of the multipath signal, and (iv) an initial estimate ^ ^ ^^ for a time delay ^^ of the multipath signal relative to the direct-path signal. 9. The method of claim 8, wherein, in accordance with a determination that the code chip edge transition is a code chip edge-up transition: the initial estimate ^^^ூ,^ corresponds to a magnitude of a last in-phase component sample I(tn) of the I(t) samples; and the initial estimate ^^^ொ,^ corresponds to a magnitude of a last quadrature component Q(tn) of the Q(t) samples. 10. The method of claim 8 or claim 9, wherein, in accordance with the code chip edge transition is a code chip edge-down transition: the initial estimate ^^^ூ,^ corresponds to a magnitude of a first in-phase component sample I(t1) of the I(t) samples; and the initial estimate ^^^ொ,^ corresponds to a magnitude of a first quadrature component Q(t1) of the Q(t) samples. 11. The method of any of claims 8-10, wherein: in accordance with the determination that the phase ^^ଷ^ௗ corresponds to ±90° within the first predefined margin: obtaining a solution for ^^ூ by solving a matrix equation
Figure imgf000065_0001
and assigning the obtained solution as the initial estimate ^ ^ ^ூ,^ . 12. The method of any of claims 8-11, wherein: in accordance with the determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: assigning the initial estimate ^^ ^ ொ,^ as zero. 13. The method of any of claims 8-12, wherein the initial estimate ^ ^ ^^,^ is 0.5. 14. The method of any of claims 8-13, further comprising: in accordance with the determination that the phase ^^ଷ^ௗ does not correspond to ±90° within the first predefined margin: using the initial estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^^^ as starting values to obtain a solution for ^^ூ, ^^^, and ^^ via a second iterative process, wherein k is a positive integer representing a number of iterations of the second iterative process and a kth iteration of the second iterative process includes: obtaining n estimated in-phase component values ( ^ ^ ^ ( ^^^ି^ ) values) for the I(t) samples based on estimates ^^ ^ ூ,^ି^ , ^^ ^ ^,^ି^ , and ^ ^ ^^ି^ from a previous (k-1) th iteration of the second iterative process; determining, for each estimated in-phase component value ^ ^ ^൫ ^^^,^ି^൯ of the ^ ^ ^ ( ^^^ି^ ) values, a corresponding in-phase component estimation error ^^ூ൫ ^^^,^൯ between (i) the respective estimated in-phase component value ^^^൫ ^^^,^ି^൯ and (ii) the respective in-phase component sample I(ti) of the I(t) samples, thereby obtaining n in-phase component estimation errors ( ^^ூ ( ^^^ ) estimation errors); solving a second set of matrix equations based on the ^^ூ ( ^^^ ) estimation errors, thereby obtaining updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ for the k th iteration; in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is within a corresponding first valid range and (ii) the updated estimate ^ ^ ^^ is within a corresponding second valid range: in accordance with a determination that the number of iterations ^^ satisfies a fifth threshold value: in accordance with a determination that a change in value of one or more predetermined first parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the second iterative process has converged; outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively; and terminating the second iterative process; in accordance with the determination that the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold, performing a next iteration of the second iterative process; and in accordance with the determination that (i) the number of iterations ^^ does not satisfy the fifth threshold value, terminating the second iterative process without producing a solution for ^^ூ, ^^^, and ^^. 15. The method of claim 14, further comprising: in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is not within the corresponding first valid range or (ii) the updated estimate ^^^^ is not within the corresponding second valid range: adjusting the initial estimate for the time delay ^^ from ^ ^ ^^ to ^ ^ ^^′ ; and executing the second iterative process using the initial estimates ^ ^ ^ூ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^ ^^′ as the starting values. 16. The method of claim 15, wherein: executing the second iterative process using the initial estimates ^ ^ ^ூ,^
Figure imgf000067_0001
and the adjusted initial estimate ^^^^′ as the starting values includes: in accordance with a determination that (i) the updated estimate ^^^^,^ is not within the corresponding first valid range or (ii) the updated estimate ^^^^ is not within the corresponding second valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, repeating the steps of adjusting the initial estimate for the time delay ^^ and executing the second iterative process; and otherwise, outputting information that no valid solution for ^^ூ, ^^^, and ^^ has been determined. 17. The method of any of claims 14-16, further comprising: after outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively, and in accordance with the determination that the phase ^^ଷ^ௗ is 0° or 180° within the second predefined margin: assigning a value of zero as the solution for ^^ொ. 18. The method of claim 17, further comprising: determining the carrier phase multipath error ^^ఌ based on the solutions for ^^ூ and ^^ொ; and correcting a pseudorandom phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. 19. The method of any of claims 14-18, further comprising: after outputting the updated estimates ^ ^ ^ூ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ூ, ^^^, and ^^, respectively, and in accordance with the determination that the phase ^^ଷ^ௗ is not 0° or 180° within the second predefined margin: starting with the initial estimate for ^^^ ொ,^ , obtaining a solution for ^^ொ via a third iterative process, wherein k is a positive integer representing a number of iterations of the third iterative process and a kth iteration of the third iterative process includes: obtaining n estimated quadrature component values (
Figure imgf000068_0001
values) for the Q(t) samples based on an initial estimate ^^^ fr th ொ,^ି^ om a previous (k-1) iteration of the third iterative process; determining, for each estimated quadrature component value
Figure imgf000068_0002
of the ^ ^ ^ ( ^^^ି^ ) values, a corresponding quadrature component estimation error ^^ொ൫ ^^^,^൯ between (i) the respective estimated quadrature component value
Figure imgf000068_0003
and (ii) the respective quadrature component sample Q(ti), thereby obtaining n quadrature component estimation errors ( ^^ொ ( ^^^ ) estimation errors); solving a third set of matrix equations based on the
Figure imgf000068_0004
estimation errors, thereby obtaining an updated estimate ^^ ^ ொ,^ for the k th iteration; in accordance with a determination that the number of iterations ^^ satisfies a sixth threshold value: in accordance with a determination that a change in value of one or more predetermined second parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the third iterative process has converged; outputting the updated estimate
Figure imgf000068_0005
the solution for ^^ொ; and terminating the third iterative process; in accordance with a determination that the change in value of the one or more predetermined second parameters between the kth iteration and the (k-1)th iteration does not satisfy a respective corresponding threshold, performing a next iteration of the third iterative process; and in accordance with the determination that (i) the number of iterations ^^ does not satisfy the sixth threshold value, terminating the third iterative process without producing a solution for ^^ொ. 20. The method of claim 19, including: computing the carrier phase multipath error ^^ఌ based on the solutions for ^^ூ and ^^ொ; and correcting the pseudorange measurement for the respective satellite in accordance with ∆ ^^; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. 21. The method of any of claims 8-20, further comprising: in accordance with the determination that the phase ^^ଷ^ௗ corresponds to ±90° within the first predefined margin: using the initial estimates
Figure imgf000069_0001
as starting values, obtaining a solution for ^^ொ, ^^^, and ^^ via a fourth iterative process, wherein k is a positive integer representing a number of iterations of the fourth iterative process and a kth iteration of the fourth iterative process includes: obtaining n estimated quadrature component values
Figure imgf000069_0002
values) for the Q(t) samples based on initial estimates ^^^ொ,^ି^ , ^^^ , and ^^^ from a previous (k-1)th ^,^ି^ ^ି^ iteration of the fourth iterative process; determining, for each estimated quadrature component value ^ ^ ^൫ ^^^,^ି^൯ of the ^ ^ ^ ( ^^^ି^ ) values, a respective quadrature component estimation error ^^ொ൫ ^^^,^൯ between (i) the respective estimated quadrature component value ^^^൫ ^^^,^ି^൯ and (ii) the respective quadrature component sample Q(ti), thereby obtaining n quadrature component estimation errors ( ^^ொ ( ^^^ ) estimation errors); solving a fourth set of matrix equations based on the ^^ொ( ^^^) estimation errors, thereby obtaining updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ for the k th iteration; in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is within a corresponding third valid range and (ii) the updated estimate ^ ^ ^^ is within a corresponding fourth valid range: in accordance with a determination that the number of iterations ^^ satisfies a seventh threshold value: in accordance with a determination that a change in value of one or more predetermined second parameters between the kth iteration and the (k-1)th iteration satisfies a respective corresponding threshold: determining that the fourth iterative process has converged; outputting the updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ொ, ^^^, and ^^, respectively; and terminating the fourth iterative process; in accordance with the determination that the change in value of the one or more predetermined parameters between the kth iteration and the (k-1)th iteration does not satisfy the respective corresponding threshold, performing a next iteration of the fourth iterative process; and in accordance with the determination that (i) the number of iterations ^^ does not satisfy the seventh threshold value, terminating the fourth iterative process without producing a solution for ^^ொ. 22. The method of claim 21, further comprising: in accordance with a determination that (i) the updated estimate ^^^^,^ is not within a corresponding third valid range or (ii) the updated estimate ^ ^ ^^ is not within a corresponding fourth valid range: updating the initial estimate for the time delay ^^ from ^ ^ ^^ to ^ ^ ^^′ ; and executing the fourth iterative process using the initial estimates ^^ ^ ொ,^ and ^ ^ ^^,^ and the updated initial estimate ^ ^ ^^′ as the starting values. 23. The method of claim 22, further comprising: using the initial estimates ^^ ^ ொ,^ and ^ ^ ^^,^ and the adjusted initial estimate ^ ^ ^^′ as the starting values: in accordance with a determination that (i) the updated estimate ^ ^ ^^,^ is not within the corresponding third valid range or (ii) the updated estimate ^ ^ ^^ is not within the corresponding fourth valid range: in accordance with a determination that a next adjusted value for the time delay would fall within a range of time values for the n pairs of correlation samples, repeating the steps of updating the initial estimate for the time delay ^^ and executing the fourth iterative process; and otherwise, outputting information that no valid solution for ^^ொ, ^^^, and ^^ has been determined. 24. The method of any of claims 21-23, further comprising: in conjunction with outputting the updated estimates ^^ ^ ொ,^ , ^ ^ ^^,^, and ^ ^ ^^ as the solution for ^^ொ, ^^^, and ^^, respectively: obtaining a solution for ^^ூ by solving a matrix equation
Figure imgf000071_0001
25. The method of claim 24, further comprising: computing the carrier phase multipath error ^^ఌ based on the solution for ^^ூ and the solution for ^^^; and correcting the carrier phase measurement for the respective satellite in accordance with the carrier phase multipath error ^^ఌ. 26. A device for mitigating the effect of a multipath-induced error in a global navigation satellite system (GNSS), comprising: one or more processors; and memory storing instructions that, when executed by the one or more processors, cause the device to perform the method of any of claims 1-25. 27. A computer-readable storage medium having stored thereon program code instructions that, when executed one or more processors, cause the one or more processors to perform the method of any of claims 1-25.
PCT/US2023/034191 2023-06-12 2023-09-29 Mitigation of multipath-induced error effects in a navigation satellite signal receiver Pending WO2024258418A1 (en)

Priority Applications (3)

Application Number Priority Date Filing Date Title
AU2023456885A AU2023456885A1 (en) 2023-06-12 2023-09-29 Mitigation of multipath-induced error effects in a navigation satellite signal receiver
CN202380099371.7A CN121488173A (en) 2023-06-12 2023-09-29 Suppression of multipath induced error effects in navigation satellite signal receivers
EP23895569.4A EP4695629A1 (en) 2023-06-12 2023-09-29 Mitigation of multipath-induced error effects in a navigation satellite signal receiver

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
US202363472574P 2023-06-12 2023-06-12
US63/472,574 2023-06-12
US18/373,913 US20240411029A1 (en) 2023-06-12 2023-09-27 Phase Multipath Mitigation in a Navigation Satellite Signal Receiver
US18/373,913 2023-09-27

Publications (1)

Publication Number Publication Date
WO2024258418A1 true WO2024258418A1 (en) 2024-12-19

Family

ID=91302290

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US2023/034191 Pending WO2024258418A1 (en) 2023-06-12 2023-09-29 Mitigation of multipath-induced error effects in a navigation satellite signal receiver

Country Status (1)

Country Link
WO (1) WO2024258418A1 (en)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060083293A1 (en) * 2004-10-18 2006-04-20 Keegan Richard G Phase multi-path mitigation
US20140376598A1 (en) * 2013-05-03 2014-12-25 Deere & Company Phase Multi-Path Mitigation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060083293A1 (en) * 2004-10-18 2006-04-20 Keegan Richard G Phase multi-path mitigation
US20140376598A1 (en) * 2013-05-03 2014-12-25 Deere & Company Phase Multi-Path Mitigation

Similar Documents

Publication Publication Date Title
US20220276389A1 (en) Satellite navigation receiver for acquisition of gnss signals
JP5957025B2 (en) User receiver positioning method
US10859709B2 (en) Satellite navigation receiver with fixed point sigma rho filter
US9515697B2 (en) Scanning correlator for global navigation satellite system signal tracking
USRE42676E1 (en) Method for determining a boundary of an information element, a system, and an electronic device
US6725157B1 (en) Indoor GPS clock
US12379508B2 (en) Satellite navigation receiver with aggregate channel digital baseband processing
US8306093B2 (en) Method and apparatus for multipath mitigation
US20060140254A1 (en) Multi-path detection method for CDMA receivers
US12066551B2 (en) Multi-constellation, multi-frequency GNSS system for interference mitigation
US11391847B2 (en) GNSS correlation distortion detection and mitigation
US7903026B2 (en) Positioning apparatus and positioning apparatus control method
US12153143B2 (en) Multi-constellation, multi-frequency GNSS system for interference mitigation
US6252546B1 (en) Method and apparatus for processing multipath reflection effects in timing systems
EP4542262A1 (en) Processing a gnss signal and wiping off a data bit based on doppler estimates
US7693211B2 (en) Fast fourier transform based phase locked loop for navigational receivers
EP3362818B1 (en) Satellite navigation receiver with fixed point sigma rho filter
US8295411B2 (en) Method and system for maintaining integrity of a binary offset carrier signal
WO2009091516A1 (en) Measurement of energy potential (signal-to-noise ratio) in digital global navigation satellite systems receivers
US20250123406A1 (en) Processing a gnss signal based on doppler estimates
US20250123407A1 (en) Processing a gnss signal based on doppler estimates
US20240411029A1 (en) Phase Multipath Mitigation in a Navigation Satellite Signal Receiver
US6618006B1 (en) Method for defining the error of reference time and an electronic device
WO2024258418A1 (en) Mitigation of multipath-induced error effects in a navigation satellite signal receiver

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 23895569

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: AU2023456885

Country of ref document: AU

Ref document number: 2023895569

Country of ref document: EP

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023456885

Country of ref document: AU

Date of ref document: 20230929

Kind code of ref document: A

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

REG Reference to national code

Ref country code: BR

Ref legal event code: B01A

Ref document number: 112025026063

Country of ref document: BR

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

ENP Entry into the national phase

Ref document number: 2023895569

Country of ref document: EP

Effective date: 20251112

NENP Non-entry into the national phase

Ref country code: DE

WWP Wipo information: published in national office

Ref document number: 2023895569

Country of ref document: EP