Enhanced Resolution Quadrature Encoder INTERFACE

 

 

 

 

 

___________________________________

 

 

A Thesis

Presented to

the faculty of the School of Engineering and Applied Science

University of Virginia

 

___________________________________

 

 

 

 

 

In Partial Fulfillment

of the requirements for the Degree

Master of Science in Electrical Engineering

 

by

 

Zach Buckner

May, 2004


 

 

 

APPROVAL SHEET

 

The thesis is submitted in partial fulfillment of the

requirements for the degree of

Master of Science in Electrical Engineering

 

 

______________________________

AUTHOR

 

 

 

This thesis has been read and approved by the examining Committee:

 

_______________________________

Thesis (or dissertation) advisor

 

_______________________________

 

_______________________________

 

_______________________________

 

_______________________________

 

_______________________________

 

 

 

Accepted for the School of Engineering and Applied Science:

 

 

_______________________________

Dean, School of Engineering and Applied Science

 

May, 2004


Abstract

Quadrature encoders produce intermediate, sinusoidal signals that must be translated into position or angle for use in control and measurement systems.  This translation is performed by an encoder interface, a system whose behavior is complicated by the distortion and noise that commonly affect the intermediate, sinusoidal signals.  Many existing encoder interface designs exist, providing a wide range of tradeoffs.  A new interface is presented which provides enhanced resolution, compensation for distortion and noise, and miniaturized digital design.

The interface was tested by modifying an existing position sensor, which provided sinusoidal quadrature signals at radians per inch.  A prototype of the sensor, modified to use the new encoder interface, proved successful.  Simulation results were obtained using a detailed model of the modified sensor.  The results showed that the sensor system is accurate to approximately 0.9 mils at 90% confidence, 1.4 mils at 99% confidence, and 2 mils at 99.9% confidence, over a wide variety of operational speeds and noise/distortion levels.  This represents a 1250%, 900%, and 600% increases in resolution, respectively, over the optimal performance of the original encoder interface used with the sensor.


Acknowledgements

I would like to thank my advisor, James Aylor, for providing the opportunity for me to work on this project and for his oversight and guidance during my graduate studies.  Thanks to Visi-Trak Worldwide and Virginia’s Center for Innovative Technology for their generous support of this research.  Thanks to Jack Vann, President and CEO of Visi-Trak Worldwide, for his effort to bring this interface to fruition and his continued commitment to the University of Virginia.  Thanks to Michael Reed, for his support, insight and creativity.  Thanks to Harry Powell, who provided much help with the circuit design and layout.   

Special thanks to my wife, Shirley, for her patience, support, and love.  To my Dad, who has given me much encouragement and extremely valuable advice over the years.  To my Mom, who has an inexhaustible supply of fun and love.

God’s Grace (as described in Romans 7 & 8 of the Bible) is the root of all hope and blessing in my life.   To God be the credit for all that is worthwhile in my work.


Table of Contents

Introduction. 1

Review of Relevant Literature. 3

Design. 5

Input and Output Signals 5

Digital Controller 5

Offset Calibration. 6

Measurement 6

Signal Selection. 7

Normalization. 7

Angle Table 9

Phase Translation. 10

Derived from d0. 11

Derived from d90. 12

Phase Register 13

Output Register 13

Phase Predictor 13

Phase Subtractor 14

Overflow Corrector 14

Output 15

Prototype. 16

Encoder Front-End. 16

Prototype Development 17

Prototype Results 25

Timing Analysis. 26

Maximum Lag Time 27

Maximum Acceleration. 27

Software Simulation. 30

Accuracy Defined. 30

Variables 30

Sensor Resolution. 31

Random Variables 31

Simulation Results 33

Conclusions. 36

References. 37

Appendices. 39

Simulation Results 39

1 Inch per Second. 40

10 Inches per Second. 41

100 Inches per Second. 42

1000 Inches per Second. 43

10,000 Inches per Second. 44

100,000 Inches per Second. 45

MSP Algorithm.. 46

main.c 46

cosine.c 47

cosine.c 48

Simulation Source Code 52

SensorSim.m.. 52

Sample.m.. 54

BatchSim.m.. 55

DumpData.m.. 57


Introduction

Quadrature encoders translate rotary or linear motion into electrical signals.  These signals are sinusoidal, ideally of the form:

                                                                                (1)

                                         ,                                     (2)

 where amp represents the common amplitudes of the signals in volts, pos represents the instantaneous position to be measured, and interval represents the distance that corresponds to a complete sinusoid period.  The basic behavior of the rotary encoder is identical, with the position variable replaced by an angular variable.  The term ‘quadrature’ refers to the fixed  phase difference between the two signals. 

Real encoders exhibit non-ideal quadrature signals with distortion and noise, which can be described according to the following equations:

                                                (3)

                  ,              (4)

where  represents additive Gaussian noise with mean 0 and variance  and off0 and off90 represent DC offsets added to the signals.  The two noise processes added to  and  are assumed to be statistically independent.

To obtain meaningful results, quadrature signals must be translated into their corresponding position variable pos.  Since pos can be much larger than a single interval, simple arccosine() or arcsine() functions are insufficient. 

In an attempt to provide a better set of tradeoffs, a new interfacing technique was developed.  It is capable of increasing resolution using the intermediary sinusoidal signals.  The degree of enhancement is a power of two, and has an upper limit dictated by the analog to digital circuitry contained in the design.  The system is simple and lends itself to an implementation that requires minimal analog circuitry and no manual calibration.  The scope of this thesis is to describe the details of this method and to provide an in-depth analysis of its performance.


Review of Relevant Literature

Many existing techniques exist to translate the intermediate sinusoidal signals into the underlying position or angle.  The most basic and lowest resolution solution is to pass both signals through zero crossing detectors, and to increment or decrement a discrete variable , based on the observed state changes.  The outputs from the zero-crossing detection stage are shown in the second and third panels of Figure 1:

Figure 1: Simple Zero-Crossing Interface

As the sinusoidal input signals complete one cycle, radians, the two-bit “zero crossing” signals transition from “11”, to “01”, to “00”, to “10”.  A finite state machine can be used to monitor these transitions, and increment or decrement the pos counter.

While simple to implement, this scheme provides only 4 state changes (zero crossings) for each interval.  This yields a maximum resolution of (interval / 4), insufficient for many applications.  In addition, this technique requires manual calibration circuitry to accurately set the zero-crossing levels to compensate for off0 and off90 distortion constants.

Many other techniques to enhance resolution have been proposed [1-12], but all suffer from a combination of these deficiencies:

·        Many require large chains of analog electronic circuitry, increasing expense and increasing the interface’s sensitivity to electromagnetic noise and environmental changes. 

·        Many circuits neglect to adequately compensate for changes in amplitude from the incoming quadrature signals. 

·        Many systems require manual calibration of analog components such as potentiometers. 

·        Many designs require the size and complexity of the circuitry to scale with the expected resolution enhancement.  For example, a further doubling of resolution requires a doubling in size and complexity of circuit components.

·        Many designs require operations that are difficult and costly to implement, such as digital division operations.


Design

