System for improving the sensitivity and stability of optical polarimetric measurements

[0002] This application claims the benefit of prior filed co-pending U. S. Provisional Application Serial No. 60/145,909 filed Jul. 27, 1999 entitled “SYSTEM FOR FUNDAMENTALLY IMPROVING THE SENSITIVITY AND STABILITY OF OPTICAL POLARIMETRIC MEASUREMENTS, WITH SPECIFIC APPLICATION” by James K. Kroeger and Alan J. Leszinske.

[0003] The invention relates to polarimetry, especially as applied to noninvasive measuring of blood glucose concentration in diabetics. It is known that this phenomenon offers the potential for developing a noninvasive blood glucose analyzer.

[0004] Diabetes is a disease which entails a large number of associated complications. Retinal deterioration leading to blindness and impaired circulation leading to limb amputation, kidney failure and heart disease are just some of the more serious complications. Many of these complications result from the large excursions in blood glucose concentrations common to diabetics due to dietary intake, inadequate exercise, genetic predisposition, and complicated by infrequent and inaccurate monitoring of the blood glucose levels. Current methods of in home monitoring of blood glucose involve the lancing or sticking of a finger and external measurement of the glucose content of the blood sample and or urine sampling by use of a litmus strip test comparing a color change relative to glucose concentrations.

[0005] Although many diabetes patients should use the “finger sticking” test to obtain blood for glucose concentration measurements four or more times per day, studies show that very few patients do this unless they absolutely have to, and many patients only do it a few times at the beginning of their treatment until they establish what they think is a pattern in their required medication schedule. They then stop the regular and frequent finger sticking tests and simply take their insulin injections or oral medications on the assumption that their body chemistry is thereafter constant. This leads to large changes in glucose concentration in the patient's blood, which in turn leads to a variety of serious medical consequences to the patient. For example, it is estimated that in 1996 there were over fifty thousand amputations of limbs due to complications of diabetes in the U.S.

[0006] Diabetics recover from cuts and bruises more slowly than do nondiabetics. This very real and basic discomfort also causes many diabetics to minimize the frequency of or altogether ignore blood glucose testing, resulting in a higher frequency of complications than otherwise would be the case. A small accurate device that could make blood glucose measurements on a non-invasive basis would be of great value to the diabetic in that it would greatly encourage frequent monitoring of blood glucose levels without pain.

[0007] It is well known that glucose in solution is an optically active material. That is, it will cause the plane of polarization of light traversing the solution to be rotated. The quantitative relationship between the amount of polarization rotation, the glucose concentration, and the optical path length of the solution has been clearly established. This is expressed mathematically as: $\left\{\alpha \right\}=\frac{100\alpha }{C*L}$ $\mathrm{or}\text{:}$ $C=\frac{100\alpha }{\left\{\alpha \right\}*L}$

[0008] Where:

[0009] α is the polarization rotation in degrees;

[0010] {α}is the specific rotation constant of glucose; ({α}=45.1 degrees per decimeter (dm) per gram per milliliter for glucose at a wavelength of 633 nanometers);

[0011] L is the path length in the solution in dm, (where 1 dm=10 centimeters (cm);

[0012] C is the glucose concentration in grams (g) per 100 milliliter of solution or g/dL. (From “Sugar Analysis”, 3rd Edition, Browne & Zerban, John Wiley & Sons, 1941, page 263.)

[0013] For the clinically meaningful glucose concentration range from 25 to 500 mg/dL (milligrams per deciliter) and a path length of 1 cm, the observed rotation ranges from about 0.00113 degrees to 0.02255 degrees at a wavelength of 633 nanometers.

[0014] It is known that human tissue has an absorption minima in the wavelength range from about 750 nanometers to 900 nanometers. Because there are no fundamental absorption processes in this region, human tissue has a reasonable optical transmission in this region of the spectrum. Light scattering by tissue remains a problem, which may limit the path length to less than 4 mm, dependent upon the type of tissue.

[0015] All of the prior art systems using crossed polarizers use only a single frequency, usually in conjunction with a null control system and a lock in amplifier that operates only at that single frequency. The prior art null compensation techniques all involve inserting a sample between the first and second polarizers and driving a Faraday modulator to reestablish the extinction condition. The problem with the prior techniques of establishing a null condition at extinction in a system using crossed polarizers is that the laser, optical modulator, and other components have parameters which drift from the time that the null condition or extinction is initially established and the time at which the sample to be measured is placed between the polarizers and an extinction condition is reestablished to determine the phase rotation caused by the sample.

[0016] According to the article “Non-Invasive Optical Glucose Sensing—An Overview” by Gerard L. Coté, PhD. Journal of Clinical Engineering, Vol. 22, No. 4, July/August 1997, a path length of 4 mm through human soft tissue (other than the eye) attenuates or scatters 95% of the signal. We conducted tests to confirm the general claims by Coté and found that both scattering and absorption are strongly wavelength dependent.

[0017] Because of the impracticality of using prior art devices and techniques to accurately measure such a small signal, the prior art use of polarimetry to measure glucose concentration levels in human tissue has been based primarily on passing light through the transparent tissue of the anterior chamber of the human eye.

[0018] The prior art fails to provide any practical, workable polarimeter system which can consistently provide accurate measurements of the glucose level in human tissue because of the inadequate sensitivity and the large degree of instability of the prior art devices. There is a strong but unmet need for a practical, reliable system which overcomes the problems of the prior art to provide a practical, reasonably priced, noninvasive system for measurement of human glucose levels.

[0019] Accordingly, it is an object of the invention to provide a device capable of consistently and accurately measuring the concentration of an optically active ingredient in a sample.

[0020] It is another object of the invention to provide a practical, economical device for noninvasive measurement of glucose levels in diabetics.

[0021] It is another object of the invention to avoid instrument instability problems which have in part prevented success of prior attempts to provide a practical system using polarimetry to noninvasively measure blood glucose levels in diabetics.

[0022] It is another object of the invention to provide a device capable of measuring an optically sensitive ingredient in biological tissue in a noninvasive manner more accurately than has been achieved in the prior art.

[0023] It is another object of this invention to provide a new very sensitive and very stable polarization spectrometer which has applications in certain types of chemical analysis.

[0024] It is another object of the invention to provide a device capable of measuring optical rotation in the presence of large percentages of more than about 95% scattered light.

[0025] It is another object of the invention to provide an improved polarimeter which is more sensitive and more stable than prior art polarimeters.

[0026] Briefly described, and in accordance with one embodiment thereof, the invention provides a system for polarimetric measurement of the concentration of a substance, such as glucose, in a sample, including a laser beam passing through a first polarizer and an optical modulator and then split into (1) a measurement beam which is analyzed and directed to a first detector coupled to a first amplifier, and (2) a reference beam which is analyzed and directed to a second detector coupled to a second amplifier. Identical multiple filtering and summing operations are performed on outputs of the first and second amplifiers to produce a first Ψ²/2 signal and a first 2βΨ signal in response to the measurement beam and a second Ψ²/2 signal and a second 2βΨ signal in response to the reference beam. The measurement beam is stabilized by a first control loop that compares the second Ψ²/2 signal to a first reference signal to produce a first error signal and a second control loop that compares the second 2βΨ signal to a second reference signal to produce a second error signal. The first error signal is multiplied by a modulation signal to produce a modulation feedback signal and adding it to the second error signal to produce a combined modulation and zeroing feedback signal. The optical modulator then is driven in response to the combined modulation and zeroing feedback signal to minimize the first and second error signals. A first value of β is computed from the first Ψ²/2 signal and the first 2βΨ signal with no sample in the path of the measurement beam, and a second value of β is computed from the first Ψ²/2 signal and the first 2βΨ signal with the sample in the path of the measurement beam. The difference between the first and second values of β is converted to a value of concentration of the optically active substance in the sample by reference to a look-up table or algorithm. Both a primarily hardware implementation of the invention and a primarily sofware/firmware DSP implementation of the invention are disclosed.

[0027] FIG. 1 is an optical schematic of a basic prior art polarimeter system.

[0028] FIG. 2 shows a modulating sinusoidal signal applied to the sin²δ curve of 90° crossed polarizers; also shown is the resulting doubled frequency output signal.

[0029] FIG. 3 is a graph comparing the sensitivity of a single frequency lock in amplifier, cos (ωt), the derivative dT/dδ, and Algorithm 1.

[0030] FIG. 4 is a graph showing the distribution of light energy detected at the output of the polarimeter between all even frequencies, the DC term, along with the total energy, sin²δ, for a range of optical modulation.

[0031] FIG. 5 is a graph depicting the distribution of light energy detected at the output of the polarimeter for even frequencies and Algorithm 2, (Ψ²/2), and the total energy, sin²δ.

[0032] FIG. 6 shows the shape of the waveform of light energy detected at the output of the polarimeter at 75° modulation.

[0033] FIG. 7 is similar to FIG. 2 with the addition of sample rotation signal β, and it also illustrates graphically the optical instabilities that are corrected by feedback loop stabilization methods.

[0034] FIG. 8 is a graph useful in comparing the sensitivity of the stabilization loop of the present invention to that of the prior art.

[0035] FIG. 9 is a graph of the known relationship of specific rotation of glucose to wavelength.

[0036] FIG. 10A is a simplified optics diagram of a presently preferred polarimeter system of the invention.

[0037] FIG. 10B is a detailed optics diagram of one preferred embodiment of the polarimeter system of FIG. 11A.

[0038] FIG. 10C is a side view diagram of a portion of the system of FIG. 10B.

[0039] FIG. 11A is a block diagram of an analog implementation of the present invention.

[0040] FIG. 11B is a block diagram of digital circuitry for producing the digital signals required as inputs to the system shown in FIG. 11A.

[0041] FIG. 12 is a block diagram of a frequency domain implementation using a digital signal processor (DSP) of the present invention.

[0042] FIG. 13 is a flow chart of the operation of the digital signal processor in the diagram of FIG. 12.

[0043] FIG. 14 is an optical diagram of an alternative embodiment to that shown in FIG. 10B.

[0044] FIG. 15 is a perspective drawing of an implementation of FIG. 10B adapted to receive an ear lobe as the sample.

[0045] FIGS. 16A-C illustrates the human eye as a sample in any of FIGS. 10A, 10B and 14 for different paths of the entrance beam as it propagates through the acqueous humor.

[0046] The components of a typical prior art polarimeter, shown in FIG. 1, include a laser 2 which is a monochromatic light source, the output of which is transmitted to a polarizer 4. The light passing from polarizer 4 is called polarized light. The polarized beam passes through an optical modulator 20 which can be used to rotate the polarized light in a totally electronic manner, as opposed to the conventional method of mechanically rotating polarizer 4 or an analyzer 7 to effect the measurement. The laser beam then passes through an analyzer 7. The second polarizer of a pair of crossed polarizers is called an “analyzer”. Analyzer 7 is said to “analyze” the polarized light from the optical modulator 20. The light emanating from analyzer 7 impinges on a suitable detector 8, such as a silicon PIN photodiode, which produces an output signal 11 that represents the amount of light that is transmitted through the entire polarimeter.

[0047] As the wavelength of light is increased, the specific rotation {α}decreases, to a value of 21.7 at a wavelength of 900 nanometers (see FIG. 9). At 900 nm and a path length of 4 mm, the observed rotation is 0.00022° at a concentration of 25 mg/dL. We have discovered that a useable system must have a basic sensitivity of the order of about 0.00001 degrees, or 10 micro degrees. To be useful, the system must also be stable to 0.00001 degrees. Using the stabilization methods described herein, we can achieve this required stability.

[0048] Prior art polarimeter 1 is capable of analyzing polarized light. Its precision is limited by the quality of polarizer 4 and analyzer 7 and the stability of all of the components, especially laser diode 2, and the capability of determining the “position” of the intensity minimum at the detector.

[0049] The following mathematical analysis will show the development of new algorithms which can be used to fundamentally improve both the sensitivity and stability of optical polarimetric measurements.

[0050] The transmission of light through two polarizers is given by:

[0051] I/I₀=K₀cos²θ+K₉₀sin²θ  Eq.(1)

[0052] where

[0053] I is the light transmitted through the polarizers

[0054] I₀is the light incident upon the polarizers

[0055] θ is the relative angle between the two polarizers,

[0056] K₀is the maximum transmission at θ=0°, and

[0057] K₉₀is the minimum transmission, or extinction, at θ=90°.

[0058] Applying the trigonometric identity

[0059] sin²θ+cos²θ=1,

[0060] cos²θ=1−sin²θ,

[0061] then

[0062] I/I₀=K₀−K₀sin²θ+K₉₀sin²θ.

[0063] and,

[0064] I/I₀=K₀−(K₀−K₉₀) sin²θ.  Eq.(2)

[0065] Using another identity:

[0066] cos(2θ)=1−2 sin²θ

[0067] or

[0068] sin²θ=1/2−1/2 cos(2θ)  Eq.(3)

[0069] then

[0070] I/I₀=(K₀+K₉₀)/2+(K₀−K₉₀)/2 cos(2θ).  Eq.(4)

[0071] Because we are operating at the 90° extinction point of the crossed polarizers, let

[0072] θ=(δ+π/2) or 2θ=(2δ+π), so

[0073] I/I₀=(K₀+K₉₀)/2+(K₀−K₉₀)/2 cos(2δ+π).  Eq.(5)

[0074] Using the identity

[0075] cos(2δ+π)=cos 2δ cos π−sin 2δ sin π

[0076] cos(2θ)=cos(2δ+π)=−cos(2δ),

[0077] then

[0078] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2 cos(2δ),  Eq.(6)

[0079] where δ is the modulation angle measured from the 90° extinction point of the crossed polarizers.

[0080] Note that applying the identity in Eq.(3), to Eq.(6) yields the mathematically identical expression

[0081] I/I₀=K₉₀+(K₀−K₉₀) sin²δ.  Eq.(7)

[0082] The series expansion of cos(2δ) is

[0083] cos 2δ=(1−2δ²+2/3 δ⁴−4/45 δ⁶+2/315 δ⁸+ . . .).  Eq.(8)

[0084] Now if we let

[0085] δ=β+Ψ_mcos(ωt),  Eq.(9)

[0086] where β represents the polarization rotation angle that we wish to measure (from extinction), and Ψ_mcos(ωt) represents sinusoidal modulation of the polarized light due to an optical modulator. ω represents the modulation frequency in radians. Ψ_mis a constant equal to the maximum angle of a given modulation from extinction. (The standard literature uses a to represent the measured polarization rotation angle and {α}to represent the specific rotation of a substance. β is used herein for the measured polarization rotation angle in order to avoid confusion.)

[0087] The binomial expansion of Eq.(9) is

[0088] δ²=β²+2βΨ_mcos(ωt)+Ψ_m²cos²(ωt)

[0089] δ⁴=β⁴+4βΨ_m³cos(ωt)+6β²Ψ_m²cos²(ωt)+4βΨ_m³cos³(ωt)+Ψ_m⁴cos⁴(ωt)

[0090] δ⁶=β⁶+6β⁵Ψ_mcos(ωt)+15β⁴Ψ_m²cos²(ωt)+20β³Ψ_m³cos³(ωt) +15β²Ψ_m⁴cos⁴(ωt)+6βΨ_m⁵cos⁵(ωt)+Ψ_m⁶cos⁶(ωt)

[0091] δ⁸=β⁸+8β⁷Ψ_mcos(ωt)+28β⁶Ψ_m²cos²(ωt)+

[0092] 56β⁵Ψ_m³cos³(ωt) +70β⁴Ψ_m⁴cos⁴(ωt)+

[0093] 56β³Ψ_m⁵cos⁵(ωt)+28β²Ψ_m⁶cos⁶(ωt)

[0094] +8βΨ_m⁷cos⁷(ωt)+Ψ_m⁸cos⁸(ωt)

[0095] Substituting this binomial expansion into Eq.(8) we have

[0096] cos 2δ={1−2β²+2/3β⁴−4/45β⁶+2/315β⁸}−{4βΨ_m+8/3β³Ψ_m−8/15β⁵Ψ_m+16/315β⁷Ψ_m}cos(ωt)

[0097] −{2Ψ_m²−4β²Ψ_m²+4/3β⁴Ψ_m²−8/45β⁶Ψ_m²}cos²(ωt) +{8/3βΨ_m³−16/9β³Ψ_m³+16/45β⁵Ψ_m³}cos³(ωt)

[0098] +{2/3Ψ_m⁴−4/3β²Ψ_m⁴+4/9β⁴Ψ_m⁴}cos⁴(ωt) −{8/15βΨ_m⁵+16/45β³Ψ_m⁵}cos⁵(ωt)

[0099] −{4/45Ψ_m⁶+8/45β²Ψ_m⁶}cos⁶(ωt) +{16/315βΨ_m⁷}cos⁷(ωt)

[0100] +{2/315Ψ_m⁸}cos⁸(ωt) − . . .  Eq.(10)

[0101] Recognizing the following series:

[0102] sin 2β=(2β−4/3β³+4/15β⁵−8/315β⁷+ . . .)  Eq.(11)

[0103] cos 2β=(1−2β²+2/3β⁴−4/45β⁶+2/315β⁸+ . . .)  Eq.(12)

[0104] We can factor Eq.(10) to produce

[0105] cos 2δ=cos 2β−2Ψ_m{sin 2β}cos(ωt) −2Ψ_m²{cos 2β}cos²(ωt) +4/3Ψ_m³{sin 2β}cos³(ωt)

[0106] +2/3Ψ_m⁴{cos 2β}cos⁴(ωt) −4/15Ψ_m⁵{sin 2β}cos⁵(ωt) −4/45Ψ_m⁶{cos 2β}cos⁶(ωt)

[0107] +8/315Ψ_m⁷{sin 2β}cos⁷(ωt) +2/315Ψ_m⁸{cos 2β}cos⁸(ωt) − . . .  Eq.(13)

[0108] For the higher power cosine terms of Eq.(13) with only one or two terms containing β, it is not evident that β follows a series. We have calculated the binomial expansion of Eq.(8) for 7 terms of the cos 2δ series which provides Eq.(13) out to the cos¹⁴(ωt) term. The β terms in cos⁶(ωt) through cos⁸(ωt) are actually alternating series of cos 2β and sin 2β.

[0109] Applying the following expressions of cosine powers as cosines of multiple angles,

[0110] cos²x=1/2(1+cos 2x)

[0111] cos³x1/4(3 cos x+cos 3x)

[0112] cos⁴x1/8(3+4 cos 2x+cos 4x)

[0113] cos⁵x1/16(10 cos x+5 cos 3x+cos 5x)

[0114] . . .

[0115] For odd powers: ${\cos }^{2n-1}x=1/\left(2^{2n-2}\right)\sum _{k=0}^{n-1}\frac{\left(2n-1\right)!}{k!\left(2n-1-k\right)!}\cos \left(2n-2k-1\right)x\left(n=1,2,\dots \right)$

[0116] For even powers: ${\cos }^{2n}x=1/\left(2^{2n}\right)\left\{\sum _{k=0}^{n-1}\frac{2n!}{k!\left(2n-k\right)!}\cos \left(2n-2k\right)x+\frac{2n!}{n!\left(2n-k\right)!}\right\}\left(n=1,2,\dots \right)$

[0117] to Eq.(13) yields

[0118] cos 2δ=(1−Ψ_m²+1/4Ψ_m⁴−1/36Ψ_m⁶+1/576Ψ_m⁸− . . . )cos 2β−(2Ψ_m−Ψ_m³+1/6Ψ_m⁵−1/72Ψ_m⁷+ . . . ) sin 2β cos(ωt)

[0119] −(Ψ_m²−1/3Ψ_m⁴+1/24Ψ_m⁶−1/360Ψ_m⁸+ . . . )cos 2β cos(2ωt) +(1/3Ψ_m³−1/12Ψ_m⁵+1/120Ψ_m⁷− . . . )sin 2β cos(3ωt)

[0120] +(1/12Ψ_m⁴−1/60Ψ_m⁶+1/720Ψ_m⁸− . . . )cos 2β cos(4ωt) −(1/60Ψ_m⁵−1/360Ψ_m⁷+ . . . )sin 2β cos(5ωt)

[0121] −(1/360Ψ_m⁶−1/2520Ψ_m⁸+ . . . )cos 2β cos(6ωt) +(1/2520Ψ_m⁷− . . . )sin 2β cos(7ωt)

[0122] +(1/20160Ψ_m⁸− . . . )cos 2β cos(8ωt) − . . .  Eq.(14)

[0123] The general form of the Bessel Function series expansion is, $\begin{array}{cc}{J}_n\left(x\right)=x^n\sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k{x}^{2k}}{2^{2k+n}k!\left(n+k\right)!}& \text{Eq.(15)}\end{array}$