This section formally documents the design of the improved quadrature encoder interface.  It outlines the role of each component in the process, and describes the flow of control and interconnections between components.

Figure 2: System Parts

Input and Output Signals

The input signals, a0 and a90 are shown in the left side of Figure 1.  These signals are described by Equations 3 and 4.  The output signals lead, trail, and valid emulate traditional quadrature outputs, as generated by the simple zero-crossing algorithm described earlier.

Digital Controller

The digital controller manages the flow-of-control for all of the components in the enclosing box in Figure 2.  For simplicity, wires that connect the digital controller to the individual components have not been drawn. 

For the prototype, the controller consisted of an MSP430 Mixed Signal Processor manufactured by Texas Instruments [14], with all of the components that are inside the box in Figure 2 implemented in source code running on the MSP430.  However, many other configurations are possible.  One promising alternative is to make some or all of the components within the box separate, discrete components, and to control the flow-of-control using a finite state machine on an FPGA.

Offset Calibration

The amplified signals, a0 and a90, have an unpredictable DC offset component due to manufacturing tolerances, placement tolerances, and variation in magnetic flux intensity.   Therefore, offsets must be corrected to a fixed level of 1.5V to assure that a0 and a90 are correctly sampled by the analog-to-digital circuitry.   Precise DC offsets are attained using a simple control loop.  For signal a0, the control loop consists of the analog-to-digital component ADC0, the averaging component Avg0, the digital-to-analog component DAC0, and the subtractor Sub0. Calibration is performed only when power if first applied to the interface, a process managed by the digital controller.  During calibration, a 0 volt signal is initially outputted by Avg0.  Thus, the output from Sub0 follows its input, a0.  The ADC0 then samples a0 at 100 kilosamples per second, and Avg0 performs a time-average of these samples over an 0.25 second interval.  The component’s output, labeled err0 in Figure 2, is generated by subtracting the resultant average value from 1.5 volts:

                                                                                                        (5)

After this stage, the subtractor will remove this error component from subsequent analog values of a0.  Hereafter, the calibrated signal s0 will have a precise DC offset of 1.5 volts.  By the same process, signal s90 will also have a correct DC offset of 1.5 volts.

Measurement

Two analog-to-digital converters, labeled ADC0 and ADC90 in Figure 2, translate analog signals within the range of 0V to 3V to 12-bit signed integer values within the range -2048 to 2047.  After conversion, calibrated signals d0 and d90 vary sinusoidally with respect to pos.  These signals are centered at a DC offset of 0, the middle of the 12-bit sampling window.

Signal Selection

After measurement, digital signals d0 and d90 are passed to signal selection component label Sel in Figure 2.  This component simply compares the absolute magnitudes of d0 and d90 and returns the signal with the lesser magnitude as its primary output, labeled p in the diagram, and returns the signal with greater magnitude as its secondary output, labeled s in the diagram.  The output signal labeled sel in Figure 2 is a digital signal generated inside the component: it is zero when d0 is the primary signal and one when d90 is the primary signal.

Normalization

The normalization component estimates the amplitude of the incoming primary and secondary signals.  According to the ideal form of the input signals described in Equations 1 and 2 and neglecting the difference in amplitudes between the two signals, the signals p and s are simple sine and cosine functions:

                                                                                                       (6)

                                                                                                      (7)

Thus, according to the trigonometric identity

                                                                                                          (8)

we can write

                                              (9)

Rewritten, the equation becomes

                                                                                                       (10)

Because of the 12-bit sampling window in the ADC hardware, signals p and s have a maximum amplitude of .  This assumptions leads to a simplified algorithm using a lookup table to perform the square root operation.  The lookup table works as follows.  Elements are stored in a read-only memory with addresses 1,2,…,2048.  The value at each memory address is equal to the square of the address.  For example the value at address 4 is 16, the value at address 5 is 25. Thus, for a given square x, the task of finding  is reduced to a 11-step binary search.

After the amplitude is estimated, the normalization component then scales the primary signal, p according to the following equation:

                                                                                                                (11)

In this equation, np is the output of the component, labeled in Figure 2.  Again, a lookup table is used to eliminate the costly division operation.  The table works as follows.  Elements are stored in a read-only memory with addresses 1,2,…,2048, each address corresponding to a possible amplitude as calculated above.   The value at each memory address is equal to:

                                                    ,                                              (12)

where  represents the round down, or floor operation.  The resulting memory value is then multiplied by p and the result is shifted right by 10 bits.  This final ‘right shift’ operation is a fast, digital technique to divide by 210.  Thus, the equation for np becomes:

                                                                  (13)

This technique maps p, a sinusoidal signal with unknown amplitude, into a signal with known amplitude, 2048.  Signals p and s are shown in Figure 3 for the case when p corresponds to a0.

Figure 3: Normalized Signals p and s

Angle Table

The angle table is a read-only memory that maps sinusoidal input signal np into a phase angle, a discrete interval on the phase number line.  The angle table follows this equation

                                         ,                                    (14)

 

where resolution is the number of subintervals within the range [0, 360º].  Constant resolution ultimately dictates the resolution of the interface; it represents the number of discernable output events for every interval.  The prototype interface, described later,  uses a resolution of 32, although any other resolution constant would suffice.

The equation listed above requires the arccosine operation, a function whose domain is [-1,1] and range is [0,180º].  The expression (np/2048) maps the signal np, which varies from[-2048,2047], to the necessary [1,1] interval.  The multiplication by (resolution / 360º) constant and the subsequent floor operation maps the result of the arccosine operation onto the integer interval [0, (resolution/2) – 1]. 

Although the above equation describes the nature of the operations, the implementation is greatly simplified through the use of a lookup table, which eliminates the costly arccosine and division operations.  The table works as follows: Elements are stored in a read-only memory with addresses –2048,-2047,…2046, 2047 corresponding to possible values of the signal np.   The value at each memory address is equal to the result of the entire equation defined above.  In other words, the entire equation is reduced to a single memory retrieval operation.

Phase Translation

The phase translation component performs two important corrections on signal ang:

1.      Because the cosine function is not one-to-one, its inverse function arccosine only resolves to the range [0, (resolution/2) – 1], corresponding to the phase interval [0, 180º].  Thus, a second step is necessary to correctly resolve the phase angle to the interval [0, resolution –1] corresponding to the entire phase interval [0, 360º].  The ambiguity of a single sample of p is highlighted in Figure 4:

Figure 4: Phase Translation.  Two possible phases for input sample=900.

2.      If signal ang is derived from signal d90 (a sinusoid that is 90 degrees ahead of d0), the resulting phase angle must be shifted forward by 90 degrees.  Figure 3 summarizes the relationships that exist between d0 and d90. 

Both of these corrections are described below. The phase translation component’s behavior is broken into several cases:

Derived from d0

If signal sel (an output from the signal selection component) is zero, then d0 was selected as the primary signal and ang was derived from d0.  As outlined earlier, this case requires that signal ang be correctly shifted into the interval [0º, 360º].   The corrected phase interval can be determined by inspecting the sign of signal d90, which appears at input s.  If d90 is greater than 0, then the phase is correctly situated in the interval , as shown in Figure 5.