[0124] and for J_n(2x), we have $\begin{array}{cc}{J}_n\left(2x\right)=\sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k{x}^{2k+n}}{k!\left(n+k\right)!}& \text{Eq.(16)}\end{array}$

[0125] Expanding Eq.(16) for n=0 to 8 yields,

[0126] J₀(2x)=1−x²+1/4x⁴−1/36x⁶+1/576x⁸− . . .  Eq.(17)

[0127] J₁(2x)=x−1/2x³+1/12x⁵−1/144x⁷+ . . .  Eq.(18)

[0128] J₂(2x)=1/2x²−1/6x⁴+1/48x⁶−1/720x⁸+ . . .  Eq.(19)

[0129] J₃(2x)=1/6x³−1/24x⁵+1/240x⁷− . . .  Eq.(20)

[0130] J₄(2x)=1/24x⁴−1/120x⁶+1/1440x⁸− . . .  Eq.(21)

[0131] J₅(2x)=1/120x⁵−1/720x⁷+ . . .  Eq.(22)

[0132] J₆(2x)=1/720x⁶−1/5040x⁸+ . . .  Eq.(23)

[0133] J₇(2x)=1/5040x⁷− . . .  Eq.(24)

[0134] J₈(2x)=1/40320x⁸− . . .  Eq.(25)

[0135] Inspecting Eq.(14), we see that the first series term is equal to J₀(2x) and the following series terms are equal to 2 J_n(2x).

[0136] repeating Eq.(6)

[0137] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2 cos(2δ)  Eq.(6)

[0138] inserting the Bessel Eq.(17-25) into Eq.(14), and the resulting cos(2δ) expansion into Eq.(6), we have

[0139] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2{J₀(2Ψ_m)cos 2β−2J₁(2Ψ_m)sin 2β cos(ωt) −2J₂(2Ψ_m)cos 2β cos(2ωt)+2J₃(2Ψ_m)sin 2β cos(3ωt)

[0140] +2J₄(2Ψ_m)cos 2β cos(4ωt)−2J₅(2Ψ_m)sin 2β cos(5ωt) −2J₆(2Ψ_m)cos 2β cos(6ωt)+2J₇(2Ψ_m)sin 2β cos(7ωt)

[0141] +2J₈(2Ψ_m)cos 2β cos(8ωt)− . . .}.  Eq.(26)

[0142] Multiplying through by I₀and (K₀−K₉₀)/2, we now have the polarimeter output light intensity I expressed in terms of K₀, K₉₀, I₀, Ψ, and β for each individual output frequency (Eq. 27). For β=0 and Ψ=0, all of the frequency terms are zero and the DC term reduces to K₉₀, the extinction of the polarizers, as expected.

[0143] I=I₀(K₀+K₉₀)/2−I₀(K₀−K₉₀)/2{J₀(2Ψ_m)cos 2β}+I₀(K₀−K₉₀){J₁(2Ψ_m)sin 2β}cos(ωt)

[0144] +I₀(K₀−K₉₀){J₂(2Ψ_m)cos 2β}cos(2ωt) −I₀(K₀−K₉₀){J₃(2Ψ_m)sin 2β}cos(3ωt)

[0145] −I₀(K₀−K₉₀){J₄(2Ψ_m)cos 2β}cos(4ωt) +I₀(K−K₉₀){J₅(2Ψ_m)sin 2β}cos(5ωt)

[0146] +I₀(K₀−K₉₀){J₆(2Ψ_m)cos 2β}cos(6ωt) −I₀(K₀−K₉₀){J₇(2Ψ_m)sin 2β}cos(7ωt)

[0147] −I₀(K₀−K₉₀){J₈(2Ψ_m)cos 2β}cos(8ωt) + . . .  Eq.(27)

[0148] The maximum clinical glucose concentration of 500 mg/dL, measured at a wavelength of 633 nm and a path length of 10 mm, yields a rotation angle of β=0.0225 degrees. (Actual human in vivo measurements will generally have longer wavelengths and shorter path lengths, each of which will reduce the maximum rotation angle at 500 mg/dL.) At the worst case angle of 0.0225 degrees, the errors in replacing the sin and cos functions are:

[0149] sin 2β with 2β, error <1×10⁻¹⁰

[0150] cos 2β with 1, error <0.5×10⁻⁷

[0151] With these approximations Eq.(27) becomes

[0152] I=I₀(K+K₉₀)/2−I₀(K₀−K₉₀)/2{J₀(2Ψ_m)}+I₀(K₀−K₉₀){J₁(2Ψ_m)2β}cos(ωt) +I₀(K₀−K₉₀){J₂(2Ψ_m)}cos(2ωt)

[0153] −I₀(K₀−K₉₀){J₃(2Ψ_m)2β}cos(3ωt) −I₀(K₀−K₉₀){J₄(2Ψ_m)}cos(4ωt)

[0154] +I₀(K₀−K₉₀){J₅(2Ψ_m)2β}cos(5ωt) +I₀(K₀−K₉₀){J₆(2Ψ_m)}cos(6ωt) −I₀(K₀−K₉₀){J₇(2Ψ_m)2β}cos(7ωt)

[0155] −I₀(K₀−K₉₀){J₈(2Ψ_m)}cos(8ωt) + . . .  Eq.(28)

[0156] An early paper by E. J. Gillham, “A High-Precision Photoelectric Polarimeter” M. A., Light Division, National Physical Laboratory, Teddington, Middlesex, Journal of Scientific Instruments, Vol. 34, November 1957, investigated the sensitivity of a polarimeter.

[0157] Multiplying Eq.(6) by I₀and taking the derivatives,

[0158] I=I₀(K₀+K₉₀)/2−(K₀−K₉₀)/2 cos(2δ)  Eq.(29)

[0159] dI/dδ=I₀(K₀−K₉₀)sin(2δ), and  Eq.(30)

[0160] d²I/dδ²=2 I₀(K₀−K₉₀)cos(2δ).  Eq.(31)

[0161] Solving Eq.(31) for d²I/dδ²=0 indicates that the maximum sensitivity of dI/dδ is at ±45° modulation from extinction as Gillham stated. Although the maximum sensitivity occurs at a modulation of ±45°, Gillham stated that precision is maximized at small modulation angles. Factors affecting precision are noise fluctuations, e.g. light intensity fluctuations and mechanical vibrations. Gillham spent several paragraphs discussing these tradeoffs and finally recommended 3° as a modulation angle.

[0162] A new sensitivity algorithm will be developed from the following mathematical analysis. Referring to the odd Bessel expansions from Eq.(17-25),

[0163] J₁(2x)=x−1/2x³+1/12x⁵−1/144x⁷+ . . .  Eq.(18)

[0164] J₃(2x)=1/6x³−1/24x⁵+1/240x⁷− . . .  Eq.(20)

[0165] J₅(2x)=1/120x⁵−1/720x⁷+ . . .  Eq.(22)

[0166] J₇(2x)=1/5040x⁷− . . .  Eq.(24)

[0167] If we multiply each odd Bessel function by its order and then sum, the result is equal to one half of the Bessel abscissa,

[0168] ½(2x)=x=1*J₁(2x)+3*J₃(2x)+5*J₅(2x)+7*J₇(2x)+ . . .  Eq.(32)

[0169] Eq.(32) can be shown if we inspect the Bessel series: $\begin{array}{ccccccc}1*J_1\left(2x\right)& =& x& -1/2x^3& +2/24x^5& -5/720x^7& +\dots \\ 3*J_3\left(2x\right)& =& & +1/2x^3& -3/24x^5& +9/720x^7& -\dots \\ 5*J_5\left(2x\right)& =& & & +1/24x^5& -5/720x^7& +\dots \\ 7*J_7\left(2x\right)& =& & & & +1/720x^7& -\dots \\ \mathrm{Sum}& =& x& +0& +0& +0& +\dots \end{array}$

[0170] For the even Bessel expansions from Eq.(17-25), excluding the J₀series,

[0171] J₂(2x)=1/2x²−1/6x⁴+1/48x⁶−1/720x⁸+ . . .  Eq.(19)

[0172] J₄(2x)=1/24x⁴−1/120x⁶+1/1440x⁸− . . .  Eq.(21)

[0173] J₆(2x)=1/720x⁶−1/5040x⁸+ . . .  Eq.(23)

[0174] J₈(2x)=1/40320x⁸− . . .  Eq.(25)

[0175] If we multiply each even Bessel function by the square of one half of its order and then sum, the result is one eighth of the Bessel abscissa squared,

[0176] 1/8(4x²)=1/2x²=1²*J₂(2x)+2²*J₄(2x)+3²*J₆(2x)+4²*J₈(2x)+ . . .  Eq.(33)

[0177] Aiain this can be shown if we inspect the Bessel series: $\begin{array}{ccccccc}{1}^2*J_2\left(2x\right)& =& 1/2x^2& -1/6x^4& +15/720x^6& -7/5040x^8& +\dots \\ 2^2*J_4\left(2x\right)& =& & +1/6x^4& -24/720x^6& +14/5040x^8& -\dots \\ 3^2*J_6\left(2x\right)& =& & & +9/720x^6& -9/5040x^8& +\dots \\ 4^2*J_8\left(2x\right)& =& & & & +2/5040x^8& -\dots \\ \mathrm{Sum}& =& 1/2x^2& +0& +0& +0& +\dots \end{array}$

[0178] For more information on Bessel formulas, see “A Treatise on the Theory of Bessel Functions”, by G. N. Watson, Cambridge University Press, 2nd edition, 1944 and “Bessel Functions For Engineers”, by N. W. McLachlan, Oxford University Press, 1955.

[0179] By using multiple bandpass filters in the time domain, or a Fast Fourier Transform (FFT) into the frequency domain, or other mathematical signal processing methodologies, we can obtain output signals that represent the magnitudes of each of the frequency coefficients in Eq.(28). For an AC coupled system, the DC term is eliminated. If we use the odd frequency terms of Eq.(28) and apply the Bessel formula Eq.(32), we obtain the following new algorithm:

[0180] Algorithm 1:

[0181] A₁=2βΨ_m=1*Z₁cos(ωt)−3*Z₃cos(3ωt)+5*Z₅cos(5ωt)−7*Z₇cos(7ωt) + . . .  Eq.(34)

[0182] where Z₁, Z₃, Z₅, Z₇, . . . represent the magnitudes of the specified odd frequencies.

[0183] A₁is a computed signal from application of this algorithm to the individual odd frequency output signals. The alternating signs arise from the frequency terms, not from Bessel formula Eq.(32).

[0184] If we use the even frequency terms of Eq.(28) and apply the Bessel formnula Eq.(33), we obtain another new algorithm:

[0185] Algorithm 2:

[0186] A₂=(1/2)ω_m²≈1*Y₂cos(2ωt)−2²*Y₄cos(4ωt)+3²*Y₆cos(6ωt)−4²*Y₈cos(8ωt)+ . . .  Eq.(35)

[0187] where Y₂, Y₄, Y₆, Y₈, . . . represent the magnitudes of the specified even frequencies.

[0188] A₂is a computed signal from application of this algorithm to the individual even frequency output signals. The alternating signs arise from the frequency terms, not from Bessel formula Eq.(33).

[0189] From the trigonometric identity

[0190] sin²δ=½−½ cos 2δ  Eq.(36)

[0191] and using only the first term of the cos 2δ expansion from Eq.(7),

[0192] cos 2δ≈1−2δ²  Eq.(37)

[0193] then

[0194] sin²δ≈δ²  Eq.(38)

[0195] repeating Eq.(9)

[0196] δ=β+Ψ_mcos(ωt)  Eq.(9)

[0197] we have

[0198] sin²δ≈δ²=β2βΨ_mcos(ωt)+Ψ_m²cos²(ωt).  Eq.(39)

[0199] Applying the identity cos 2x=2 cos²x−1 to Eq.(39),

[0200] sin²δ≈δ²=β²+2βΨ_mcos(ωt)+Ψ_m²/2+Ψ_m²/2 cos(2ωt).  Eq.(40)

[0201] If we look at the modulated light output of two optically crossed polarizers at a PIN photodiode current-to-voltage operational amplifier output, we can see the physical significance of the four terms of Eq.(40) on an oscilloscope:

[0202] β²: As the analyzer is moved away from extinction, we can see the DC component of the extinction voltage rise due to this term.

[0203] 2βω_mcos(ωt): A sinusoidal term, at the modulating frequency, that represents the product of our desired measurement, β, and the sinusoidal amplitude, Ψ_m. 2Ψ_mcos(ωt) is the carrier for β. Not easy to see except for large β because it is mixed with a double frequency term. It becomes clear if the output of a bandpass filter centered at frequency wt is observed as β is varied.

[0204] Ψ_m²/2: Because the cos(2ωt) term cannot swing below extinction (by the definition of extinction), this DC component “lifts” the cos(2ωt) term above extinction by exactly ½ of the cos(2ωt) magnitude. In an AC coupled system, both this term and β²term are not used.

[0205] Ψ_m²/2 cos(2ωt): This term corresponds only to the Faraday modulation angle, Ψ_m. It is at double the modulating frequency and is independent of β. The independence is clear if the output of a bandpass filter centered at frequency 2ωt is observed as β is varied.

[0206] Now Ψ_mcos(ωt) is the electronic signal applied to the optical modulator which provides an optical rotation of Ψ_m′ cos(ωt). After passing through the sample and the analyzer, the signal is detected (typically with a PIN photodiode and a current-to-voltage operational amplifier). The AC amplifier output is a composite of the two sinusoids: 2βΨ_m″ cos(ωt) and Ψ_m^″2/2 cos(2ωt). The relationship of Ψ_mto Ψ_m′ to Ψ_m″ is a function of the optical modulator characteristics, the PIN detector, its amplifier, and the feedback resistor (and all their temperature coefficients). But it is important to recognize that all of these influences are common to both frequencies 2βΨ_mcos(ωt) and Ψ_m^″2/2 cos(2ωt). For the remainder of this document, we use Ψ_malone as both an ″ input and output variable in place of Ψ_m′ or Ψ_m″. Ψ_mcan be viewed as the carrier for the measured polarization rotation information, β.

[0207] Repeating Eq.(26)

[0208] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2{J₀(2Ψ_m)cos 2β−2J₁(2Ψ_m)sin 2βcos(ωt)

[0209] −2J₂(2Ψ_m)cos 2βcos(2ωt)+2J₃(2Ψ_m)sin 2βcos(3ωt) +2J₄(2Ψ_m)cos 2βcos(4ωt)−2J₅(2Ψ_m)sin 2βcos(5ωt)

[0210] −2J₆(2Ψ_m)cos 2βcos(6ωt)+2J₇(2Ψ_m)sin 2βcos(7ωt) +2J₈(2Ψ_m)cos 2βcos(8ωt)− . . .},  Eq.(26)

[0211] and recognizing the Bessel relationships,

[0212] cos(2Ψ)=J₀(2Ψ)−2J₂(2Ψ)+2J₄(2Ψ)−2J₆(2Ψ)+2J₈(2Ψ)+ . . .  Eq.(41)

[0213] sin(2Ψ)=2J₁(2Ψ)−2J₃(2Ψ)+2J₅(2Ψ)−2J₇(2Ψ)+ . . . ,  Eq.(42)

[0214] Eq.(41-42) are the sources of the alternating signs of the frequencies in Eq.(27-28). Summing the frequencies according to the Bessel relationships in Eq.(41-42), Eq.(26) can be simplified to a closed form expression,

[0215] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2{cos(2Ψ)cos(2β)−sin(2Ψ)sin(2β)}  Eq.(43)

[0216] for β=0, which eliminates all of the odd frequencies, Eq.(43) reduces to:

[0217] I /I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2 cos(2Ψ)  Eq.(44)

[0218] which is identical to Eq.(6)

[0219] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2 cos(2δ),  Eq.(6)

[0220] further, the cos term of Eq.(44) can be described as

[0221] cos(2δ)=cos(2Ψ)=J₀(2Ψ)+Σ(even frequencies)

[0222] The total light energy, cos(2δ), is then equal to the sum of the DC term, J₀(2Ψ), and the even frequencies. FIG. 4 shows how the light is split between the DC term and the even frequencies. FIG. 5 shows Algorithm 2 and the even frequencies, with the DC term eliminated by AC coupling. The computed Algorithm 2 signal, Ψ_m²/2, is equal to the even frequencies for small modulation angles, Ψ_m, and then increases for high Ψ_mas the square function overtakes the sin²δ function. The even frequency energy broad maxima at 77° indicates perhaps a practical modulation limit for these Algorithms.