Figure 5:  Phase Translation.  Resolving the two possible phases for input sample 900 using the sign of the secondary signal.

If d90 is less than 0, then the phase should be corrected to be in the interval , as shown in Figure 3.   These two cases lead to the following equation:

                                                                              (15)

In this equation, the symbol act represents the output of the phase translation unit. 

As an example, consider the situation illustrated in Figure 4.  The figure depicts a normalized sample of 900 appearing at signal np.  The angle table will map this to a phase angle of 64º, although it is possible that that the correct angle is 296º (since there are two possible phase angles that correspond to value 900).  To determine the correct angle, the value of the secondary signal s is inspected.  This is shown in Figure 5.  If the value of s is positive, then the correct angle is 64º.  If it is negative, the correct angle is 360º – 64º = 296º.

Derived from d90

If sel equals one, then ang was derived from d90.  The correct phase interval can be determined by inspecting the sign of signal d0, which appears at input s.  If d0 is less than or equal to zero, then only the 90 degree correction is necessary.  If d0 is greater than zero, then two corrections are necessary.  Signal ang is first mapped into the interval .  Additionally, the signal is shifted forward by 90 degree correction.  These two cases lead to the following equation:

                             (16)

where % is the modulo division operator.  After these corrections, signal ph1 correctly represents the phase angle encoded as an unsigned integer between the range [0, resolution -  1].

Phase Register

Following each sensing iteration, the output of the phase register (ph2) is set to match the output of the phase translation component (act).  Thus, the phase register always holds the phase value from the previous iteration. 

Output Register

Following each sensing iteration the Output Register is set to match the final, digital output value, labeled out in Figure 2.  Thus, the delta register always holds the discrete change in position from the previous iteration.  This value is labeled outold in Figure 2.

Phase Predictor

The phase predictor uses the phase from the previous iteration (ph2) and the total movement from the previous iteration (del2) to predict the phase that will result from the current sensing iteration.  The phase predictor consists of a subtractor and a modulo division unit, and yields the following behavior:

                                                                                     (17)

Phase Subtractor

This phase subtractor component determines the amount of angular movement for the current iteration, relative to the predicted phase angle.    This movement is simply the difference between the actual phase value, act, and the phase value for the phase predictor, pred.  The resulting digital output from this component is the signal labeled m  in Figure 2. 

Overflow Corrector

Phase angle transitions from 3590 to 00 or from 0º to 359º yield an apparent movement of –359º or 359º, respectively, though an actual movement of only 1º or –1º, respectively.  The following technique corrects this erroneous overflow/underflow condition:

                                                                 (18)

where m is the apparent movement labeled in Figure 2.