[0223] In the previous description of the single series term analysis, the DC component, Ψ_m²/2, “lifts” the composite output waveform above extinction. As the modulation angle Ψ_mincreases, the composite output waveform becomes more and more peaked requiring more and more DC component (shown in FIG. 4) to “lift” the waveform above extinction. FIG. 6 shows the composite output waveform containing only even frequencies at 75 degrees modulation.

[0224] Repeating Eq.(43),

[0225] I/I₀=(K₀+K₉₀)/2−(K₀−K₉₀)/2{cos(2Ψ)cos(2β)−sin(2Ψ)sin(2β)},  Eq.(43)

[0226] the derivative with respect to Ψ is,

[0227] dI/dΨ=I₀(K₀−K₉₀){sin(2Ψ)cos(2β)+cos(2Ψ)sin(2β)}.  Eq.(45)

[0228] Applying the identities

[0229] sin xcos y=½{sin(x+y)+sin(x−y)}  Eq.(46)

[0230] cos xsin y=½{sin(x+y)−sin(x−y)}  Eq.(47)

[0231] we get

[0232] dI/dΨI₀(K₀−K₉₀)sin(2Ψ+2β)  Eq.(48)

[0233] Repeating Eq.(30)

[0234] dI/dδ=I₀(K₀−K₉₀)sin(2δ)  Eq.(30)

[0235] again for β=0, Eq.(48) is the same as Eq.(30) showing that the odd frequencies are a derivative function of Ψ.

[0236] FIG. 3 shows a plot of Algorithm 1, the derivative dI/dδ, and the Bessel coefficient J₁(2Ψ) (which peaks at 52.746°). Summing all of the odd frequencies of Eq.(28) yields a curve identical to the first derivative in FIG. 3 with a maximum at 45°. This means that all standard polarimeters, using lock-in amplifiers at a single frequency, cos(ωt), are actually traversing a Bessel J₁(2Ψ) curve with a peak at 52.75°, not at 45°. This single frequency J₁(2Ψ) cos(ωt) curve is approximately equal to dI/dδ for small modulation angles.

[0237] Standard polarimeters gain very little increase in sensitivity as the modulation angle exceeds about 30° and of course the sensitivity actually begins to decrease at high modulation angles. Now we see from FIG. 3 that by using Algorithm 1, maximum sensitivity is no longer limited to 52.75° but continues to linearly increase beyond 52.75°.

[0238] In taking measurements with a polarimeter, we first adjust the instrument to A₁=0=2βΨ_m, or β₀=0. We then insert a sample into the optical path to get a value for the sample, A_1S=2β_SΨ_S, where the subscript “S” indicates a measurement with a sample. If the initial baseline value is not zero, we cannot get a correct value for A_1Sby subtracting the nonzero baseline. This is because Algorithm A₁is a derivative, a rate of change signal, and the system gain changes as the initial baseline value moves away from zero. Gillham solved this problem by adding a second Faraday rotator as a compensator in a null feedback loop to insure that all measurements were made at null. As stated earlier, all of the more recent researchers have continued to use a null feedback loop.

[0239] If we evaluate the coefficients of the odd and even frequencies at maximum β=0.02255°=0.000394 radians, (for a glucose concentration of 500 mg/dL at 633 nm wavelength and 10 mm path length) and at various modulation angles, we have Table 1, shown below.

[0240] In Table 1, the J₁(2Ψ) and J₂(2Ψ) initial terms are normalized to 100% and the higher order terms are scaled to those values. Remembering that we multiply higher frequency magnitudes by 3, 5, and 7 to correct the J₁(2Ψ) magnitude, (4, 9, and 16 to correct the J₂(2Ψ) magnitude), we could have quantization errors when the correction terms greatly exceed the first term. A practical maximum for modulation may be 60° to 75° (75° is also near the maximum light energy point for the even frequencies shown in FIG. 5).

[0241] Gillham was correct; sensitivity alone without sufficient precision is not useful. From FIG. 7 we see the output waveform after passing through a sample with rotation β (β shown by horizontal open arrows, not to scale). It is apparent that the stability of the output signal is dependent upon the stability of the light intensity, I₀, through the crossed polarizers (intensity variation is indicated by the vertical double filled arrows on the left side of the plot). It is also clear that variation in the null point, (extinction), also directly affects the output signal β (null variation is shown by the horizontal double filled arrows in the center).

[0242] Most standard laser diodes have a photodetector integrated with the laser diode that senses the laser intensity at the backside of the device. The photodetector output is used to control the laser diode power supply to provide constant intensity light output. With this standard control system, the intensity specification I₀is approximately 1%.

[0243] In a Russian paper by Korolev et al., “Stabilization of the emission amplitude of semiconductor laser diodes”, American Institute of Physics, 1990, they use an external electro-optic amplitude modulator with polarizers crossed at 45° as a feedback system. With this system, laser amplitude fluctuations were reduced 25 to 30 db.

[0244] In our system, we already have an optical modulator and 90° crossed polarizers, (δ=0°). In addition we already have a signal proportional to intensity: Ψ_m²/2 from Algorithm 2. Rather than directly stabilizing the light intensity fluctuations, if we stabilize Ψ_m²/2, then our Algorithm 1 output signal β (from 2βΨ_m) will be stable. Note that Ψ_mis the carrier for β. To stabilize Ψ_m, we use the reference channel Ψ_m²/2 signal from a two channel optical system as the input to a control system. The control system then varies the amplitude of the optical modulation to compensate for changes in light amplitude. This is done by varying the amplitude of the sinusoidal signal driving the optical modulator.

[0245] There is a fundamental optical difference between modulation at 0° (our system), and modulation at 45° (Korolev system). A modulator operating at 45° is restricted to small modulation values around 45°. The optical sensitivity of the modulation at 0° is inherently greater and in addition our sensitivity algorithm permits much greater modulation (see FIG. 8). We use the parabolic variable, Ψ_m²/2, instead of the linear variable, Ψ_m, as our control signal which also improves our control system sensitivity.

[0246] Optical modulators exhibit drift due to various causes. For example, we have observed fluctuations in the null point when modulation current is first applied to a Faraday rotator optical modulator. We believe this is due to temperature induced stress variations in the Faraday glass with the stress causing polarization rotation changes. Similar temperature or mechanically induced polarization changes can occur in any of the optical components. Because of the stress induced drift caused by environmental temperature changes the optical test cell requires close monitoring and adjustments need to be recognized during the course of the measurement series. A stabilized environment and a calibration of test cells is also needed. Because Ψ_mhas been stabilized by the intensity feedback loop, we can now use the reference channel 2βΨ_msignal for a second feedback system. This loop drives a DC current through the Faraday coil to provide an opposite rotation to the stress induced rotation, thus maintaining our baseline at null, (extinction).

[0247] A two channel optical system typically will have two detectors with an analyzer for each. We already have an initial polarizer common to both channels. We depend upon extinction in the measurement channel, but our β control system is maintaining extinction only in the reference channel. The analyzers in the two channels will be mechanically adjusted as near as possible to the same null point. However they will never be exactly the same, plus the user may jar the instrument causing a mechanical shift which changes the measure channel zero. Instead of controlling the reference channel signal at zero (extinction), we can introduce an offset in that control system. This offset will be set to force the measurement channel to zero (extinction). This will also give us a means of self checking the instrument: if the offset exceeds some predetermined limit, a message can be displayed that the instrument needs service and calibration.

[0248] It may be possible to employ additional electronic signals for the adjustment of the changes that occur in the temperature stress induced changes to the glass measurement cell or to the device proper. The utilization of additional physical structures, ie: filters, phase shift materials or additional constant verdet glass coils and the necessary signal gathering and processing to provide corrective adjustments by current variation or the utilization of mathematical formula to manipulate the data as a correction methods may be feasible.

[0249] An alternative method of balancing the two optical channels is to use a compensating Faraday coil in the measurement channel. It can be used to zero, or null, the measure channel prior to inserting a sample, or it can be used in the null feedback loop described by Gillham.

[0250] For a given Faraday material, both the Verdet constant and the temperature coefficient of the Verdet constant are known. By measuring the temperature of the compensating Faraday element, we can accurately control the current in the Faraday coil to provide a very precise rotation of the polarized light in the measure channel. This is then used as a calibration factor for the instrument for that specific measurement.

[0251] We have seen references by other researchers (our experience also) to letting the system “warm up” for a period of time before taking data. Both intensity and β fluctuations are fundamentally thermal processes that start at <5 Hz and progress over a long period of time to <0.01 Hz. We need to take readings in less than 1 minute after initial power on thus making stabilization control systems a critical element of a practical instrument.

[0252] Gillham first recommended 3° as a modulation angle. Several more recent researchers used 3° modulation (“Noninvasive Glucose Monitoring of the Aqueous Humor of the Eye: Part1. Measurement of Very Small Optical Rotations”, B. Rabinovitch, W. F. March, R. L. Adams, Diabetes Care, May-June, 1982.); 3° modulation (“Microdegree Polarimetry using a Diode Laser for Glucose Detection”, M. J. Goetz, M. D. Fox, R. B. Northrop, IEEE 0-7803-0902, 1992); and 1° modulation (“Noninvasive Glucose Sensing Utilizing a Digital Closed-loop Polarimetric Approach”, B. D. Cameron, G. L. Coté, IEEE 0018-9294, 1997).