The adder that follows the overflow correction stage is used to compute the total change in position for the current iteration (which is equal to the change for the last iteration, plus the change due to acceleration, labeled mc in Figure 2.

Output

Figure 6: State Machine for Output Signalling

The output component generates three output signals labeled lead, trail, and valid in Figure 2.  These signals are toggled at the end of each sensing iteration.  First, the lead and trail signals are toggled to emulate a traditional quadrature output (which is usually generated by zero crossing detection circuits from sinusoids in quadrature).  In this invention, a simple finite state machine controls the lead and trail signals, as described in Figure 5.  The signal mc determines the number of transitions to follow along the state machine, at 100 nanoseconds per state.  If mc is positive, then “forward” transitions are followed;  if mc is negative, then “backward” transitions are followed. 

After this signaling is complete, the valid signal is raised high for 100 nanoseconds, then returned back to zero.  This signal is used to indicate to decoding logic that the output signaling has finished.  While the output component could be implemented as a separate external component, this invention implements it in software within the MSP430 digital controller.


Prototype

A physical prototype was built and tested to validate the interface’s design principles and lay the groundwork for subsequent improvements.  The prototype consisted of a modified Hall effect position sensor, originally manufactured by Visi-Trak Worldwide, an industrial sensing and controls company based in Cleveland, OH. 

Encoder Front-End

The analog front-end generates the sinusoidal signals characteristic of quadrature encoders.  A schematic of this front-end is shown below:

Figure 7: Quadrature Encoder Front-End for Visi-Trak Sensor

The hydraulic ram has regularly spaced bands of magnetic material around the girth of the ram.  These bands are spaced at intervals of 1/20th of an inch.  Four Hall effect magnetic sensors, labeled Hall1, Hall2, Hall3, and Hall4 in Figure 7, are mounted in close proximity to the ram.  These sensors are spaced at precise distances with respect to one another, such that the first sensor is at 0 degrees, the second at 90 degrees, the third sensor at 180 degrees, and the fourth sensor at 270 degrees.  These angular measurements are relative to the variable interval first introduced in Equation 1, where 360 degrees is equal to one interval, as in the following equation:

                                                                                                          (19)

As the ram moves, the four Hall sensors are fixed in position, such that the alternating bands of magnetic material pass beside the Hall effect sensors as the rod moves along it’s long dimension.  Thus the Hall effect sensors provide four time-varying, sinusoidal Hall voltages, with precise phase relationships.  If variable pos is assigned to be the linear position of the ram (where the origin, pos=0, is assigned to be directly overtop a band), then the following relationships exist between p and the hall voltages:

                                                                                                                        (20)

                                                                                                (21)

                                                                                            (22)

                                                                                            (23)

Hall voltages h1 and h3 are passed through a differential amplifier (labeled Amp0 in Figure 7), which provides both amplification and common-mode noise rejection.  The output voltage for Amp0, labeled a0, follows the equation:

                                                                                            (24)

Thus, the output can be expressed

                                                                                                       (25)

Likewise, signals h2 and h4 are connected to Amp90, yielding:

                                                                                                      (26)

Prototype Development

The first stage of the prototype development was to identify the hardware that would serve as the Digital Controller, as labeled in Figure 2.  The Digital Controller is the most significant element of the design, due to the complex flow of control requirements of the designs.  Specific requirements included:

·        Must be able to store the time sequence behavior (either as compiled machine code or as finite state machine) in non-volitile memory.

·        Must have additional nonvolatile ROM to store the lookup tables for performing the arccosine and square-root operations.

·        Must support integer multiplication and addition

·        Must be small enough to be retrofitted to an existing sensor design with a limiting board dimension of .5 inch (12.7mm).

Several alternatives were considered, but the MSP430 Mixed Signal Processor manufactured by Texas Instruments [14] proved to be the best fit.  In addition to the design requirements, the chip was low cost and provided adequate analog-to-digital conversion onboard.

After initial experimentation with the MSP430 chip using the evaluation kit provided by Texas Instruments, development proceeded in several stages.  The first stage used a breadboard, containing many of the components that would later appear in the full-scale prototype.  However, because the design requires accurate positioning of the Hall effect sensor elements, it proved necessary to reuse the solder pads from the existing Visi-Trak sensor circuit board.  The traces from these solder pads were cut, and wires extended from the board to connect to breadboard circuitry.  The following diagram shows the harness containing the modified Visi-Trak sensor:

Figure 8: Test Harness for First Generation Prototype

Not shown in this picture is the array of laboratory equipment needed to sustain its operation, including the breadboard, 2 DC power supplies (one for the Hall effect elements, one for the MSP430 processor), and the MSP430 evaluation board.  Although this early stage revealed several additional design challenges, it nonetheless demonstrated correct behavior. 

·        Following this successful breadboard prototype, a schematic was created for a full-scale prototype.  This schematic is shown in Figure 9:

Figure 9: Schematic of Full Scale Prototype

 

 

Due to the multifaceted use of many of the components selected for the prototype, the schematic has a one-to-many correspondence with the system parts identified in Figure 2 and Figure 7.  Specifically:

·        The amplifier stage of the schematic includes both Amp0 and Amp90 components from Figure 7, as well as the subtractor elements Sub0 and Sub90 from Figure 2. 

·        The PWM filters shown in the schematic represent the digital to analog converters from Figure 2.  The MSP430 does not provide analog outputs, but instead provides a flexible, on-processor timer circuit.  This timer is used to drive Pulse-Width Modulated (PWM) digital outputs, such that:

                                                ,                                          (27)

where supply voltage represents the voltage that powers the MSP430 processor, 3.4V.  The PWM filters are simple low-pass filters (with a time constant near 1/4 second) designed such that the output voltage is a close approximation to the desired voltage.

The MSP430 performs the functionality for all of the components listed in the “Digital Controller” group in Figure 2.  The MSP430 has an onboard nonvolatile flash memory that contains the source code for each component.  This source code is listed in the Appendix.

After layout, the design files were sent to a manufacturer that specializes in flexible circuit board production.  While flexible circuit board is not required for the encoder interface, it was the base material used for previous generations of the Visi-Trak sensor.  To enable the prototype to use the same packaging canister as its predecessor, the existing board outline and flexible film were used.  The film chosen for the sensor was Kapton, a flexible and durable film manufactured by DuPont [13].  The following picture shows the resulting, populated circuit board:

Figure 10: Populated Full-Scale Prototype

This board is the same physical size as the existing sensor, allowing it fit the existing stainless steel canister used to package it.  The tape visible on the bottom leads of the MSP430 (the rightmost IC) protects the leads from accidental contact with the stainless steel canister after installation.  These leads were not used in the design, and because of extremely tight constraints on the board height, these leads direct abutted the board outline.  The tape is an added measure of protection, but should be unnecessary – all of the leads on the bottom are digital I/O lines that are configured at startup to a high impendence input mode.

The sensor has fold lines generated by cross hatching the copper on the bottom layer of the board.  These fold lines allow the sensor to be compressed to roughly 40% of its original length.  The folded sensor is shown in the following figure, with a quarter to emphasize the sensor’s scale:

Figure 11: Folded Full-Scale Prototype

Initial experimentation with the prototype revealed three design errors:

1.      No pads were included on the circuit board to allow the MSP430 processor to be reprogrammed after it was soldered to the board.  This was a mere design oversight, but greatly complicated the debugging and testing of the sensor.

2.      The calibration signals from the PWM filters in Figure 9 were unbuffered before entering the amplifier inputs sens1ref and sens2ref.  The high source impedance from the filters failed to meet the requirements outlined in the LT1167 Specifications.  The experiments with the breadboard prototype failed to reveal this error because the amplification stages and calibration stages were tested individually, but never together.

3.      A PWM filter input line was inadvertently connected to the wrong pin from the MSP430.

These errors were overcome through a series of tedious assembly steps.  Traces on the board were broken using an Xacto knife and a microscope, and additional traces were added using soldered, fine gage wire.  During the debugging stage, these wires connected to off-board buffering circuitry and an external, reprogrammable socket for the MSP430.  This stage is shown in the following diagram:

Figure 12: Prototype Sensor with External Buffer and Processor Socket

After the initial debugging stage, the prototype was miniaturized by affixing the external buffering circuitry directly on the back surface of the board and eliminating the external processor socket.  This enabled the prototype sensor to fit entirely into the original canister packaging, as shown in the following diagram:

Figure 13: Fully Assembled Prototype Sensor with Target Bar

Prototype Results

The prototype correctly interfaced the Hall effect sensor front end, and reliably and accurately interpolated to resolutions 32 times the underlying interval resolution.  The sensor passed a variety of laboratory tests, and proved successful on a full-sized hydraulic ram at the Visi-Trak corporate headquarters.

Due to the earlier mentioned design error, however, it was impossible to change the program on the on-board DSP without reassembling a new sensor.  These prototypes required extensive manual assembly, making it difficult to test the performance of the sensor over a variety of parameter states.  This necessitated the development of an accurate software simulation, which will be covered in a separate section.


Timing Analysis

Much of the behavior for the prototype sensor is defined by the source code running on the mixed signal processing unit of the interface.  This source code is executed sequentially in an infinite loop, thus the total loop time and the duration of individual sections of this loop have a critical impact on:

·        The time lag between sampling and output signaling

·        Maximum acceleration that the piston can undergo without introducing cumulative error.

·        Realistic software simulation

The timing was analyzed using the assembled source code from the MSP algorithm.  The source code was broken into functional subsections, which are indicated in the source code listing that appears in the Appendix.  

The number of cycles for each assembly instruction varies according to the instruction type and the addressing mode.  Using the data tabulated in the MSP430 User’s Manual [14], the total number of clock cycles for each subsection was tallied.  This tally is summarized in Table 1.  Because the output time depends on the number of output steps per iteration and therefore speed, the table shows the output signaling times at 200 inches per second. 

The table also includes the time expended in each subsection per iteration.  These times are calculated by multiplying the number of cycles by the maximum clock period, 8 MHz.  The final column in the table lists the percentage of total iteration time that each subsection represents.


 


Operation

Cycles

Time (ms)

Time (%)

Amplitude

Setup, Add

12

1.5

2.77%

Multiply

44

5.5

10.16%

Square Root

165

20.6

38.11%

Sample and Process

Read Channels

40

5.0

9.24%

Determine Magnitudes

26

3.3

6.00%

Compare Magnitudes

13

1.6

3.00%

Phase Translation

Setup

10

1.3

2.31%

Multiply

22

2.8

5.08%

Shift

40

5.0

9.24%

Check Bounds

10

1.3

2.31%

Over/Underflow

8

1.0

1.85%

Determine Direction

4

0.5

0.92%

Output (Quadrature)

Setup

11

1.4

2.54%

Loop Maintenance

6

0.8

1.39%

Signaling

11

1.4

2.54%

Activate Valid Line

11

1.4

2.54%

Totals

 

433

54.1

100.00%

Table 1:  Breakdown of Processing Time, Per Section, Per Loop Iteration

Maximum Lag Time

The lag time is the time that elapses between reading the channels and activating the valid line.  Lag time is undesirable – it represents the ‘staleness’ for each reported position.  According to Table 1, the lag time is 42.2 us. As with all of the times reported in Table 1, the output signaling time is reporting for the sensor traveling at 200 inches per second.  Lag times increases with piston speed. 

Maximum Acceleration

The total loop time defines the maximum acceleration that the piston can undergo without introducing cumulative error.  This relationship is due to the discrete time sampling that is performed each loop iteration.  The algorithm uses samples from each iteration to detect phase changes in the underlying sinusoidal signals.  These phase changes are used in turn to estimate the causal change in position.

Specifically, the algorithm uses the difference between the predicted digitized phase, labeled pred in the system diagram, and the actual digitized phase, labeled act in the system diagram.  The difference represents the change in position that is in excess of that predicted change, where the predicted change is based on the velocity from previous iteration. The following equation lies at the heart of this detection:

                                                                              (28)

 

In other words, the change in position for the current iteration is equal the change in position for the previous iteration plus , a step factor due to acceleration.  The following equations relate acceleration and velocity to  and :

                                                                                  (29)


                                                                          (30)

For each iteration,  is known (it is the end result from the previous iteration) and is measured.  This quantity is measured by inspecting the difference between the actual, sampled phase for the current iteration and the predicted phase.  The predicted phase is calculated under the assumption of no velocity:

                                            (31)

 

Thus,

                                          (32)

where the overflowCorrection() function was defined earlier. The phase[pos] is periodic (with a period equal to cosineres) due to the nature of the arccosine operation:

                           phase[pos] = phase[pos + n*resolution], for all integers n                     (33)

and the phase function limited the range [0, resolution – 1].  Thus, it is impossible to meaningfully discern large changes in phase.  Thus, it must be assumed that difference in phases, and thus , will only change in small steps less than or equal to resolution / 2.

                                                                                   (34)

Since relates to acceleration according to Equation 30:

                                                                                                            (35)
Software Simulation

Accuracy Defined

Software simulation was performed to characterize the prototype sensor’s accuracy. The term accuracy refers to the correctness of the position estimate that is provided by each iteration of the interface algorithm.  Specifically, the error for each discrete iteration n is:

                                                                                      (36)

, where posest[n] refers the estimated position and posact[n] refers to the true position.  This error[n] is a random process.  Thus, over a large number of successive iterations of the algorithm, the errors approximate a distribution.  Let this distribution of errors have a probability distribution function called errorpdf(error).  Then let the cumulative distribution function be

                                                                              (37)

The accuracy is defined to be the error value such that

                                                                                      (38)

where confidence is commonly 90%, 99%, 99.9%, etc.  For example, with a confidence level of 99.9% and an accuracy of 1.8 mils, one can expect that 999 out of 1000 estimates provided by the interface will have an error less than 1.8 mils.

Variables

Accuracy is a function of several key variables.  These variables represent the most influential design variables, manufacturing errors, and environmental noise that affect the sensor.  Specifically, these include:

·        Sensor Resolution

·        Amplitude for Signal A

·        Amplitude for Signal B

·        Offset for Signal A

·        Offset for Signal B

·        Phase Error between A and B

·        Noise Level – Bandlimited white noise injection

Each of these variables will be discussed in the sections that follow.

Sensor Resolution

While this constant can be set in the source code for the simulation, for the analysis presented herein, the constant was held at 32.  This resolution provides a theoretical optimal accuracy of .8 mils, a value appropriate for many positioning applications.

Random Variables

Each sensor includes conditioning and calibration components designed to provide a fixed amplitude for each sinusoidal signal, with no DC offset. However, due to installation error and manufacturing defects – the amplitude and offset for a particular sensor instance may not match the desired values.  The effects of these errors were first introduced in Equations 3 and 4 as amp0, amp90, off0, and off90.  In addition, there may exist some deviation from the desired 90˚ phase shift between the quadrature signals, due to fact that the discrete Hall effect elements are manually placed and soldered on the circuit board.  Finally, because the sensor will operate in a harsh industrial environment, additive Gaussian noise may be present on each intermediate, sinusoidal signal.  This noise component was first introduced in Equations 3 and Equation 4 as errnoise. The factor errnoise  represents the variance for a 0-mean, normal distribution used to generate the IID random process, with units of volts.

These magnitudes of these distortions are modeled as random variables, amp0, amp90, off0, off90,errphase, and errnoise.  Because of the limited number of manufactured versions of the improved sensor, there is little data available to characterize these random variables.  However, it is natural to conclude that these random variables are Gaussian random variables, with a mean value equal to the desired value for each corresponding parameter.  It is further assumed that these random variables are statistically independent.  Finally, the standard deviation for each random variable was estimated based on limited experimental observations (estimates based on confidence intervals).   The following table summaries the mean and variance for each random variable:

Error

Mean

StDev

Unit

Amplitude

1.4

.05

Volts

Offset

0

.05

Volts

Phase

0

1

Degrees

Noise

0

.02

Volts

Table 2:  Parameters for Distortion Random Variables

 

 

 


Simulation Results

The simulation was executed for 1000 iterations for each speed, where each characteristic variable (amplitude error, offset error) was chosen according the distributions documented earlier.  Speeds were chosen according to the exponential sequence {1 in/sec, 10 in/sec, 100 in/sec, 1,000 in/sec, 10,000 in/sec, 100,000 in/sec}.  For each speed, a final distribution of errors is obtained, as shown in the following figure for a speed of 1 in/sec:

Figure 14: Probability Distribution Function of Position Error
for Speed of 1 Inch per Second

This is, in turn, corresponded to a cumulative distribution function:

Figure 15: Cumulative Distribution Function of Position Error
for Speed of 1 Inch per Second

This inverse of this function can be used to derive accuracy as a function of confidence, as in the following diagram:

Figure 16: Accuracy Versus Confidence for Speed of 1 Inch per Second

This series of plots is included for all speeds in the Appendix.  The following figure provides a summary of important details from these plots.  The accuracy, as defined in Equation 38, is plotted for three confidence levels (90%, 99%, and 99.9%):

Figure 17: Accuracy versus Confidence for Simulation

This shows that accuracy is not heavily dependent on speed.  The sensor is accurate to 1 mil across all speeds at 90% confidence.  It is accurate to 1.6 mils across all speeds at 99% confidence.  It is accurate to 2.2 mils across all speeds at 99.9% confidence.  This represents a 1250%, 900%, and 600% increases in resolution, respectively, over the optimal, noise-free performance of the original encoder interface used on previous generations of the Visi-Trak sensor.


Conclusions

The encoder interface described in this paper was conceived, modeled, prototyped, and simulated over the course of this research, from October 2002 to December 2003.  The interface represents a major improvement over existing interface designs.  It provides enhanced resolutions over existing encoder interfaces and robust compensation for noise and distortion of the intermediate sinusoidal signals.  Its small size and automatic calibration allow the interface circuitry to be packaged within the housing of the encoder, thus reducing the susceptibility to noise and eliminating the need to package and wire the interface separately.  Because the interface is entirely digital, the output signaling specification is easy to modify, lending the design to further output enhancements like wireless, Ethernet, or proportional voltage outputs.

A provisional patent application for the interface has been filed, and plans to extend the patent protection internationally are underway.  A final design for sensor, which overcomes the errors identified during the development of the full prototype, is complete and currently awaiting manufacture.


References

[1]     W. Papiernik, “Evaluation arrangement for a digital incremental transmitter,” U. S. Patent 4,587,485, May 6, 1986.

[2]     M. Taniguchi, “Signal processing apparatus for pulse encoder with A/D conversion and clocking,” U. S. Patent 4,972,080, November 20, 1990.

[3]     E. L. Griebeler, “High resolution position sensor circuit,” U. S. Patent 5,023,239, April 30, 1991.

[4]     S. Ishii, M. Tsukiji, T. Nishimura, and K. Ishizuka, “Device having signal interpolation circuit and displacement measuring apparatus comprising the device,” U. S. Patent 5,067,089, November 19, 1991.

[5]     S. Ishii and T. Eguchi, “Object displacement detection,” U. S. Patent 5,235,406, August 10, 1993.

[6]     A. J. Santos and M. E. LaCroix, “Resolution multiplying circuit,” U. S. Patent 5,442,313, August 15, 1995.

[7]     N. Kawamata, “Method of and apparatus for detecting an amount of displacement ,” U. S. Patent 5,719,789, February 17, 1998.

[8]     P. L. M. Heydemann, “Determination and correction of quadrature fringe measurement errors in interferometers,” Applied. Optics., vol. 20, no. 19, October 1981.

[9]     N. Hagiwara and H. Murase, “A method of improving the resolution and accuracy of rotary encoders using a code compensation technique,” IEEE Tranactions on. Instrumentation and Measurement., vol. 41, pp. 98–101, February 1992.

[10]   J. R. R Mayer, “High resolution of rotary encoder analog quadrature signals,” IEEE Transactions on Instrumentation and  Measurement., vol. 43, pp. 494–498, June 1994.

[11]   K. K. Tan, H. X. Zhou, and T. H. Lee, “New Interpolation Method for Quadrature Encoder Signals,” IEEE Transactions on Instrumentation and Measurement, Vol. 51, No. 5, October 2002.

[12]   Motion Engineering, “SIM4 Scale Interpolation Module,” Application Notes 206, Rev. C, April 2000, http://www.motioneng.com/pdf/scale_intrp_c.pdf

[13]   DuPont High Performance Materials, "Kapton® polyimide film," March 2004, http://www.dupont.com/kapton/

[14]   Texas Instruments Technical Staff, MSP430x1xx Family User’s Guide, Texas Instruments, 2003.


Appendices

Simulation Results

The following sections contain complete results for the software, across the range of piston speeds: 100,101,102,103,104, and 105 inches per second.  Each section contains the probability distribution function of the position errors, the cumulative distribution function of the position errors, and a summary plot that describes how error increases as a function of confidence level.  In each of the confidence plots, confidence is shown on a logarithmic scale.


1 Inch per Second

 

10 Inches per Second

 

 

100 Inches per Second

 


1000 Inches per Second

 

10,000 Inches per Second

 


100,000 Inches per Second

 

MSP Algorithm

main.c

#include "main.h"

 

void algorithm(void);

 

long offset;

int zeroA;

int zeroB;

 

void calibrate(void) {

  long i;

  long aveA = 0;

  long aveB = 0;

  long iter = 500000;

 

  LCDDISPLAY("Calibrating...", 14);

  LCDSECONDLINE();

  LCDDISPLAY("Move Sensors", 12);

  

 

  for (i=0; i<iter; i++) {

    aveA += ADC12MEM0;

    aveB += ADC12MEM1;

  }

  zeroA = aveA / iter;

  zeroB = aveB / iter;

 

  LCDCLEAR();

}

 

void main(void) {

 

 

  // temporary variable used during calibrations

  int calib;

 

  // only use the watchdog timer if we're using the LCD

  #ifdef LCD

    WDTCTL = WDT_ADLY_1000;               // WDT 1s/4 interval timer

    IE1 |= WDTIE;                         // Enable WDT interrupt

  #else

      WDTCTL = WDTPW + WDTHOLD;             // Stop watchdog timer

  #endif

 

  // setup the clock for fast as she'll go

  DCOCTL = DCO0 + DCO1 + DCO2;

  BCSCTL1 |= RSEL0 + RSEL1 + RSEL2;

 

  LCDINITIALIZE();

 

 

  // PWM ---------------------------- 

   

    TACTL = TASSEL1 + TACLR;              // SMCLK, Clear Tar

    CCR0 = 4095;                                // PWM Period

    CCTL1 = OUTMOD_7;                     // CCR1 reset/set

    CCR1 = 2048;                          // CCR1 PWM duty cycle

    CCTL2 = OUTMOD_7;                     // CCR2 reset/set

    CCR2 = 2048;                          // CCR2 PWM duty cycle

    P1DIR |= 0x0C;                        // P1.2 and P1.3 output

    P1SEL |= 0x0C;                        // P1.2 and P1.3 TA1/2 otions

    TACTL |= MC0;                         // Start Timer_A in up mode 

 

  // ADC ----------------------------

  //

 

    P6SEL |= 0x88;                        // Enable A/D channel A0

    ADC12CTL0 = ADC12ON+SHT0_2+MSC;       // Turn on ADC12, set sampling time

    ADC12CTL1  = SHP + ADC12SSEL_2;      // Use sampling timer, use MCLK

// for conversion timer.

    ADC12CTL1 |= CONSEQ_3;                // Use consecutive mode 

    ADC12MCTL0 = INCH_7;

    ADC12MCTL1 = INCH_3 + EOS;

    ADC12CTL0 |= ENC;                     // Enable conversions 

    ADC12CTL0 |= ADC12SC;                 // Start conversion

 

  // Port ----------------------------

  //

 

    P4DIR |= 0xfF;

 

  // Go ----------------------------

  //

 

    // PWM circuits has just been activated.  Do two calibrations,

    // the first simply acts as a delay so allow the filtered PWM

    // voltage to settle.

    calibrate();

    calibrate();

 

    // NOTE channel1 uses CCR2, vice  versa

    calib = CCR1 - (zeroB - 2048);

    if (calib < 0) calib = 0;

    if (calib > 4095) calib  = 4095;

    CCR1 = calib;   

    calib = CCR2 - (zeroA - 2048);

    if (calib < 0) calib = 0;

    if (calib > 4095) calib  = 4095;

    CCR2 = calib;

 

    // Do another round of calibration.

    calibrate();   

    calibrate();

   

    // NOTE channel1 uses CCR2, vice  versa

    calib = CCR1 - (zeroB - 2048);

    if (calib < 0) calib = 0;

    if (calib > 4095) calib  = 4095;

    CCR1 = calib;   

    calib = CCR2 - (zeroA - 2048);

    if (calib < 0) calib = 0;

    if (calib > 4095) calib  = 4095;

    CCR2 = calib;

   

    // Delay, letting the filtered PWM voltage settle

    calibrate();   

    calibrate();

    

             

    // LCD requires the watchdog timer -- Enable interrupts

    #ifdef LCD

      _EINT();                             

    #endif

 

  algorithm();

 

}


cosine.c

#include "main.h"

#ifdef COSINE

 

#define READ \

    a = ((int) ADC12MEM0) - zeroA;\

    b = ((int) ADC12MEM1) - zeroB; \

    if (a < 0) { aMag = a * -1; } else { aMag = a; } \

    if (b < 0) { bMag = b * -1; } else { bMag = b; }

   

extern const unsigned int div2[];

extern const long squares[];

extern const unsigned char cosine[];

extern long offset;

extern int zeroA;

extern int zeroB;

extern long fastMultiply(int a, int b);

 

 

void algorithm(void) {

  signed char last = 0;

  signed char delta = 0;

  signed char signal = 0;

  signed char current = 0;

  int arrayPos;

  long amp = 800;

  long ampSQ;

  int i;

  int l;

  int r;

  unsigned char quadrature = 0;

  int a, b;

  int aMag, bMag;

 

  READ;

 

  while (1) {

     

   

    ampSQ = fastMultiply(a,a);

    ampSQ += fastMultiply(b,b);

   

    // This next block of code performs a cheap square root

    if (ampSQ < squares[1024])  {

      l = 0; r = 1023;

    } else {

      l = 1024; r = 2047;

    }

   

    if (ampSQ < squares[l+512])  {

       r -= 512;

    } else {

      l += 512;

    }

 

    if (ampSQ < squares[l+256])  {

       r -= 256;

    } else {

      l += 256;

    }

   

    if (ampSQ < squares[l+128])  {

       r -= 128;

    } else {

      l += 128;

    }

   

    if (ampSQ < squares[l+64])  {

       r -= 64;

    } else {

      l += 64;

    }

       

    if (ampSQ < squares[l+32])  {

       r -= 32;

    } else {

      l += 32;

    }

   

    if (ampSQ < squares[l+16])  {

       r -= 16;

    } else {

      l += 16;

    }

   

    if (ampSQ < squares[l+8])  {

       r -= 8;

    } else {

      l += 8;

    }

   

    if (ampSQ < squares[l+4])  {

       r -= 4;

    } else {

      l += 4;

    }

   

    if (ampSQ < squares[l+2])  {

       r -= 2;

    } else {

      l += 2;

    }  

   

    if (ampSQ == squares[l])  {

      amp = l;

    } else {

      amp = l+1;

    }

   

    READ;

       

    // We use either signal a or signal b as the "lookup" source,

    // depending on which signal is smaller (smaller is better)

    if (aMag < bMag) {

 

      arrayPos = (fastMultiply(a,div2[amp]) >> 10) + 2048;     

      if (arrayPos < 0) arrayPos = 0;

      if (arrayPos > 4095) arrayPos = 4095;

 

      current = cosine[arrayPos];

      if (b < 0) current = COSINERES - 1 - current;

   

    } else {

      arrayPos = (fastMultiply(b,div2[amp]) >> 10) + 2048;   

      if (arrayPos < 0) arrayPos = 0;

      if (arrayPos > 4095) arrayPos = 4095;

 

      current = cosine[arrayPos];

      if (a > 0) current = COSINERES - 1 - current;

      current = (current + (COSINERES / 4));

     

      // the previous operation could bump this above the resolution.

      // This is a cheap way to perform a modulo divide by COSINERES

      current &= COSINERES - 1;

     

    }

 

    delta = current - last;

 

    // check for overflow or underflow

    if (delta > (COSINERES / 2)) {

      delta -= COSINERES;

    } else {

      if (delta < -(COSINERES / 2)) {

        delta += COSINERES;

      }

    }

   

   

    // output

    if (delta < 0) {

      for (i = 0; i > delta; i--) {

          quadrature = ((quadrature & 0x01) << 1) | ((~quadrature & 0x02) >> 1); 

          P4OUT = quadrature;

      }

    }

    if (delta > 0) {

      for (i = 0; i < delta; i++) {

          quadrature = ((~quadrature & 0x01) << 1) | ((quadrature & 0x02) >> 1);

          P4OUT = quadrature;

      }

    }   

    P4OUT |= 0x04;

    _NOP();

    P4OUT &= ~0x04;

   

   

                

    last = current;

  }

}

 

#endif
Cosine  Constants

Unlike the previous code written in the C language, the following code is written in Java.  Once executed, it writes a large collection of tables to the screen  or standard output.  This resulting text is then saved in a C file called cosineconstants.c.  These are the tables necessary to compile cosine.c, listed above.

public static void main(String[] args) {

String s = ("unsigned char cosine[] = {");

for (double i=0; i < 4096; i++) {

double y = (i - 2048d) / 2048d; // returns 0 to pi

double x = Math.floor( (Math.acos(y) / Math.PI) * 16d );

s+=  ((int) x) + ", ";

 

// if this is a newline

if (((i + 1) % 64 == 0) & (i != 0))

s += "\n";

}

System.out.println(s);

 

 

s = ("float divisors[] = {\n");

for (float i=0; i < 2048; i++) {

String c = "" + (int) ((2048f / i) * 1024);

for (int j = c.length(); j < 10; j++) { c += " "; }

s+= c + ", ";

 

// if this is a newline

if (((i + 1) % 64 == 0) & (i != 0))

s += "\n";

}

System.out.println(s);

 

s = ("const long squares[] = {\n");

for (long i=0; i < 2048; i++) {

String c = "" + i*i;

for (int j = c.length(); j < 10; j++) { c += " "; }

s+= c + ", ";

 

// if this is a newline

if (((i + 1) % 64 == 0) & (i != 0))

s += "\n";

}

System.out.println(s);

}


Simulation Source Code

SensorSim.m

function errors = sensorSim(cosineres,amp1,amp2,off1,off2,phase,noise,stepSize, simEnd, errorHist, speed, sineAmp, timeConstant);

 

% units

us = 10^-6;

 

% calculated and initialized

numSteps             = floor(simEnd / stepSize);

pos                  = 0;          % the actual position

cPos                 = 0;          % the simulated / calculated position

last                 = 0;          % phase angle for last iteration

delta                = 0;                                                                

%output: the actual output histogram

errors               = zeros(1, length(errorHist));%

sineConstant = sineAmp*simEnd/(pi*2);

expConstant = speed/timeConstant;

 

% perform an initial sample

[sig1, sig2] = sample(pos, amp1, amp2, off1, off2, phase, noise);

 

for step = 1:numSteps

   time = stepSize*step;

  

   amp = floor(sqrt(sig1^2 + sig2^2));

  

   %exponential + sinusoid

   pos = double(speed*time  -  speed*exp(time*timeConstant)/timeConstant   -    sineConstant*cos(4*pi*time/simEnd) + sineConstant + expConstant);

  

   %sinusoidal velocity

   %pos = double(speed*time - sineConstant*cos(2*pi*time/simEnd) + sineConstant);

     

   %ramp

   %pos = (speed / (2 * (simEnd)))*time^2;

    

   %constant

   %pos = speed*time;

  

   [sig1, sig2] = sample(pos, amp1, amp2, off1, off2, phase, noise);

  

   if (abs(sig1) < abs(sig2))

      p = sig1/amp;

      if (p  > 1)

         p = 1;

      elseif (p < -1)

         p = -1;

      end

     

      current = round(cosineres*acos(p) / (2*pi));

     

      if (sig2 < 0)

         current = cosineres - current;

      end

      current;

      current  = mod(current, cosineres);

     

   else

      p = sig2/amp;

              if (p  > 1)

         p = 1;

      elseif (p < -1)

         p = -1;

      end

     

      current = round(cosineres*acos(p) / (2*pi));

       

      if (sig1 > 0)

              current = cosineres-current;

      end

       

      current = (current + (cosineres / 4));

      current = mod(current, cosineres);

       end   

   

    % the predicted change in actual (> resolution) units is equal to the

    % delta from the last iteration.

    predicted = delta;

   

    % using the cosineRes units, figure out how far off our guess is.

        actualOverPredicted  = current - mod(last + predicted, cosineres);

     

    % check for overflow or underflow

    if (actualOverPredicted  > (cosineres/ 2))

      actualOverPredicted  = actualOverPredicted  - cosineres;

     else

      if (actualOverPredicted  < -(cosineres / 2))

        actualOverPredicted  = actualOverPredicted  + cosineres;

     end

    end  

     

    % correct the prediction accordingly. 

    delta = predicted + actualOverPredicted;

     

   cPos = cPos + delta;

   last = current;

  

   error = (cPos / (20*cosineres)) - pos;  

       errors = errors + hist(abs(error), errorHist);

     

end


Sample.m

function [sig1, sig2] = sample(pos, amp1, amp2, off1, off2, phase, noise)

 

sig1 = amp1*cos(2*pi*20*pos) + off1 + randn(1)*noise;

sig2 = amp2*cos(2*pi*20*pos - (pi/2) - (phase*2*pi/360)) + off2 + randn(1)*noise;

 

% emulate 12 bit sampling

sig1 = floor(sig1*4096/3);

sig2 = floor(sig2*4096/3);

 

if (sig1  > 2047)

   sig1 = 2048;

elseif (sig1 < -2048)

   sig1 = -2048;

end

 

if (sig2  > 2047)

   sig2 = 2047;

elseif (sig2 < -2048)

   sig2 = -2048;

end


BatchSim.m

% units

us = 10^-6;

 

% constants

iterations    = 100;

ampM                 = 1.4;

ampSd         = .05;

offSd         = .05;

phaseSd       = 1;

noiseSd       = .02;

timeConsant   = -40;

 

cosineres     = 32;

stepSize             = 50*us;

maxAccel             = 250 / stepSize;

 

simEnd        = .2;

errorMax      = ((1 / 20) / 32) * 5;

errorStep     = errorMax / 1000;

errorHist     = 0:errorStep:errorMax;

 

confData      = zeros(6,4);

 

for s = 0:5

  

    speed = 10^s;  

   sineAmp = speed / 10;

  

   % this should set the time constant so as to ramp up to the target

   % speed with maximum acceleration

   timeConstant = -1 * maxAccel / speed

  

  

    % output: the histogram bins for logging errors

    numErrorSteps    = floor(errorMax / errorStep);

    errorsSum        = zeros(1,numErrorSteps + 1);

 

    for i=1:iterations

  

       if(mod(i-1, 10) == 0)

              i

       end

  

       amp1 = randn(1)*ampSd + ampM;

       amp2 = randn(1)*ampSd + ampM;

       off1 = randn(1)*offSd;

       off2 = randn(1)*offSd;

       phase = randn(1)*phaseSd;

       noise = randn(1)*noiseSd;

  

       errors      = sensorSim(cosineres,amp1,amp2,off1,off2,phase,noise,stepSize,simEnd,errorHist, speed, sineAmp, timeConstant);      

       errorsSum = errorsSum + errors;

     

    end

 

    % dump the data and plots

    errorPDF  = errorsSum / (iterations*simEnd/stepSize);

    paramStr = strcat('res',num2str(cosineres),'_ampM',num2str(ampM),'_ampSd',num2str(ampSd),'_offSd',num2str(offSd),'_phaseSd',num2str(phaseSd),'_noiseSd',num2str(noiseSd),'_speed',num2str(speed));

    paramStrV =

strcat('Res=',num2str(cosineres),', AmpM=',num2str(ampM),', AmpSD=',num2str(ampSd),', OffSD=',num2str(offSd),', PhaseSD=',num2str(phaseSd),', NoiseSD=', num2str(noiseSd),', Speed=',num2str(speed));                 

    conf = dumpData(errorPDF, errorHist, paramStrV, paramStr);

    save(strcat('data_', paramStr, '.mat'), 'errorPDF', 'conf');

  

   confData(s+1, 1) = speed;

   confData(s+1, 2:4) = conf;

  

end


DumpData.m

% errorPDF  - normalized PDF of error distribution

% errorHist - the x-axis labels for the distribution

% paramStrV - Verbose String that describes the plots

% paramStr  - Condensed String that describes the plots

function conf = dumpData(errorPDF, errorHist, paramStrV, paramStr);

 

pdf = errorPDF;

cdf = cumsum(pdf);

 

f9 = find(cdf >= .9);

d9 = errorHist(f9(1));

f99 = find(cdf >= .99);

d99 = errorHist(f99(1));

f999 = find(cdf >= .999);

d999 = errorHist(f999(1));

conf = [d9 d99 d999];

 

figure(1);

bar(errorHist, pdf);

       set(gca, 'xlim', [0 errorHist(length(errorHist))]);

       title(sprintf(strcat('Probability Distribution Function of Position Error\n', paramStrV)));

   xlabel('|Position Error|, mils (.001 in)');

   ylabel('Probability Density');

   generateFigure(strcat('pdf_', paramStr));

  

figure(2);

plot(errorHist, cdf);

       set(gca, 'xlim', [0 errorHist(length(errorHist))]);;

       title(sprintf(strcat('Cumulative Distribution Function of Position Error\n', paramStrV)));

   xlabel('|Position Error|, mils (.001 in)');

   ylabel('Probability');

   generateFigure(strcat('cdf_', paramStr));

  

%figure(3);

%clf

%time = 0:50*10^-6:2*(1/(20*speed));

%cossample = amp1*cos(2*pi*20*speed*time) + off1 + randn(size(time))*noise;

%sinsample = amp2*sin(2*pi*20*speed*time) + off2 + randn(size(time))*noise;

%plot(time,cossample, time, sinsample);

%      Title(sprintf(strcat('Example Sample Path, 2 Cycles\n', paramStrV)));

%   XLabel('Time, Seconds');

%   YLabel('Voltage');

%   print ('-depsc2', '-tiff', strcat('samplePath_', paramStr, '.eps'));

  

figure(3);

plot(1./(1 - cdf), errorHist);

set(gca, 'XLim', [0, 10^4]);

       title(sprintf(strcat('Sensor Accuracy versus Confidence\n', paramStrV)));

   xlabel('Confidence');

   ylabel('|Position Error|, mils (.001 in)');

   set(gca, 'XScale', 'log');

   set(gca, 'XTickLabel', ['90%   '; '99%   '; '99.9% '; '99.99%']);