[0253] Gillham was the first to propose and use a Faraday coil as a compensator in a closed loop feedback system. The feedback system removes the rotation of the sample and continually operates at a null point (extinction). The value of the sample rotation then is the DC current required to return to null.

[0254] Rabinovitch similarly used Gillham's feedback compensator method with the statement that “higher sensitivity is achieved if the method is made to be null-point, dependent on a feedback mechanism.” Likewise, Goetz and Cote used feedback systems to null. A null feedback system acts to minimize the rate of change of sensitivity (d²I/dδ²), however it does not maximize sensitivity.

[0255] However, by using the Algorithms developed above, and by optically modulating between 60° and 75°, we can increase sensitivity by 20 to 25 times what other researchers have achieved. As noted above, human tissue with a path length of less than 4 mm (including skin), reduces the glucose rotation angle signal by more than 95%, due to scattering. Thus, we have sufficient sensitivity to make this measurement. It is important to recognize that the level of modulation as indicated above is in excess of that needed for measurement in transparent human tissue. I.e. the anterior chamber of the eye which is primarily comprised of aqueous humor. Aqueous humor is comprised of substances which do not have the properties of, say for instance, skin which our and other researchers have demonstrated scattering and absorption problems. The lack of defined structures, red blood cells and platelets provide for this site as a currently preferred measurement location. Since the eye is transparent, modulation levels of less than 1 degree to approximately 5 degrees should be adequate for “eye” measurements, in which the human eye is the sample 22 as shown in FIGS. 10A, 10B, and 14. As the transparency of the sample decreases, high modulation levels are needed to achieve the necessary signal strength for accurate measurements.

[0256] The preferred embodiment for our signal processing is a digital signal processor (DSP) implementing a Fast Fourier Transform (FFT) which calculates the magnitudes of all of the frequency components of our signal. Instead of an FFT, model-based spectral analysis (MBSA) could be used which has the advantage of shorter sampling times for the same frequency resolution. (An alternative implementation would use bandpass filters for each of the frequencies with precision rectifiers to derive magnitude data. This would require 10 to 20 times the number of I.C.'s and 10 to 20 times the cost of a DSP implementation.)

[0257] Increasing the size (i.e., number of bins) of the FFT, (64,128,256, . . . ) increases the selectivity for each of the frequency bins, thereby reducing noise. The signal to noise ratio increases by 3 db each time the size of the FFT is doubled. A practical DSP for this application would be limited by computation speed to a 256 or 512 FFT.

[0258] Because the same DSP will generate the modulation sinusoidal signal and time the input A/D data sampling, we can insure that we have a signal sampling interval exactly equal to an integer number of samples per cycle of the signal at our photodetectors. Further, because all of the frequencies that we need to analyze are integer harmonics of our modulation frequency, they will all have an integer number of cycles in our sampling period. This is the ideal case for FFT analysis as there is no ‘leakage’ of data from one bin to another.

[0259] Instead of just increasing the size of the FFT, which is limited by computational power, we can average FFT's. This method of coherent averaging has the same requirement as above: the phase of the sampled input signal must be the same at the beginning of each measured sample set. Coherent averaging then preserves the amplitude of signals that are synchronous (coherent), while reducing the variance of the noise. We can realize a signal-to-noise improvement proportional to the square root of the number of FFT's averaged (see “Understanding Digital Signal Processing”, Richard G. Lyons, Addison-Wesley, 1997).

[0260] Recognizing the synchronous nature of the output signals of our polarimeter and using these methods of digital signal processing we can noninvasively detect the polarized light rotation of glucose in human tissue even with 95% signal loss due to scattering and absorption.

[0261] Referring now to FIG. 10A,.a presently preferred embodiment of the optics portion of the noninvasive glucose measurement system of the present invention is designated by numeral 10. Optic system 10 includes a laser 2, which can be a laser producing light in the 750 to 900 nanometer range in which human tissue has minimum absorption. (Additionally, visible lasers in the 600-700 nm range could be used in transparent tissue such as the anterior chamber of the eye.) Laser 2 produces a beam 2A that passes through a first polarizer 4. The polarized beam 2B emanating from polarizer 4 passes through an optical modulator 20, which presently is preferred to be a Faraday rotator. It is important to recognize that the utilization of other methods and effects may be utilized in the conceptual approach such as voltage other than current, phase shift as apposed to angular measurements and other materials such as PLZT other than Faraday glass.

[0262] For example, there are several devices that can serve as the optical modulator 20. A Faraday rotator is one such device. A solid state Kerr cell could also be used as an optical modulator. A Kerr cell is a system in which the electric field is applied transversely to the propagating light beam. The medium of the Kerr cell can range from liquids, such as nitrobenzene, to crystalline materials, such as KDP (potassium di-hydrogen sulfate) or amorphous solids, such as PLZT (lead lanthanum zirconate-titanate). A Pockels cell, in which the electric field is applied along the propagation direction in an appropriate crystalline medium, again such as KDP, could also be used as an optical modulator.

[0263] Optical modulator 20 has two terminals 20A and 20B which receive a modulating current produced by the circuitry 100A shown in FIG. 11A or the circuitry 100B shown in FIG. 12. The modulated laser beam 2C emanating from Faraday rotator 20 impinges on a non-polarizing beam splitter 12 transmitting a beam 2E and reflecting a beam 2D. Beam 2D passes through a second analyzer 7A and emerges therefrom as a beam 2H which then impinges on a first detector 8A. (Polarizer 4 and analyzers 7 and 7A can be polarizers which have a high extinction ratio, such as POLAROID sheet polarizer HN22 or CORNING POLARCOR polarizer.) The transmitted beam 2E passes through an in situ human tissue sample 22 such as an ear lobe, or eye if present, and emerges from the sample as beam 2F which passes through a third analyzer 7 as beam 2J, which impinges on a second detector 8.

[0264] FIGS. 10B and 10C show a preferred, more detailed implementation of the polarimeter system of FIG. 10A. The portion including laser 2, polarizer 4, optical modulator 20, and beam splitter 12 are the same as in FIG. 10A. Similarly, beam 2C in FIG. 10B is reflected by beam splitter 12 to provide beam 2D which passes through analyzer 7A as shown in FIG. 10B to produce beam 2H, which impinges on reference channel detector 8A.

[0265] However, as shown in FIGS. 10B and 10C, glass element 18 incorporating complementary reflections to cancel reflection-induced distortions in the optical rotation and having a known temperature-dependant Verdet constant is used as electro-optic control element to guide the light of the measurement channel through analyzer 7. Beam 2F in FIG. 10B emanates from sample 22 in the measurement channel and is reflected by inclined surface 18A of glass guide 18 to produce beam 2V, which is guided to inclined surface 18B to produce beam 2Z. Beam 2Z passes through analyzer 7B to produce beam 2J, which impinges on detector 8 of the measurement channel. A compensating coil 19 is used to zero the measurement channel electronically, by driving compensating coil 19 so as to produce extinction of the laser beam passed by the crossed polarizers with no sample in its measurement path.

[0266] The temperature of the glass guide 18 is measured by means of a temperature sensor 21 (FIG. 10B), the output 21A of which is amplified and used to adjust the current supplied to compensating coil 19 according to the temperature coefficient of the Verdet constant of the glass guide 18. As shown in subsequently described FIG. 12, the analog output produced on conductor 21A by temperature sensor 21 can be digitized by an A/D converter 21B and supplied as a digital input to digital signal processor 82, which takes into account the temperature of glass guide 18 in computing the amount of drive current to be supplied by optical compensating coil driver 87 to the terminals 19A and 19B of compensating coil 19. Because the temperature coefficient of the Verdet constant is linear, compensation is accomplished by applying the slope of the temperature coefficient to the measured temperature and then correcting the current applied to the coil around the Faraday glass so that the rotation of the polarized light remains constant with temperature changes. By maintaining this constant rotation of the polarized light, we have a method of calibrating the entire instrument every time it is powered on. The simple approach of a look up table is adequate in this instance. However, development of appropriate algorithms may be included in the eventual design upon determination of device components and identification of targeted measurement bands or results.

[0267] The use of glass element 18 in the measurement channel allows use of a smaller and more convenient mechanical design to measure glucose concentration in the human body. (See FIG. 15.) The inclined surfaces 18A and 18B of glass guide 18 are inclined at 45 degrees relative to the beam 2V and, using a known technique, are rotated 90 degrees relative to each other about the measurement beam path so as to provide complementary reflections that cancel the 5 distortions of the optical rotation resulting from the two reflections from inclined surfaces 18A and 18B. The use of glass element 18 and compensation coil 19 in the measurement path also allows calibration of the polarimeter with no sample present in the measurement path. This is accomplished by providing a precise known current (for example, one generated from a reference source) through compensation coil 19. That causes a certain amount of rotation of the light passing through glass guide 18, and a calibration number representative of that amount of rotation can be produced in response to the detector output and stored as a calibration number. The stored calibration number then can be used every time a measurement of a sample is made to determine the amount of rotation caused by the sample.

[0268] A single Faraday glass element with integral complementary reflections is not required for the function described above; two separate mirrors incorporating complementary reflections and a separate Faraday glass compensation element can be an alternate implementation. Either metal or dielectric mirrors can be used.

[0269] In accordance with the present invention, the circuitry required, in effect, to stabilize the physical components shown in FIG. 10A and to implement Algorithms 1 and 2 can be provided as an analog circuit including various bandpass filters required for each of the terms in Eq.(28) or by the “frequency domain” circuit of FIG. 12 including a digital signal processor (DSP) which computes Fast Fourier Transforms (FFTs) to obtain the frequency coefficients represented in above Algorithm 1 and Algorithm 2 to obtain β_m.

[0270] Referring next to FIG. 11A, time domain modulation circuit 100A includes first detector 8 and second detector 8A previously shown in FIG. 11. Each detector can be a PIN diode such as a HAMAMATSU 55821-0 silicon PIN photodiode. Current through first PIN detector diode 8 is connected to the summing input of a current-to-voltage amplifier 25, which preferably is an operational amplifier having a FET input stage, such as an LT1113 low noise dual JFET op amp, available from Linear Technology Corporation (LTC). The output of amplifier 25 is connected to the input of an AC coupled amplifier 26, which can be one half of an LT1013 op amp from Linear Technology Corporation. The output of AC amplifier 26 is connected to the inputs of four bandpass filters 28-1, 28-2, 28-3, and 28-4. Bandpass filter 28-1 has a gain of 10 a center frequency of 1 kilohertz. Bandpass filter 28-2 has a gain of 1 and a center frequency of 2 kilohertz. Bandpass filter 28-3 has a gain of 30, thereby implementing the 3×multiplication for Algorithm 1, and a center frequency of 3 kilohertz. Bandpass filter 28-4 has a gain of 4, thereby implementing the 4×multiplication of Algorithm 2, and a center frequency of 4 kilohertz. Each of bandpass filters 28-1,2,3,4 can be a LTC1164-8 switched capacitor or 8 pole bandpass filter. These bandpass filters are respectively clocked by the four signals 71 shown in FIG. 11B, wherein each bandpass filter is clocked at a frequency 100 times its center frequency.

[0271] The outputs of bandpass filters 28-2 and 28-4 are connected to inputs of conventional full wave precision rectifiers 30-2 and 30-4, respectively. The outputs of precision rectifiers 30-2 and 30-4 are connected to the (+) and (−) inputs of an analog summing amplifier 31, which can be an LT1014 op amp from LTC. The LT1014 op amp together with a feedback capacitor also form low pass filter 32. The output of low pass filter 32 is connected to conductor 33, and produces a signal Ψ_m²/2, which is Algorithm 2. The outputs of bandpass filters 28-1 and 28-3 are connected to inputs of synchronous rectifiers 30-1 and 30-3, respectively.

[0272] Synchronous rectifiers 30-1 and 30-3 can be implemented by means of a CD74HCT4053 triple SPDT analog switch from Harris Semiconductor and two LT1014 op amps from LTC, operated to provide a gain of +1 or −1. The rectification is done according to the 1 kHz and 3 kHz digital sync signals rather than at the zero crossing points of the rectified signals. The output of synchronous rectifier 30-1 is connected to the (+) input of summing amplifier 35, and the output of synchronous rectifier 30-3 is connected to the (−) input of summing amplifier 35, which is identical to summing amplifier 31. The output of summing amplifier 35 is connected to the input of low pass filter 36, which is identical to low pass filter 32. The output of low pass filter 36 is connected to conductor 37, on which the signal 2β_mΨ_m, which is Algorithm 1, is produced. Synchronous rectifier 30-1 is clocked by a 1 kilohertz digital sync signal on conductor 41. Synchronous rectifier 30-3 is clocked by a 3 kilohertz digital sync signal on conductor 42.

[0273] Similarly, the second PIN detector diode 8A produces a current representative of the reference channel beam 2H of FIG. 10A, which is converted to a voltage by current-to-voltage amplifier 25A and is separated from its DC components by AC amplifier 26A, the output of which is connected to the inputs of bandpass filters 38-1, 2, 3, and 4 (which are identical to above described bandpass filters 28-1, 2, 3, 4, respectively). Synchronous rectifiers 30-1 and 30-3 and precision rectifiers 30-2 and 30-4 are “repeated” as synchronous rectifiers 40-1 and 40-3 and precision rectifiers 40-2 and 40-4 connected as shown. Connections of summing amplifiers 41 and 45 are identical to those of previously described summing amplifiers 31 and 35.

[0274] The connections of low pass filters 42 and 46 are identical to the connections of low pass filters 32 and 36, respectively. However, the output of low pass filter 42 produces a signal Ψ_R²/2 on conductor 43, which is connected to the input of a PI (proportional, integral) loop control circuit 50, which can be implemented as an LT1013 op amp with an input resistor and a feedback loop including a capacitor in series with a resistor. If a derivative function is needed, it can be easily implemented with the same op amp and additional resistors/capacitors. Another input of PI loop control circuit 50 is connected by conductor 53 to a first potentiometer 52, which produces a stable, selectable voltage on conductor 53. The voltage on conductor 53 is divided down from a reference voltage V_REFproduced by a stable voltage reference circuit 51. The output of PI loop control circuit 50 on conductor 57 is applied to the (+) input of an analog multiplier circuit 62, the other input of which is connected by conductor 61A to the output of a low pass filter 61. Analog multiplier 62 can be an AD633 available from Analog Devices, Inc. Potentiometer 52 can be adjusted to set the magnitude of the AC modulating signal applied to terminals 20A and 20B of optical modulator 20.

[0275] The input to low pass filter 61 is connected to conductor 60, which receives a 1 kilohertz squarewave signal from the digital logic circuitry 66 of FIG. 11C. Low pass filter 61 is an 8 pole switched capacitor filter having a frequency of between 1.2 and 1.5 times the 1 kHz fundamental frequency applied to conductor 60, and can be an LTC1069-1, commercially available from Linear Technology Corporation. Since the first harmonic of a 1 kHz squarewave is at 3 kHz, the output of low pass filter 61 is a sine wave at 1 kHz, by rejecting all higher harmonics.

[0276] The output of low pass filter 46 produces the signal 2β_RΨ_Ron conductor 47, which is connected to the input of a second PI loop control circuit 56 that is identical to first PI loop control circuit 50. A second input of PI loop control circuit 56 is connected by conductor 55 to a second potentiometer 54, which produces a second reference voltage on conductor 55 based on the reference voltage V_REFproduced by reference voltage circuit 51. The output of second PI loop control circuit 56 is connected by conductor 58 to one input of an analog summing circuit 64, the other input of which is connected by conductor 62A to the output of multiplier 62. The output of summing circuit 64 is connected by conductor 64A to the input of a power amplifier 65, the outputs of which are connected to conductors 20A and 20B of optical modulator 20 of FIG. 2. Potentiometer 54 can be adjusted to provide an offset to compensate for differences between the measurement channel circuitry and the reference channel circuitry. This is accomplished by adjusting reference channel potentiometer 54 to adjust extinction of the measurement channel before sample 22 is introduced.

[0277] The digital circuitry required to produce signals 71 and 72 (FIG. 11B) for the bandpass filters and low pass filters in FIG. 11A is straightforward, but is shown in FIG. 11B for completeness. A crystal clock circuit 67 produces a 2.4 megahertz clock signal applied to the divide by 16 counter 68, which produces four output signals of frequencies 150, 300, 600, and 1200 kilohertz as shown. The 300 kilohertz signal provides one of the bandpass filter clocking signals 71. The 2.4 megahertz clock signal is divided by 3 by divider 69 to produce the 800 kilohertz signal that is input to a divide by sixteen counter 70, which produces 50 kilohertz, 100 kilohertz, 200 kilohertz and 400 kilohertz signals as shown. The 100, 200, 300 and 400 kilohertz clock signals are connected to provide the center frequencies of the above mentioned bandpass filters.

[0278] The 150 kilohertz output of counter 68 is applied to an input of divide-by-five counter 73A, the MSB output of which is connected to the input of divide-by-five counter 73B, the MSB output of which is connected to the input of divide by two counter 73C.

[0279] The outputs of the three counters 73A, B, and C are decoded by a decoder circuit 74, the inputs of which can be adjusted by suitable jumper switches to allow digital phase adjustment of the output signal applied to an input of flip-flop 75 in increments of ±1.8°. The output of flip-flop 75 is connected to conductor 42 and provides the 3 kilohertz digital sync signal required by the circuit of FIG. 11A. The 50 kilohertz output of divide-by-sixteen counter 70 is applied to an identical sequence of counters 76A, B, and C, the outputs of which are decoded by a decoder 77 identical to decoder 74 to produce an output applied to the input of flip-flop 78, wherein the phase of output decoder 77 can be adjusted in increments of ±1.8°. The output of flip-flop 78 is connected to conductor 41 to produce the above mentioned 1 kilohertz digital sync signal.

[0280] Referring next to FIG. 12, a frequency domain implementation 100B of the optical modulator driver circuit includes above-described PIN diodes 8 and 8A and their associated current-to-voltage amplifiers and AC amplifiers are shown. The measurement channel includes an analog-to-digital converter 80 having its input connected to the output of AC amplifier 26 and its output connected to a digital data bus of a DSP integrated circuit 82, which can be a TMS 320VC33 32-bit floating point DSP, from Texas Instruments. Similarly, the reference channel includes an analog-to-digital converter 80A having its input connected to the output of AC amplifier 26A and having its digital output connected to the digital bus 83 of DSP 82. The digital bus 83 of DSP 82 is connected to the input of a digital-to-analog converter 84. The analog output of digital-to-analog converter 84 is connected to the input of power amplifier 65. DSP 82 also is connected to a display 85.

[0281] Digital bus 83 also can be connected to the digital input of the second digital-to-analog converter 86, the analog output of which produces a DC drive to an optical compensating coil driver circuit 87 the outputs of which are connected to terminals 19A and 19B of compensating coil 19 of FIG. 10B, if that feature is used. If compensating coil 19 is used, it can be controlled by digital signal processor 82 to provide a value of the DC drive sufficient to zero the measurement channel.

[0282] According to the present invention the frequency domain implementation of the polarimetry instrument is the most preferred embodiment, because it can be implemented using only about a half dozen integrated circuits, in contrast to the approach of FIGS. 11A and 11B, which requires more than forty integrated circuits to provide an implementation of only the first four frequency terms of Eq.(28).

[0283] For DSP 82 to derive suitable values of Ψ_mand β_mto accurately implement the first four frequency terms of equation 28, FFTs of “size” 128 to 512 are preferred. Analog-to-digital converters 80 and 80A can be 50 kilohertz devices, having a digital word width of 16-24 bits. Analog-to-digital converters 80 and 80A need to sample the analog inputs at twice the highest frequency of the highest harmonic that is needed, according to the Nyquist requirement.

[0284] Referring to the flow chart of FIG. 13, DSP 82 first collects digital data from both analog-to-digital converter 80 of the measurement channel and analog-to-digital converter 80A of the reference channel and synchronously generates sine wave data on digital bus 83 which is sampled by the input of digital-to-analog converter 84. The hardware interrupt of DSP 82 is clocked every 8 μs to synchronously generate the 1024 Hz sine wave data and to sample analog-to-digital converters 80 and 80A every 32 μs.

[0285] Then DSP 82 goes to block 92 of the flow chart of FIG. 13 and inputs time domain data from analog-to-digital converters 80 and 80A, and then computes FFTs every millisecond for both the measurement channel and the reference channel, using commercially available software that can be downloaded from Texas Instruments' web site, File: CZCX0032.ASM, entitled “132 Point DIT Rasix-2, Complex FFT”, to compute the FFTs.

[0286] The FFT algorithm transforms the output data from analog-to-digital converters 80 and 80A into amplitudes of each harmonic, i.e., into the frequency domain. Thus, DSP 82 goes to block 93 of the flow chart of FIG. 13 and computes the values produced by Algorithm 1 and Algorithm 2, using the frequency harmonic magnitudes from the FFT frequency bins to obtain values of ψ_m²/2, Ψ_R²/2, 2β_mΨm, and 2β_RΨ_R. The FFT bins can be considered to be analogous to the precision rectifier outputs of FIG. 11A, which then are summed and filtered to provide the Ψ_m²and 2β_mΨ_mterms.

[0287] Next, the operating program of DSP 82 goes to block 94 of the flow chart and computes Z transform feedback loop coefficients to, in effect, simulate the operation of PI loop control circuits 50 and 56 of FIG. 11A. The Z transform is a well known technique which can be used in implementing a control loop in digital form, and permits DSP 82 to adjust the loop coefficients to provide the needed loop stability. The output of the two loops are input to a digital multiplier simulated by DSP 82, one being provided to a multiplication function analogous to block 62 of FIG. 11A and the other being input to a summing function analogous to summer 64 of FIG. 11A. The sine wave generated in block 1 of FIG. 13 then is multiplied by the output of the first simulated control loop and the results are added to the output of the second simulated control loop to adjust the sine wave magnitude and offset and then output the sine wave data to the digital-to-analog converter 84, providing the stabilized digital data representing the desired optical modulator input. The operations of blocks 92, 93, and 94 are performed every 8 to 10 milliseconds, i.e., roughly a hundred times per second.

[0288] Next, DSP 82 goes to block 95 of FIG. 13 and coherently averages a number of FFTs in order to increase the signal-to-noise ratio as previously explained.

[0289] Then DSP 82 goes to block 96 and applies suitable calibration factors for a particular patient which may be needed to calibrate the raw data and make suitable temperature compensation adjustments as necessary. In block 97 of FIG.13, the difference in rotation (β_S−β₀) then is converted to a glucose level in the patient's tissue and displayed in a suitable digital display.

[0290] FIG. 14 shows an alternate implementation of the optics shown in FIG. 10B, wherein the beam 2B from polarizer 4 is collimated by a collimator 13A and then introduced into an optical fiber 14A. Fiber 14A is wrapped around a torus to form a fiber optic modulator 20. (Fiber optic modulators such as modulator 20 or compensating coil 19-1 are described in “An All Fiber Optic Device For Dynamic Rotation of the Polarization Plane of a Light Wave”, by A. Irace et al., IEEE Journal 0-7803-31-09 dated May 1996, pages 684-686.) Numeral 14B illustrates that optical fiber exiting from fiber optic modulator 20 and passing into a conventional fiber optic splitter 12A from which optical fibers 14C and 14D diverge. Optical fiber 14C terminates in a collimator 13B, from which beam 2E emerges and then passes through sample 22, if it is present in the measurement channel. The beam 2F emerging from the sample 22 enters a collimator 13C, which introduces beam 2F into an optical fiber 14E. Optical fiber 14E can pass directly into a collimator 13E, or can make a number of loops around a torus to form compensating coil 19-1 as shown. In any case, beam 2Z emerges from collimator 13E and passes through analyzer 7. Beam 2J emerges from analyzer 7 and impinges on detector 8.

[0291] Optical fiber 14D terminates on collimator 13D, from which beam 2Y emerges and passes through analyzer 7A, and emerges therefrom as beam 2H which impinges on reference channel detector 8A.

[0292] FIG. 15 shows an implementation of the polarimeter system shown in FIG. 10B contained in a housing 101 having a recess 102 to allow an ear lobe or fold of skin to be introduced as sample 22 into the measurement path. A sample-receiving space 106 into which the ear lobe or fingertip is inserted is bounded by an extension 103 from recess 102. A moveable section 104 is positioned in recess 102 and can be adjusted in the directions of arrows 105 to determine the width of the sample-receiving space 106. With the ear lobe in space 106, section 104 is adjusted to press the ear lobe between the opposed sides of section 104 and extension 103. This provides a measurement of the width (path length) of the sample through which the measurement beam passes. The width of the sample is used in the determination of glucose concentration from the amount of optical rotation of the measurement beam as it passes through the sample. A variety of linear transducers can be used to provide a measurement of path length. It will be necessary to provide the appropriate drawing and description of the method to be employed in making measurement through the anterior chamber of the eye. Due to individual differences in epidermis thickness and melanin concentration, calibration of the instrument to a specific user may be required. Beam 2E passes through the ear lobe (i.e., sample 22) and emerges therefrom as beam 2F, which is reflected by surface 18A of glass guide 18 to surface 18B and emerges therefrom as beam 2Z.

[0293] The structure shown in FIG. 15 could be modified to guide the beam 2F through a portion of the anterior chamber of a human eye which constitutes the sample, and the optical rotation caused by the tissue of the anterior chamber accurately indicates the blood glucose concentration of the patient. In this case an optical method of measuring the path length of the beam across the cornea is used to measure path length.

[0294] As mentioned above, the sample 22 can be the aqueous humor in the anterior chamber of the human eye, as shown in FIGS. 16A-16C. The entrance beam 2E in FIGS. 16A-C can be the entrance beam 2E shown in any of FIGS. 10A, 10B, or 14, and the exit beam 2F shown in FIGS. 16A-C can be the exit beam 2F shown in any of FIGS. 10A, 10B or 14. The human eye 150 includes an aqueous humor 150, a cornea 154, an iris 153, and a lens 158 (shown in FIGS. 16B and 16C). The entrance beam 2E enters the left side of cornea 154. In FIG. 16A, the entrance beam is refracted to produce beam 152, which propagates to a spot 155 of iris 153. Beam 152 is reflected from spot 155 as reflected beam 156. Reflected beam 156 propagates through aqueous humor 151 to the right side of cornea 154, exits the cornea, and is refracted to produce exit beam 2F.

[0295] FIG. 16B shows a different path of the beam transmitted through the aqueous humor 151, wherein entrance beam 2E is refracted as beam 152 which propagates to and is reflected from the lens 158. Beam 152 is reflected as beam 156. Reflected beam 156 propagates through aqueous humor 151 to the right side of cornea 154, exits the cornea, and is refracted as exit beam 2F.

[0296] FIG. 16C shows another path of the beam transmitted through aqueous humor 151, wherein the entrance beam 2E enters the cornea 154 and is refracted to produce beam 162. Refracted earn 162 propagates through aqueous humor 151 parallel to iris 153 to the right edge of cornea 154, and is refracted as it leaves cornea 154 to produce exit beam 2F.

[0297] While the invention has been described with reference to several particular embodiments thereof, those skilled in the art will be able to make the various modifications to the described embodiments of the invention without departing from the true spirit and scope of the invention. It is intended that all elements or steps which are insubstantially different or perform substantially the same function in substantially the same way to achieve the same result as what is claimed are within the scope of the invention. For example, the technique described above could be used on drawn blood samples, other fluid samples, and tissue.