CN1790916B - Calibration system for a linearity corrector using filter products - Google Patents

Abstract

A calibration system for calibrating a linear corrector using a sum of filter products, and a method of calibrating the linear corrector, are disclosed. The calibration system includes first and second signal generators for introducing a test signal into a signal processing system, such as an ADC. An acquisition memory and processor are provided for acquiring and analyzing the output of the signal processing system, and then programming the filter coefficients into the linearity corrector. The calibration method analyzes intermodulation and harmonic components collected from the signal processing system and then looks for amplitude and phase responses for the filter. The magnitude and phase responses are then used to determine a set of filter coefficients.

Description

Calibration system using linearity corrector of filter product

cross reference

This application claims priority to US Provisional Application No. 60/625,372, filed November 4, 2004. the

Copyright Notice

Pursuant to 37 C.E.R.1.71(e), applicant represents that a portion of this disclosure contains material that is and claims copyright protection, such as, but not limited to, source code listings, screen shots, user interfaces, or user instructions, or any submitted Other aspects for which copyright protection is or may be useful in any jurisdiction. The copyright owner has no objection to the facsimile reproduction of any of the patent documents or patent disclosures, such as appearing in the Patent and Trademark Office patent files or records. All other rights are reserved and all other reproduction, distribution and creation of derivative works based on the content, public display and public performance of this application or any part thereof are prohibited by applicable copyright laws. the

technical field

The present invention relates to linearity error correction, and more particularly to a linearity corrector that uses filter products to reduce or eliminate distortion produced by signal processing systems such as analog-to-digital converters (ADCs). the

Background technique

Reducing the distortion produced by the ADC increases spurious-free dynamic range (SFDR), which is useful for systems that use ADCs to acquire data, such as spectrum analyzers and other electronic measurement instruments. the

Modern high-speed ADC designs use a deep clock pipeline to help accurately convert an analog input into a sampled digital representation in a series of refined steps. ADC designers put a lot of effort into removing significant sources of nonlinearity in the analog processing circuitry. However, it is generally difficult to remove all error sources. Designers try to eliminate the most obvious problems in the circuit until computerized modeling such as SPICE modeling shows that the converter meets specification. Linearity can be improved by using techniques such as reducing the dynamic range of a non-linear device or using feedback around it. However, some circuit topologies have inherent distortion mechanisms that cannot be completely removed. the

Pipelining also provides the opportunity for internal digital and analog circuit activity to modulate internal analog signal processing. In many such cases, the self-adjustment of the input signal as a linear function of itself or its derivatives produces residual non-linear distortion. This resulted in some low-level distortion that was difficult to remove. This modulation can occur through internal power distribution. In this case, the number of circuit channels capable of producing voltage modulation on that supply rail can be quite high. Simulating this effect complicates device modeling and slows down computer simulations. For the first order, these effects on supply modulation will increase almost linearly, so they can be modeled as linear finite impulse response (FIR) filters. the

Modulation occurs at one or more points in analog signal processing, which correspond to multiplications. In pipeline ADCs, modulation typically occurs in high-gain analog amplifiers between conversion stages. In this case, the harmonic and intermodulation distortion is typically characterized by the presence of 2nd and 3rd order distortion terms, with very little higher order distortion occurring. the

Previously proposed solutions were based on Volterra filters. The impulse response of the ADC can be many clock periods, for example 64 clock periods can be used. A correction system using a Volterra filter will require a similar response length. In a 3rd order distorted Volterra system, this results in a filter of approximately (N ³ )/6 taps, which for a corrected system with a response length of 64 would result in an order of approximately 50000 taps. Filter systems with such a huge number of taps were too complex and expensive to be realized in practical systems at the time.

Another solution has been proposed in other applications related to correcting distortion in loudspeakers, which uses a filter structure that approximates some aspects of the Volterra filter. Figure 1 shows a version of this scheme with 1st order correction and 3rd order correction. First order compensation is provided by filter 12(h ₁ ). Third order compensation is provided by multiplying the output of filter 14 with the output of filter 16 using multiplier 18, filtering the output of multiplier 18 using filter 20, and multiplying the output of filter The output of 20 is multiplied by the output of filter 22 and finally the output of multiplier 24 is filtered using filter 26 . By adding the output from the first order compensation and the third order compensation using adder 28, linear cubic compensation can be provided. The 3rd order compensation for the system shown in Figure 1 implements the following equation,

It is described as a general 3rd order nonlinear filter structure. This implementation uses filter 20 after multiplier 18 and filter 26 after multiplier 24 . Once the linear cubic compensation is obtained, it is subtracted from the output of the unknown system being compensated. This requires that the corrector have full access to the raw signal input to the unknown system, which is not available when the raw signal is not digital. While it can be a useful sub-solution of the Volterra filter, it also has the disadvantage of not being suitable for systems with a good linear frequency response. Filters immediately following this multiplier cannot distinguish between original components and aliased components caused by the nonlinear effects of previous multiplications. While additional filtering after this multiplier can provide some correction to the frequency-dependent magnitude and phase response in the signal path, aliasing does not allow this filter to reproduce the original The phase and magnitude response differences between the aliased and aliased components are corrected for. the

A legacy issue in linear compensation systems involves calibration. These systems would require solving a system of filter coefficients that is nonlinear with respect to the output. For any calibration scheme that can be applied to the system, solving for more coefficients will require more computation. the

Details and improvements over existing solutions are discussed in more detail below. the

Contents of the invention

If this distortion mechanism can be modeled by recovering the coefficients of an equivalent distortion filter in an ADC, the ADC output can be passed through a digital processing network that distorts the signal in essentially the same way, and then subtracts the distortion to obtain reduce or eliminate this ADC distortion. Although it is impossible to completely eliminate all ADC distortion, this method improves the spurious-free dynamic range (SFDR) of the ADC. For example, an ADC with an SFDR of 80dB can be improved by a factor of 15dB, depending on the characteristics of the ADC. This improvement also removes some filters from the previously proposed layout, thereby simplifying the design for use in systems with a relatively smooth linear frequency response. This simplification can be alternated with longer filters in the correction system, resulting in improved performance with the same throughput. This improvement is of interest for precision measurement applications, such as those associated with spectrum analyzers, oscilloscopes, or other measurement instruments that use ADCs. the

To achieve these advantages, a suitable calibration system and method are provided. The calibration system includes: a first analog signal generator that generates a first analog signal, a second analog signal generator that generates a second analog signal, and connected to the first analog signal generator and the second analog signal generator, An analog adder for adding the first and second analog signals and providing an analog output. This analog output is input to an ADC to obtain a digital output. A linearity corrector, which performs a summation of the filter products, is connected to the digital output of the ADC. During normal operation, the linearity corrector will provide a corrected output with reduced distortion compared to the output of the ADC. An acquisition memory is connected to the digital output; and a processor for controlling the first signal generator and the second signal generator is connected to read the acquisition memory data. The processor is also coupled to the linearity corrector for programming filter coefficients for the filter product. the

A method of calibrating a linearity corrector is also provided. The calibration method is accomplished by inputting a first waveform and a second waveform to a signal processing system having an output. The acquisition memory and processor acquire and analyze intermodulation and harmonic components of the first and second waveforms as each is output from the signal processing system. By determining the minimum of the objective function based on the distortion model and the intermodulation and harmonic components, the magnitude and phase response of the filter used in the linearity corrector connected to the output of the device is found. This magnitude and phase response is then used to compute a set of filter coefficients. the

Description of drawings

Figure 1 (Prior Art) is a block diagram of an arrangement of a prior art linearity corrector for compensating loudspeakers. the

2 is a block diagram of a filter product based linearity corrector including compensation of 1st order, 2nd order and 3rd order distortion. the

Figure 3 is a block diagram of a filter product based linearity corrector including compensation of 1st and 3rd order distortions. the

4 is a block diagram of a filter product based linearity corrector including compensation of 1st order, 3rd order and 4th order distortion. the

Figure 5 is a block diagram of a compensated ADC. the

Fig. 6 is a block diagram of an ADC system including a compensation and calibration system. the

Figure 7 shows the basic calibration process. the

Figure 8 shows additional sub-steps of the calibration process. the

Figure 9 shows additional sub-steps of the calibration process. the

Figure 10 shows additional sub-steps of the calibration process. the

Detailed ways

This specification includes pages 21-51 of the example source code and its annotation file for illustrating the specific implementation of the specific embodiment of the present invention, and the file name is Slavin.txt. the

As mentioned above, the previously proposed solutions are based on Volterra filters. However, since Volterra filters are very large and difficult to implement with ADCs, a solution is needed that utilizes a more manageable filter design while also reducing some of the remaining major distortions. Taking the Volterra filter as a starting point, a generalized nonlinear filter system can be defined mathematically as:

             (方程1) 

(equation 1)

Wherein, N is the impulse response length of the filter, and k is the order index of the filter. the

For example, if n=3, we have the DC value (h ₀ ), the sum of the linear FIR filter term at k=1, the 2nd order distortion filter at k=2 and the 3rd order filter at k=3 . Thus, for n=3, the Volterra filter can be expressed as:

                 (方程2) 

(equation 2)

The Volterra filter coefficients are linear with the output y, so in theory, it is possible to use the training data to find a set of h. Some of the product terms happen to be permutations of the same set of input samples, so for each order index k the number of distinct values in the set h is given by the following equation, which corresponds to the number of filter taps:

                (方程3) 

(Equation 3)

Unfortunately, this impulse response of a pipelined ADC system can be quite large, so N can be very large, resulting in a large number of taps for k=3. For example, if the impulse response of the pipelined ADC system is 64 clock cycles, so that N=64, then the number of taps required will be approximately 44000. Additional taps are needed for other order filters, if present. the

Existing embodiments of ADC linearity correctors rely on a subset of the Volterra filter system. A subset of this Volterra filter system is characterized by:

                  (方程4) 

(equation 4)

For 1≤k≤n, the system order n defines a set of product orders. the

Although the number and value of taps are unknown, the distortion model is assumed to be of this form. The correction model has the same form, except that the order and filter length are chosen in advance based on experiments with a particular ADC configuration. Calibration involves finding the filter taps. Note that in general the filter taps are different for each filter. For system order n=3, ignoring the h ₀ (DC) term, we get:

(equation 5)

Embodiments of such a structure may be characterized as having a product of N-tap filters implemented using filters for each linear convolution. the

FIG. 2 shows an embodiment of a linearity corrector 100 for implementing Equation 5. Referring to FIG. The output of a signal processing system such as an ADC is provided as input to the linearity corrector 100 . Each of the linear convolutions, as provided in Equation 5, is implemented using filters 102-112. This filter can be implemented as a FIR filter. The first order term corresponds to filter 102 . In an alternative embodiment, this first order term is obtained by replacing filter 102 with a fixed delay approximately equal to half the length of the other filters. In another embodiment, the first order term is obtained by replacing filter 102 with a combination of fixed delays and filters such that the total delay is approximately half the length of the other order filters. The second order term is implemented by multiplying the outputs of filter 104 and filter 106 using multiplier 120 to produce a second order filter product. The third order is achieved by multiplying the outputs of filter 108, filter 110, and filter 112 using multiplier 122 to produce a third order filter product. The output from filter 102 is then added to the outputs of multipliers 120 and 122 using adder 124 to provide a simple sum of the filter products as output. The term "simple sum" as used herein refers to the operation of adding the values of the multiplier without additional filtering between the multiplier and adder 124 . This simplicity is achieved by connecting the multiplier directly to the adder. The term directly connected (or directly connected) is used herein to mean that there are no filters or other processing elements on the path, and there may be registers or other elements in the path that do not alter the essential characteristics of the signal data on the path. The output is now a compensated signal with reduced non-linearity produced by the signal processing system, eg ADC. It should be noted that embodiments of the invention eliminate the filter after the multiplier provided in prior art solutions. While this would require the use of a filter with extra taps compared to the prior art shown in Figure 1, it enables the use of this filter to compensate for some distortion without the need for additional compensation from the filter following the multiplier. For example, if the prior art of Fig. 1 used an all-pass output filter delayed by half a clock cycle (a so-called sm(x)/x or sinc(x) filter), this filter could be incorporated into the filter. By only using a filter before the multiplier, it is possible that the filter product system will be able to better distinguish between original and aliased components. the

While in some embodiments the amount of computation can be reduced by using filters of different lengths, the use of longer length filters increases the number of variables that need to be resolved during calibration, which slows down the calibration algorithm. For hardware implementations, filters of longer length would also require additional delays to match the filtered signal delay. Therefore, in this embodiment of the corrector 100, all filters are of equal length. the

FIG. 3 shows an embodiment of a linearity corrector 100 designed to compensate for 1st and 3rd order distortions without compensating for 2nd order distortions. In some applications, the 2nd order distortion is not significant enough to justify including 2nd order compensation. As shown in FIG. 3 , third order compensation is provided by multiplying the outputs of filter 108 , filter 110 , and filter 112 using multiplier 122 to produce a third order filter product. A simple sum of the 3rd order filter product and the 1st order filter product can then provide a compensated signal with reduced or eliminated 1st order and 3rd order distortion. the

As shown in Figure 4, an embodiment of the corrector 100 including 4th order compensation can be provided. Using the general form of Equation 4, for system order n=4, and neglecting the h ₀ (DC) term, we have:

(equation 6)

As shown in FIG. 4 , this 4th order term may be performed by multiplying filter 140 , filter 142 , filter 144 , and filter 146 together using multiplier 148 . Furthermore, the filter product is directly summed with other order compensations using adder 124 without any intermediate filter between the multiplier 148 and adder 124 . As is clear from the preceding examples, one of ordinary skill in the art will be able to solve Equation 4 to place it in a similar form to Equations 5 and 6 for any desired order, which can then be used as generally taught herein A simple sum of filter products is implemented. As shown in Figure 4, the 2nd order compensation is removed to allow for calibration as discussed in detail below. This configuration will provide compensation in those environments where order 4 is of greater importance. Similarly, compensation for 4th and 5th order can be provided and 2nd and 3rd order eliminated to allow for calibration as described below. the

The justification for the validity of the proposed decomposition of the usual Volterra form, and the corresponding filter product structure, lies in the understanding that for each product order k there is a single self-modulation mechanism, not in the fact that any stochastic Volterra filter system is capable of such on the possibility of the way being decomposed.

Various embodiments of linearity corrector 100 may be implemented using dedicated hardware, such as an FPGA or ASIC, or using a general purpose processor running software. Currently, FPGAs or ASICs are useful for performing real-time corrections, although software running on general-purpose processors is useful for post-acquisition corrections. In the future, it will also be possible to utilize software on a general-purpose processor for real-time corrections. the

Although in some embodiments, the linearity corrector 100 is used to compensate the signal output from the ADC, in other embodiments, the structure of the linearity corrector 100 can be integrated in the same package (packing) as the ADC, or may Integrated on the same chip as this ADC to form a compensating ADC. The compensated ADC 190 is shown in FIG. 5 . It includes an ADC module 192 that contains various circuits that convert analog signals into digital signals. The digital output of the ADC module 192 is input to a linearity corrector 100 which may be implemented with the above teachings. The output of the linearity corrector is a compensated output with reduced distortion. This combined structure provides a compensated ADC. the

In order to properly optimize the linearity corrector described above, it is necessary to calibrate the linearity corrector to determine appropriate filter coefficients for each filter. Unlike the general Volterra filter, the filter product output of the corrector shown in Figure 2-4 is not linear with respect to its coefficients, so in general, the calculation of the filter coefficients is a nonlinear optimization problem. the

Therefore, a linearity corrector system with calibration and a calibration method are proposed. As shown in FIG. 6, an ADC 200 equipped with a linearity corrector 100 and a calibration circuit. An analog select switch 202 is provided to control the input into the ADC. The analog selection switch 202 is controlled by a processor 204 . During calibration, processor 204 controls the analog select switch 202 to allow a calibration frequency signal to enter the ADC 200 input. Processor 204 also controls digital selector switch 206 to control whether the digital output directly from the ADC 200 goes to acquisition memory 208 during calibration. The processor 204 also controls a first calibration frequency generator 210 and a second calibration frequency generator 212 . During calibration, the outputs of the first calibration frequency generator 210 and the second calibration frequency generator 212 are combined in a summing circuit 214 before being input to the ADC 200. In an embodiment of the calibration method, the analog calibration signal and analog summation into the ADC have significantly less distortion than the ADC being calibrated. the

During calibration, the acquisition memory 208 acquires uncorrected ADC output for calibration. Process data acquired multiple times over a range of calibrated frequency pairs. The collected data is processed by processor 204 . the

During normal operation, the analog input is fed to the ADC 200 through the analog select switch 202, and the output of the ADC is corrected using the linearity corrector 100. The basic operation of the corrector has been described above in connection with Figures 2-4. Filters with this linearity corrector can be loaded by processor 204. Other processing can then be provided using the second processor 216, for example the corrected output can be processed like a direct output to the ADC. The embodiment shown in FIG. 6 shows acquisition memory 208 for calibration and other processing. In alternative embodiments, separate memory may be used for other processing. Similarly, the embodiment shown in Figure 6 shows a second processor 216, in other embodiments processor 204 can provide this additional processing and control the calibration process. the

As shown in Figure 7, there are four basic steps to calibrating a bank of filters within a filter product in a linearity corrector. For each order of product above order 1, this step is performed using that order of product as a parameter, eg order 2 or order 3. Step 310 acquires and analyzes the intermodulation (IM) and harmonic components of the input frequency pair. Step 340 uses the results of step 310 to find the magnitude and phase response of the filters in the correction network. Step 360 uses the magnitude and phase responses found in step 340 to design each filter according to the desired magnitude and phase response using a nonlinear phase filter design algorithm. Step 380 programs the filter within the linearity corrector 100 to correct for ADC distortion by setting the filter taps. the

The basic signal acquisition and analysis of step 310 can be accomplished by selecting the base data record length as provided in step 410 in FIG. 8 . In step 420 a set of calibration frequencies is selected. Step 430 generates two sinusoidal inputs based on the nominal calibration frequency, and the calibration index selected in step 420 . Step 440 calculates a set of harmonic and intermodulation frequencies for each calibration center frequency index. Step 450 tests for aliasing that may cause frequency components to overlap. Step 460 either i) performs a set of discrete Fourier transforms (DFT) at exactly the calculated harmonic and intermodulation frequencies, or ii) performs a fast Fourier transform (FFT) that produces results for the complete set of possible frequencies , then select the output at the calculated harmonic and intermodulation frequencies. the

In step 410, a base data record length L is selected for DFT analysis. To improve efficiency, in some embodiments of the calibration method, L is chosen to be a power of two. The choice of L corresponds to the fundamental analysis frequency F=f _s /L, where f _s is the sampling rate. All frequencies f _m are then expressed as harmonics (integer multiples) m of the fundamental analysis frequency F, ie f _m =mF. The angular frequency is expressed as:

ω _m = 2πf _m = 2πmF = 2πmf _s /L (Equation 7)

At step 420, a set of calibration frequencies is selected that will be used in subsequent steps to design the filter response. These frequencies are chosen as harmonic multiples of F and are regularly spaced throughout the Nyquist band where the input frequency can occur, but always excluding the extreme of 0 (DC) and f _s /2.

At step 430, for each calibration index c with nominal calibration frequency ω _c , two sinusoidal inputs are generated at angular frequencies ω _c-1 = ω _c - ω ₁ and ω _c+1 = ω _c + ω ₁ corresponding to At:

x _c (t) = p _c-1 sin(ω _c-1 (t+λ _c-1 ))+Q _c+1 sin(ω _c+1 (t+λ _c+1 )) (Equation 8)

where λ _c is the delay at frequency index c. In an embodiment of the calibration method, the magnitudes {P _c-1 , Q _c+1 } are similar to improve IM calibration results. Any small differences can be removed in a later correction stage. In an embodiment of the calibration method, the calibration frequencies are chosen such that they use a large part of the ADC input dynamic range. Then, an integer multiple of M groups of L adjacent A/D output samples is collected. Every Mth sample is averaged to obtain low noise y(t), t=0...L-1, where L is the base data record length provided in step 410 above. Simple averaging of this type is very fast, allowing fast processing of large numbers of values (could be millions of samples or more) for ML. L should be chosen such that it is not too small or that the analysis frequencies would be too far apart. Large values of L can generate more calibration frequency data than needed, causing subsequent algorithms to run more slowly. A value of L=2 ¹⁰ =1024 has been successfully used.

At step 440, for each calibration center frequency index c (where 0<c<L), a set of harmonic and intermodulation frequencies is generated. The DC term is ignored because DC errors can occur for many other reasons. For example, a second-order system with IM and harmonic distortion frequencies at {ω _c-1 −ω _c+1 , 2ω _c-1 , ω _c-1 +ω _c+1 , 2ω _c+1 ,}, corresponding to the exponent : m={2, 2c-2, 2c, 2c+2}. In fact, it is possible to calibrate using only the sum-frequency term m={2c-2, 2c, 2c+2}. The advantage of using only the sum term of frequencies and not the difference term is that the former is not generated by any low order distortion filter products. This includes the original input signal of what may be referred to as a 1st order (linear) "distorted" system.

It should be noted that in sampled systems, IM and harmonic distortion will contain aliased frequency components. Therefore, at step 450, a test is performed to handle the case where aliasing results in overlapping frequency components, since the measurement will not agree with the distortion model. For a given c, generate a set of possible harmonic and IM terms (except DC). In some embodiments, this set will contain all possible harmonic and IM terms except the DC term. For each frequency of indices _mi from list m in step 440, since sampling theory cannot distinguish aliasing components with frequency interval L, they are convoluted with each other as follows:

w(m _i )=m _i mod L (equation 9)

Also, within this wrapping interval L, aliasing cannot be distinguished between w(m _i ) and Lw(m _i ). This aliased index remapping is the lower range interpretation, given by:

alias(m _i )=min(w(m _i ), Lw(m _i )) (equation 10)

where min() is the minimum value of its input value. Therefore 0≤alias(m)<L/2. This remapping of m _i to alias(m _i ) is only used to detect any distortion frequency collision for a given integer m _i and L. Collision detection is mostly easy to implement by generating all harmonic and IM frequency indices, then remapping each frequency index as described above in Equation 9 and Equation 10, and numerically sorting the resulting aliased frequency index values. Then, adjacent frequency indices with the same remapped frequency index value are found. If a conflict is found, no measurement is performed at this value of c. If L is a power of 2, and c is even, then {c-1, c+1} is always relatively prime. Thus, collisions only occur at aliased frequencies. Even then, conflicts are usually rare. For example, in a 3rd order system, collisions only occur when c=L/4. In this case, the intermodulation product is located at 2(L/4-1)+(L/4+1)=3L/4-1, so its aliasing frequency is at L-(3L/4-1)= L/4+1, this frequency is also the input frequency. Similarly, (L/4-1)+(2L/4+1)=3L/4+1, which aliases onto /4-1, is also the input frequency.

In one embodiment of the calibration method, a real input DFT of a set of L-samples is performed at the selected harmonic and intermodulation frequencies, as shown in step 460 . The real and imaginary components at each frequency m are obtained using a single frequency DFT on this averaged real sample sequence:

X _m = Re(Z _m )

Y _m = Im(Z _m ) (Equation 11)

where Re and Im are the real and imaginary parts of the complex input. The form of Equation 11 is for the normalized sampling frequency from Equation 7 provided in connection with step 410 above. In an alternative embodiment, a Fast Fourier Transform (FFT) may be used instead of multiple DFTs to obtain real and imaginary components at each frequency m. the

In this distorted system, each filter has an unknown magnitude and phase response at this stage. In order to find these responses, a mathematical model of the distortion for the system is provided. Based on this model provided, nonlinear optimization can be performed. The self-modulation model for the Nth order product represents n unknown filter responses as modified inputs for each product term:

(Equation 12)

This self-modulation model can be rewritten in an optional general form as:

(Equation 13)

in:

                                  (方程14) 

(Equation 14)

and

B _i,c+1 = A _i,c+1 /A _i,c-1 (Equation 15)

Assuming a sufficiently large L, if the two pairs of frequencies C-1 and C+1 used for calibration are close enough, then for a given filter, the delay at each frequency can be assumed to be approximately the same, then Equation 15 can be simplified to :

             （方程16) 

(Equation 16)

For this nth order distortion system, the number of variables is equal to 2n+1. The sine wave input amplitudes to the ADC are closely matched for each pair of calibration frequencies used. Otherwise, the system cannot distinguish the difference between their magnitudes at two different frequencies and the difference in the magnitude response of each filter. the

If the magnitude response through each filter is approximately the same at these two different calibration frequencies, then a more simplified self-modulation model can be given as follows:

              (方程17) 

(Equation 17)

By choosing an ADC that fits one of these distortion models, as expressed in Equation 13, Equation 16, or Equation 17, a set of phase and magnitude solutions can be found for each filter channel at each frequency. Any equations representing this distortion model can be extended to provide a set of intermodulation and harmonic frequency components corresponding to those found at step 440 above. This reduction can be expressed as cosine (real part) and sine (imaginary part) terms at each frequency. Once the expansion and collection of frequency terms is obtained, along with the set of DFTs determined at step 460, the difference between the real and imaginary parts of these expressions and the DFT result can be obtained symbolically The sum of squares, which can be obtained from the DFT or FFT as described above. This results in a large equation for the 2nd order distortion, and an even larger equation for each higher order. the

because of phase And the delay λ is related to the angular frequency ω, that is, by:

(Equation 18)

Substituting ω _m from Equation 7, and using normalized frequencies (ie, f _s =1), rewrites any Equation 13, Equation 16, or Equation 17 in the model equations in frequency exponential form. For example, Equation 16 can be rewritten as:

         (方程19) 

(Equation 19)

For example decomposing into real and imaginary parts with order n=2, we extend Equation 19 and collect terms to obtain terms at four non-DC intermodulation and harmonic frequency indices: {2, 2c, 2c-2, 2c+2 }:

(Equation 20)

            (方程21) 

(Equation 21)

(Equation 22)

(Equation 23)

where the set of coefficients {C, S} is determined as:

For equation 20:

For equation 21:

For equation 22:

For equation 23:

Now, estimate the sum of squares of the differences between these coefficients and their corresponding DFT or FFT components from Equation 11:

(Equation 24)

Here Equation 24 may represent the objective function or error. For the second-order case, the unknown quantity to be searched can be expressed as a vector:

x _c ={K _c , B _1,c+1 , B _2,c+1 ,λ _1,c ,λ _2,c } (Equation 25)

Once this objective function has been provided, the magnitude and phase response of each filter can be determined. Referring now to FIG. 9, one embodiment of the present calibration method is shown. As shown in step 540, the objective function is first minimized. After minimizing the objective function, a post-solution correction is performed on the magnitude, as shown in step 560 . The corrected solution is now decomposed between the filters provided for each stage, as shown in step 580 . the

To minimize this objective function, as shown in step 540, will require solving a system of nonlinear equations. A class of Newton-based algorithms can be used to solve these nonlinear equations and minimize the objective function for each desired order. Possible Newton-based algorithms include Modified Newton (MN), Gauss-Newton (GN) and Quasi-Newton (QN). At each search iteration k, these methods estimate the numerical gradient vector (g _k ) of the objective function and the numerical inverse of the Hessian (S _k ) of the objective function, or for the GN and QN methods the Hessian of the objective function inverse approximation. The gradient vector of the objective function and the inverse of the Hessian of the objective function can be described as follows:

    在x＝x_k                 (方程26) 

at x = x _k (equation 26)

At x = x _k (equation 27)

By estimating the gradient vector and the inverse Hessian, the search direction (d _k ) that can be used to minimize the objective function can be obtained, where the search direction is given by:

d _k = -S _k g _k (Equation 28)

的梯度向量定义为这样一个向量，其第i个分量是该函数关于其第i个未知变量的偏导数。因此，它是这个阶段的符号结果而非数值。对于对应于方程16的具有5个未知量的第二阶情况，该目标函数由方程24提供：  The objective function

The gradient vector of is defined as a vector whose ith component is the partial derivative of the function with respect to its ith unknown variable. Therefore, it is a symbolic result of this stage rather than a numerical value. For the second order case with 5 unknowns corresponding to Equation 16, this objective function is given by Equation 24:

The (i, j)th element of the Hessian matrix is the partial derivative of the objective function with respect to its i-th and j-th variables. Again, for the 2nd order case, the Hessian matrix is given as:

The matrix elements are therefore entries for the unknown variable x. The partial derivative operator is commutative, so the Hessian matrix is always square-symmetric. the

The corresponding gradient and Hessian matrix can be found for this objective function. Although solving the derivatives of such complex expressions is very tedious, computerized symbolic computations can be performed by using symbolic mathematics software tools such as the symbolic derivative function in Mathematica. Alternatively, the derivative f′(x) of any function f(x) can be solved numerically using well-known approximations:

(Equation 29)

For suitably small values of Δ, this computational precision is assumed to be available. Symbolic or numerical differentiation can be used. Symbolic differentiation may be provided by symbolic mathematics software, and such processing provided by such software uses a general-purpose processor for processor 204 . It should be noted that if symbolic differentiation has been obtained, in some embodiments of the present calibration method these symbolic results can be encoded as functions using a lower level language such as C to improve the performance of this operation in processor 204 . In an alternative embodiment, numerical differentiation may be used by the processor instead of an encoded implementation providing the symbolic solution. the

As mentioned above, the inverse Hessian can be solved directly or iteratively by using Newton-based algorithms. For example, in one embodiment of the present calibration method, an approximation of the inverse of the Hessian matrix is determined iteratively using a quasi-Newton algorithm. There are multiple subcategories of QN iterations that can be used. For example, the BFGS (Broyden, Fletcher, Goldfarb, Shanno) iteration formula can be used. In one embodiment of the calibration method, the BFGS iterative formula is used with the Fletcher's Line Search algorithm. The (k+1)th BFGS iteration of the inverse Hessian S _k+1 is given by:

(Equation 30)

where the column vector p _k , the scalar q _k and the square matrix r _k are given by:

p _k =S _k γ _k

(Equation 31)

S _k (square matrix) is the kth approximation of the inverse Hessian, and S ₀ =I is the identity matrix. It should be noted that the superscript T is transpose, eg: is the transpose of δ _k . and a column vector:

γ _k =g _k+1 -g _k (Equation 32)

where g _k is the gradient vector obtained earlier from (Equation 26). The new solution vector estimate xk ₊₁ for the next iteration is given by the correction column vector _δk applied to the previous _xk :

x _k+1 = x _k +δ _k (Equation 33)

To find the unknown _δk for each iteration, a constrained search along the search direction _dk is used. Once x _k is found using Equation 33 and d _k is found using Equations 28 and 30, then, a non-negative scalar α _k is found that minimizes:

e ₂ (x _k+1 )=e ₂ (x _k +δ _k )=e ₂ (x _k +α _k d _k ) (Equation 34)

This is a one-dimensional optimization problem in α _k regardless of the dimensionality of the original problem. Thus, the solution of this particular type of equation is known as a line search.

The current solution may estimate a vector x _k , and its corresponding search direction d _k is given as input to the Fletcher line search algorithm to give a scalar value α _k . Then, the column vector correction for x _k in Equation 33 is given by:

δ _k =x _k+1 -x _k =α _k d _k (Equation 35)

It is well known that finding the global minimum of the objective function is a difficult problem in general. The minimum found always depends on the objective function, but can also depend on the initial condition x ₀ ° If the distortion model fits reasonably well, the convergence is usually fast and reliable. However, if the true distortion mechanism includes multiple distortion products or some other mechanism that does not follow the distortion model, the solution will diverge. In this case, a solution is advanced, and the error function equation 24 decreases in each iteration, but one of the unknowns in the vector x diverges exponentially towards infinity. A limit on the number of BFGS iterations can reliably detect this kind of problem. However, reaching this limit also slows down the calibration algorithm.

If, after noise reduction averaging, divergence is seen, a divergence failure flag can be set to indicate the need for further diagnostics, or failure of the ADC. the

However, once the minimum solution is approached, the algorithm converges quadratically, so it converges quickly to repeatable results. This property can be used to explore the solution space of the model of Equation 16. A solution can be rejected as a duplicate if it converges to a previously found solution. This process continues until a statistically satisfactory decision has been obtained that no further minimum can be found. the

In some embodiments of the present calibration method, there may be other adjacent solutions at angles having adjacent integer multiples of π. These duplicate solutions were identified and rejected as trivially relevant or "uninteresting". the

One way to find the solution is to start x with a random value of B in Equation 19, using:

B = e ^{u(Random-1/2)} (Equation 36)

可以在以下范围中随机产生：  Random is a random real value from 0 to 1. A value of u=0.1 can have good results for real ADC data. each corner

Can be randomly generated in the following ranges:

                  (方程37) 

(Equation 37)

The angle can then be converted to a delay in Equation 19 at a known frequency using Equation 18. the

In another embodiment, a solution for a known ADC design is found by restricting the initial range of random angles used as a function of the calibration index c to the range that worked in previous calibrations. This restriction on the initial range speeds up the search for the solution. the

Even if the starting value is constrained by at least Equation 37, the final value may exist outside these ranges into the closest adjacent iterative solution. This is because the start point may exist on the other side of the repeated "peak" between the repeated minima. However, due to the similarity of the repeated function "landscapes" at the interval π, no repeated solution will be found even further. In one embodiment of the present calibration method, nearest neighbor solutions are found and duplicate solutions are excluded. the

For each calibration frequency, once the first solution angle is found, the central solution (i.e. the angle closest to 0) can be easily found. In one embodiment of the present calibration method, the found angles are directly verified as central solutions (within ±π/2). In an alternative embodiment, multiple=floor(angle/π+0.5) is used to find the closest integer multiple of π to the central solution region, and then the multiple of π is subtracted from the first solution angle (solution angle) Go to bring the corner into the center area. In some applications, a solution close to ±π/2 may fluctuate with noise between +π/2 or -π/2 in an otherwise smooth change in angle with a calibration frequency. These fluctuations can be corrected by looking for angular continuity in later processing. Then, a second run of the BFGS algorithm can quickly find a more accurate central solution by using the substituted first solution as a good starting point.

From this central solution, the nearest neighboring solutions to this central solution can be found systematically. As in the process of finding the central solution above, the repeated solution can be found by adding all combinations of {-π, 0, π} to each corner and using the resulting vector as the starting point for the BFGS solution search. Again, because each starting point is very close to the correct value, the BFGS algorithm converges very quickly. There are 3 ² −1=8 adjacent solutions for the 2nd order search and 3 ³ −1=26 adjacent solutions for the 3rd order search.

After finding these solutions, a search at other random starting positions can be used to find any more solutions, if they exist. If a solution matches one of the previously found solutions, it is rejected. As mentioned earlier, the process of verifying duplicates is made easier by the high repeatability of the solutions found by the BFGS algorithm with signed gradients. the

If a new solution does not match any found so far, its central solution can be found, and all its neighbors can be found using the same method as before. In one embodiment of the calibration method, fruitful search results reset a timeout counter (tmeout counter), while failures to find new solutions decrement the timeout counter. If this timeout reaches zero, the search for other solutions stops. The central solution that gives the best error measure according to Equation 24 is chosen as the final solution for the calibrated frequency index. the

In one embodiment of the present calibration method, if a solution has been found at the previous calibration frequency index c _prev , it can be used as the starting point for the next search at c. If it fails to converge, the aforementioned slow random search method is used.

If no solution is found for a particular c, a second pass can be made using the solution at c _next as a starting point. If a solution is still not found, then patching or interpolation between adjacent values can be used. If no solution can be found for too many calibration value indices, then the ADC is faulty and a fault flag can be set.

If an acceptable set of calibration frequencies has been found for all calibration frequencies, then the magnitude and phase response of each channel into the product of the self-modulation system is characterized by the set of selected discrete calibration frequencies. With a uniform distribution of these frequencies, the response for each channel can be represented by an FIR filter as in Figures 2-4. In one embodiment of the calibration method, the coefficients of each filter are now iteratively found. the

In another embodiment of the calibration method, the set of vectors _xc from Equation 25 are calibrated before they are passed to the filter design, as shown in step 560 . In an embodiment of the calibration method, the magnitude of the calibration input will not significantly affect the filter response over a reasonable calibration range. As can be seen from Equation 14, this dependence is achieved by using can be eliminated by dividing _Kc , where n is the order index of the product. Also, in some embodiments, any difference in the amplitudes of the two calibration sinusoidal inputs of Equation 8 should also have negligible effect on the filter response. In this case, not only the relative response of the distortion filter, but also the relative magnitude of the calibration sinusoidal input affects the measurement. Therefore, in these embodiments, the relative magnitudes P _c-1 and Q _c+1 are used for correction. The ratio:

R _c =Q _c+1 /P _c-1 (Equation 38)

gives the amount by which the B parameter is overestimated in x. This parameter of x in equation 25 is therefore corrected for the nth order calibration to give a new vector X _as :

            (方程39) 

(Equation 39)

At any calibration frequency, the combined magnitude response of the filter product system is a function of the magnitude response of each filter. Thus, we can make the coefficients of one filter twice as large and the other twice smaller, with the same result. Variations in the magnitude response of the system can be significant between calibration frequencies. We can achieve this across all filters in the filter product system by making each filter in the product have the same 1/nth root response from the first term in Equation 39 for the nth order distortion corrected system The average dispersion tracks the magnitude response to varying loads such that:

                    (方程40) 

(equation 40)

For an nth order distorted system, there are n pairs of amplitude and frequency responses in each vector _Xc of frequency index c. Arbitrarily specifying a pair of magnitude and phase responses to one of the n filters may produce a non-smooth magnitude response, phase response, or both for each filter as a function of frequency. A non-smooth response can result in a filter design that does not track the expected response well for the chosen filter length. Thus, the magnitude and phase responses are distributed among the available filters at each stage, as shown in step 580 . Because a sufficiently dense set of calibration frequencies usually results in a high degree of correlation between the delayed responses to adjacent calibration indices, in one embodiment of the calibration method, by using the calibration cycle period (1/f _c- Modulo half of ₁ ), test possible permutations of the new set relative to the previous set at frequency index c _prev to obtain a better fit between this new set of delays and the previous set of delays. In case of an odd number of slips, a delay corresponding to half a cycle is applied to one of the filters. In embodiments of the calibration method, such as those employing the model described in Equation 16, each phase value has a corresponding magnitude, so rearranging the phases also rearranges the corresponding magnitudes. If two phases cross, they can be grouped into the minimum phase value for one filter and the maximum phase value for the other filter. The resulting kink on each phase map is not as detrimental to the resulting filter quality as a discontinuity, but results in a feasible filter design, provided enough filter taps are provided.

Through the previous calibration steps, an appropriately grouped set of magnitude and phase responses is obtained for each filter path. The filter design step 360 in FIG. 6 designs each filter using the magnitude and phase responses calculated in step 340 to give an approximation of the calculated response. In the case of radio frequency (RF) applications, the middle part of the upper Nyquist band of the ADC is used without calibration for the outer frequency regions in the band. This means that the external frequency is not constrained. As shown in Figure 10, in an embodiment of the calibration method, an iterative process is used to accomplish the filter design and find the filter taps. In an alternative embodiment, a simple inverse FFT can be used. the

As shown in step 610, a filter length is selected for each channel. In one embodiment of the calibration method, the same filter length is used for all channels. In alternative embodiments, different filter lengths are used for some or all channels. In this alternative embodiment, an additional delay may be required to compensate for this different filter length. The filter length can depend at least in part on the expected performance and the ADC being calibrated. By way of example, a full filter length of 63 taps or a full filter length of 127 taps are used. the

At step 620, an FFT length is selected and the magnitude and phase responses are converted to a filter frequency response vector in complex Cartesian coordinates. The FFT length should be longer than the filter length selected in step 610 . For example, in an embodiment of the present calibration method, a real FFT length of 512 samples is chosen. Because the FFT is longer than the filter length, values other than frequencies corresponding to filter taps are initialized to zero here. the

At step 630, an inverse FFT of the filter frequency response vector is performed to obtain a set of real filter taps, but also taps extended to the FFT length, eg 512. In the time domain, taps beyond this chosen filter length are forced to zero. the

At step 640, a forward FFT is performed to go back to the frequency domain. Some regions of the frequency response vector should be unconstrained. The FFT bins at the frequency of interest are rewritten back to the original values for the expected magnitude and frequency response, leaving the unconstrained frequency bins unchanged. the

Steps 630 and 640 are repeated to iteratively converge to the desired solution. Once the filter tap provided at step 630 has become stable from the previous iteration, the filter tap is found. the

Once the filter taps are obtained from step 630, the hardware can be programmed as provided at step 380 in FIG. The programmed filter product system is now capable of correcting for ADC harmonic and intermodulation distortion. the

It will be apparent to those of ordinary skill in the art that many changes may be made to the details of the above-described embodiments of the present invention without departing from the basic principles thereof. Accordingly, the scope of the invention should be determined by the appended claims.

Slavin.txt

Mathematica code for ADC Correction Algorithm. 

(＊Copyright 2005 Tektronix Inc.＊)

(＊by Keith Slavin＊)

(*COPYRIGHT NOTIFICATION

＊Pursuant to 37 C.F.R.1.71(e), applicant notes that a portion of this disclosure

contains material that

＊is subject to and for which is claimed copyright protection, such as, but not

limited to, source code listings,

*screen shots, user interfaces, or user instructions, or any other aspects of this

submission for which copyright

＊protection is or may be available in any jurisdiction. The copyright owner has no

objection to the facsimile

＊reproduction by anyone of the patent document or patent disclosure, as it appears

in the Patent and Trademark office

＊patent file or records. All other rights are reserved, and all other

reproduction, distribution,

＊creation of derivative works based on the contents, public display, and public

performance of the application

＊or any part thereof are prohibited by applicable copyright law. 

*)

BeginPackage["correction`", "Global`"]

(*frequency Indexing: i1=i-1, i2=i+1*)

(*

＊design filters to correct ADC output.lowF，hiF are normalized

＊frequencies from 0 to 0.5(of fs).'data'is a 2-Darray of acquired

＊ADC data in/login/keith/adc_correction/xcal_data6/summary.To load，

＊type:<<Z:adc_correction\xcal_data6\summary; The result is in "data". 

＊Type: ProcessData[data, 0.1, 0.4, 127, 10^-12] to get the filter

＊coefficients in the form: 

*{{2nd_orde_f1_list, 2nd_order_f2_list}, 

＊{3rd_order_f1_list, 3rd_orde_f2_list, 3rd_order_f3_list}} 

*)

ProcessData[caldata_, calFreqLo_, calFreqHi_, correctFilterLen_, error_]: =

Module[

{order, responses2, correctionFilters2, responses3, correctionFilters3,

correctionFilters},

order=2;

Print[″Generate′rsponses2′:″];

responses2=ProcessCalData[ca]data, order, calFreqLo, calFreqHi, error];

(＊design the correction filters from the cal responses＊)

print["Generate'correctionFilters2':"];

correctionFilters2 = DesignFilters[responses2, correctFilterLen];

order=3;

Print[″Generate′responses3′:″];

responses3 = ProcessCalData[caldata, order, calFreqLo, calFreqHi, error];

Print[″Generate′correctionFilters3′:″];

correctionFilters3 = DesignFilters[responses3, correctFilterLen];

Return[{correctionFilters2, correctionFilters3}];

]

Page1

Slavin.txt

(*

＊ Run a calibration on randomly generated distortion filters, and then correct

＊ the same distortion on a multi-sine input signal. 

＊ (＊Note: 0.2...0.3 range puts 3rd order IM in 0.1...0.4fs

＊ and 2nd order IM from 0...0.1 and from 0.4...0.5fs.*) 

* testFreqLo = 0.2;

* testFreqHi=0.3;

＊ testFreqStep＝0.004; (＊step in fs units from testFreqLo to testFreqHi＊)

＊ calFreqLo＝0.1; (＊low calibration range in fs units＊)

＊ calFreqHi＝0.4; (＊high calibration rage in fs units＊)

＊ distortFilterLen＝63; (＊odd length of all distorting filters＊)

＊ spread＝4.0; (＊spread(samples) of random distorting filters＊)

＊ analysisLen＝1024; (＊samples in DFT analysis of cal frequencies＊)

* freqPairs=128; (*frequency pairs in cal(must divide analysisLen)*)

＊ correctFilterLen＝63; (＊correction filter lengths(must be odd)＊)

＊ testLen＝32768; (＊samples used in test signal analysis＊)

*

＊ these three parameters are in LSB units.They make

＊ noise+distortion+original tones with similar amplitudes to the

＊ AD6645 ADC:

*

* distortionAmpl = 0.04;

* tonePairAmpl = 15000.0;

* noiseAmpl = 2.0;

*)

Demo[]:＝

Module[

{} 

MainTest[0.2, 0.3, 0.004, 0.1, 0.4, 63, 4.0, 1024, 128, 63, 32768,

0.04, 15000.0, 0.01];

]

MainTest[testFreqLo_, testFreqHi_, testFreqStep_, calFreqLo_, calFreqHi_,

distortionFilterLen_, spread_, analysisLen_, freqPairs_,

correctFilterLen_, testLen_,

DistortionAmpl_, tonePairAmpl_, noiseAmpl_]:＝

Module[

{margin, testdata, order,

filter2, filter3, filters,

caldata, error,

responses2, responses3,

correctionFilters2, correctionFilters3, correctionFilters,

distorted, corrected},

Print[″Version0.3″];

(*generate second-order random distortion filters*)

print["Generate'filter2':"];

order=2:

filter2 = GenFilters[distortFilterLen, order, spread];

(*generate third-order random distortion filters*)

Print["Generate'filter3':"];

order=3;

filter3 = GenFilters[distortFilterLen, order, spread];

(*combine filter data*)

filters = {filter2, filter3};

(*generate multi-tone test data record*)

Page2

slavin.txt

Print[″Generate′testdata′:″];

margin=(distortFilterLen-1)+(correctFilterLen-1);

testdata=GenTest[testLen, mrgin,

... testFreqLo, testFreqHi, testFreqStep];

(*generate a set of cal data tone pairs*)

print["Generate'caldata':"];

cal data = DistortCal[filters, analysisLen, freqPairs,

distortionAmpl, tonePairAmpl, noiseAmpl];

(*process the cal data to get filter responses*)

error=10^-12;

correctionFilters=

ProcessData[caldata, calFreqLo, calFreqHi, correctFilterLen, error];

(＊distort the test data using the original filters＊)

print["Generate'distorted':"];

distorted=DistortData[testdata, filters,

distortionAmpl, tonePairAmpl, noiseAmpl];

(*plot the distorted signal*)

PlotunwindowedFFT[Take[distorted, testLen]];

Print["Generate'corrected':"];

corrected = CorrectData[distorted, correctionFilters];

(*plot the corrected signal*)

PlotUnwindowedFFT[Take[corrected, testLen]];

]

(*aliases data set by first finding data mod n, and then

wrapping data>n/2 to n-data＊)

AliasF[data_, n_]:＝

Module[

{len, modData, i},

len=Length[data];

modData = Mod[data,n];

Return[Table[Min[modData[[i]], n-modData[[i]], {i, 1, len}]];

]

(＊Expand distortion product＊)

AllTerms[i_, order_] :=

Module[

{t, tt, j, q0, q1, q2, q3, q3a, q4, q5, q6, q7, q8, terms, len},

If[order=0,

Return[{0}];

];

q0＝TrigReduce[(Cos[(i-1)tt]+Cos[(i+1)tt])^order];

(＊extract tt variable outside as a common factor for each cos term＊)

q1=q0/.Cos[a_]:>Cos[Apart[a]];

(＊expand into separate terms＊)

q2 = Expand[q1];

q3 = Table[q2[[j]], {j, 1, Length[q2]}];

(*find the cosine terms only, discard their coefficients, and convert

any DC term to 0＊)

q4=Table[If[NumericQ[q3[[j]]], 0, q3[[j]]/(q3[[j]]/.tt->0)],

      {j, 1, Length[q3]}];

(*eliminate the Cos[]around each term*)

q6=q4/.Cos[a_]:>Apart[a];

(*form(k*i+m)*tt->k*i+m for integer k, m*)

q7=q6/.a_*tt:>a;

Page3

Slavin.txt

q8=q7;

(＊make i term positive:if k is -ve，then its derivative is -ve＊)

terms=Table[If[D[q8[[j]], i]<0, -q8[[j]], q8[[j]]], {j, 1, Length[q8]}];

Return[Expand[terms]];

]

(*

* * *

＊S is an approximation to the inverse of the Hessian, G.

*

＊BFGS iteration is:

＊ T T T T

＊ ga S ga dx dx dx ga S+S ga dx

＊ S[k+1]＝S+(1+-------)------ -(------------------) 

＊ T T T T

＊ dx ga dx ga dx ga

＊ T T T T T T T T

＊ (dx ga+ga S ga)dx dx dx ga S+S ga dx

＊ ＝S+----------------------- -(------------------) 

＊ T T 2 T T

*

dx ga

*

＊ T T T

＊ If: v1＝s ga, v2＝dx ga, v3＝v1 dx ＝S ga dx

＊ then:

＊ T T T

* (v2+ga v1)dx dx (v3+v3)

＊ S[k+1]＝S+----------------- - -------- 

* 2

* (v2) v2

*

＊Note: SetPrecision[] is used to compensate for the gradual loss in

＊accuracy that Mathematica imposes on iterative calculations. 

*)

BFGS[x0_,vars_,gradients_,objectiveFunc_,epsi_,maxIterations_]:=

Module[

{len, k, x, S, g, err, ds, a, ax, dx, xnew, g1, ga, v1, v2, v3, part0, part1, part2,

eTable, rTable, dxdxT, max, subst, residual},

len=Length[x0];

k=0;

(*convert 1-D list to a 2-D column vector*)

x = Transpose[SetPrecision[{x0}, 39]];

(*S[k] is an approximation to the inverse of the modified Hessian of

e2[x]at x＝x[k]＊)

S＝IdentityMatrix[len]; (＊S is a square matrix and always symmetric＊)

g = GradFunc[x, vars, gradients]; (*returns a gradient column vector*)

err=Norm[g]; (*returns sqrt of sum of squares of elements of g*)

eTable = Table[0, {maxIterations}];

rTable = Table[-40.0, {maxIterations}];

While[(err>epsi)&&(k<maxIterations),

   ds＝-S.g; (＊ds is a direction column vector＊) 

(*scalar'a'line search result*)

ax = FletcherLineSearch[x, ds, vars, gradients, objectiveFunc];

a = SetPrecision[ax,39];

dx＝a＊ds; (＊dx is a column vector＊)

  xnew＝x+dx; (＊new column vector result＊)

g1 = GradFunc[xnew, vars, gradients]; (*column vector result*)

    ga = g1-g; (*also a column vector*)

(*dx.Transpose[dx] is a square symmetric matrix*)

dxdxT = dx.Transpose[dx];

v1=s.ga; (*also a column vector*)

Page 4

Slavin.txt

v2[=(Ttanspose[ga].dx)[[1]][[1]]; (*v2 is scalar*)

If[v2==0.0,

    Print[″BFGS error:v2=″, v2]; 

    Print[″Iteration:″,k]; 

    Print[″v2=″, v2]; 

    Print[″g=″,g]; 

    Print[″S=″,S]; 

Print[″Determinant[S]=″, Det[S]];

    Print[″Inverse[S]=″, Inverse[S]]; 

    Print[″ds=″,ds]; 

    Print[″a=″, a]; 

    Print[″dx=″,dx]; 

    Print[″xnew=″,xnew]; 

    Print[″g1=″,g1]; 

    Print[″ga=″,ga]; 

Abort[];

];

v3＝v1.Transpose[dx]; (*v3 is a square matrix＊)

(*Transpose[ga].v1 is a scalar, so part0 is scalar*)

Part0=(v2+(Transpose[ga].v1)[[1]][[1]])/(v2^2);

  part1＝par0＊dxdxT; (＊part1 is square symmetric＊)

(＊any matrix added to its transpose is always square symmetric＊)

Part2=(v3+Transpose[v3])/v2; (*part2 is square symmetric*)

(＊S(old and new)is square symmetric＊)

S = SetPrecision[S+part1-part2, 39];

x = SetPrecision[xnew, 39];

g = SetPrecision[g1, 39];

Subst=-Table[vars[[j]]->x[[j]][[1]],{j,1,len}];

residual = Log[objectiveFunc/.subst];

err = Norm[dx];

(*Print[″iteration=″,k,″,dx=″,err];*)

(*Print[″x=″, Flatten[x]];*)

k+=1;

eTable[[k]]=err;

rTable[[k]]=residual;

If[residual<-40.0,

Break[];

];

];

If [k==maxIterations,

(*Print[″BFGS Failed to converge.Final vector:″,x];*)

Return[];

];

(*Print[″Total iterations=″, k,″, Line error=″, err];*)

(*Print[″Finished:x=″, Flatten[x]];*)

Return[Flatten[x]];

]

(＊＊)

CorrectData[data_, filters_]:=

Module[

{maxOrder, len, c, sum, j, nrSamples, samples, t},

maxOrder = Length[filters];

len=Length[filters[[1,1]]];

c = (len-1)/2;

sum=0;

For[j=1, j<=maxorder, ++j,

  sum + = GenDistortion[data, filters[[j]]];

];

Page5

Slavin.txt

nrSamples = Length[sum];

Print[″nrSamples=″, nrSamples];

samples = Take[data, {c+1, nrSamples+c}]-sum;

Return[samples];

]

(*make real tap filter have the required DC'delay'(in sample units,

where 0 is the first tap or sample.Usually DCdelay＝Length[tans]/2

is used*)

DCadjust[realtaps_, DCdelay_]:＝

Module[

{newtaps, taps, i, sum1, sum2, offset, div},

taps = Length[realtaps];

sum1 = Sum[irealtaps[[i+1]], {i, 1, taps-1}];

sum2 = Sum[realtaps[[i+1]], {i, 0, taps-1}];

div = taps * (DCdelay-(taps-1)/2);

If[div==0.0,

Print[″DCadjust divide by zero.″];

Abort[];

];

offset＝(sum1-sum2＊DCdelay)/div;

Print[″offset=″, offset];

newtap=realtaps+offset;

Return[newtaps];

]

(*find the'delay'of a filter at DC*)

DCdelay[realtaps_]:＝

Module[

{sumc, sumci, i, taps},

(*calculate the dc delay without Arg[]aliasing.*)

taps = Length[realtaps];

  sumc = Sum[realtaps[[i+1]], {i, 0, taps-1}];

If[sumc==0.0,

Print[″DCdelay divide by zero.″];

Abort[];

];

sumci=Sum[realtaps[[i+1]]i, {i, 1, taps-1}];

Return[N[sumci/sumc]];

]

Delays[list_,n_]:＝

Module[

   {i, L, len, reals, angle, diff, offset, items, complex,

delay, radiansPerSample, prevAngle, newDelay},

L = Length[list];

len=n/2;

Items=len+1;

Reals=Join[list, Table[0, {i, 1, n-L}]];

complex = Take[Fourier[reals, {FourierParameters->{-1, 1}}], items];

delays = Table[0, {items}];

Delays[[1]]=DCdelay[list];

For[i=1, i<items, ++i,

If[Abs[complex[[i+1]]]<10^-12,

Delays[[i+1]]=delays[[i]],

(*else*)

` radiansPerSample=N[Pi*i/len];

delay = delays[[i]];

     prevAngle = delay * radiansPerSample;

angle = Mod[Arg[complex[[i+1]]], pi, prevAngle-Pi/2];

Page6

Slavin.txt

newDelay=angle/radiansPerSample;

Delays[[i+1]]=newDelay;

];

];

(* possible oversampling to interpolate*)

    For[i=1, i<len, ++i, 

If[Abs[complex[[i+1]]]<10^-12,

Delays[[i+1]]＝(delays[[i]]′+delays[[i+2]])/2;

];

];

Return[delays-(L-1)/2];

]

(*

＊ design filters to correct ADC output. 

* DesignFilters[] works on the output of ProcessCalData

*)

DesignFilters[data_, correctionFilterLen_]:=

Module[

{fi, samples, freqs, start, end, skip, i, j, freqIndex, ampl, temp, cycle,

  refDlyA, refDlyB, dly, qq, K, logsize, logFreqs, mags, parts, filters,

    nrVars, gainA, gainB, delays, totalGains, xx, sortList, sorted

 prev, next, permute, undef, modified, odd},

xx=data;

freqs＝Length[xx]-1; (*input is oflength 2＾m+1*)

  logFreqs = Log[2, freqs];

If [! IntegerQ[logfreqs],

Print[Length[xx]-1 is not a power of 2.″];

Abort[];

];

logSize=logFreqs+3;

samples = 2^logSize;

skip=samples/(2*freqs);

Print[″skip=″, skip];

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″];

undef=-2;

start=end=undef;

For[fi=0, fi<=freqs, ++fi,

i=fi+1;

If[Length[xx[[i]]]==0,

If[(start!=undef)&&(end==undef),

end=fi-1;

];

 , (*else*) 

If[start==undef,

start=fi;

];

];

];

Order=0; (*default is set later*)

Print[″start=″, start, ″, end=″, end];

For[fi=start, fi<=end, ++fi,

i=fi+1;

  freqIndex = fi*skip;

`` cycle=N[samples/freqIndex];

If[fi==start,

nrVars = Length[xx[[i]]];

Order=(nrVars-1)/2;

page7

Slavin.txt

odd=False;

, (*else*)

(*a previous result exists*)

Prev=Table[xx[[i-1, order+1+j]], {j, 1, order}];

  next = Tabl[xx[[i, order+1+j]], {j, 1, order};

{permute, modified, odd} = MatchMod[prev, next, cycle];

xx[[i]]=Flatten[{xx[[i,1]],

Table[xx[[[i,1+permute[[j]]]],{j,1,order}],

modified}];

];

If[xx[[i,1]]<0.0

xx[[i,1]]=-xx[[i,1]];

odd=! odd;

];

If[odd,

If[xx[[i, order+2]]>0.0,

xx[[i, order+2]]-=cycle＊0.5,

xx[[i, order+2]]+＝cycle＊0.5

];

];

(*Print[″end:xx[[″,i,″]]=″,xx[[i]]];*)

];

(＊initialize uncalibrated regions to allow later matrix transpose＊)

For[fi=0, fi<start, ++fi,

i=fi+1;

xx[[i]]=Table[0,{2*order+1}];

];

For[fi=end+1, fi<=freqs, ++fi,

i=fi+1;

xx[[i]]=Table[0,{2*order+1}];

];

qq = Transpose[xx];

K=qq[[1]];

Print[″K:″];

PlotFreq[K, start, end];

For[j=1, j<=order.++j,

Print["Amplitude", j];

PlotFreq[qq[[1+j]], start, end];

];

For[j=1, j<=order, ++j,

Print[″Delay″,j];

PlotFreq[qq[[order+1+j]], start, end];

];

(*gainA^(order-1)*gainB＝K, gainA applies to allfilters except one

which uses gainB＊)

totalGains = Take[K, {start+1, end+1}];

gainA＝Abs[totalGains]＾(1/order);

gainB＝totalGains/(gainA^(order-1));

delays=Table[Take[qq[[j+order+1]], {start+1, end+1}], {j, 1, order}];

filters=Table[Fj]terDesign[If[j==order, gainB, gainA], delays[[j]

    freqs, start, correctionFilterLen], {j, 1, order }];

For[j=1, j<=order, ++j,

ListPlot[filters[[j]], {Plotjoined->True, PlotRange->All}];

];

mags = Table[Magnitude[filters[[i]], samples], {i, 1, order}];

parts=Table[Take[mags[[j]],{1,1+freqs*skip,skip}],{j,1,order}];

Print[″Original Combined amplitude response:″];

PlotFreq[K, start, end];

Print[econstructed combined response:″];

Page8

Slavin.txt

ampl = Product[parts[[i]],{j,1,order}];

PlotFreq[amp], start, end];

Print[″Amplitude response difference:″];

PlotFre[K-amp],start,end];

refDlys=Table[GenDelay[fi]ters[[j]],2*freqs],{j,1,order}];

For[j=1, j<=order, ++j,

Print[″Reconstructed delay response″,j,″:″];

PlotRawDelay[refDlys[[j]]];

];

Return[filters];

]

(*single frequency DFT with gain normalization, returns a complex value*)

DFT[data_,m_]:=

Module[

{len, sum, t},

len=Length[data];

sum=0;

For[t=0, t<len, ++t,

   sum+＝(Exp[I＊2＊pi＊m＊t/len] data[[t+1]]); (＊1-based array indexing＊)

];

Return[2sum/len];

]

(*

＊Generates a broad-spectrum set of distorted calibration data. 

*'filters' = {filter2, filter3, ...}

*where filter2={{<taps>}, {<taps>}

＊wherefilter3＝{{<taps>},{<taps>}，{<taps>}} 

＊and each list<taps>has the same length′len′(must be odd). 

＊For example, before running Distortcal[], with

*'distortFilterLen=63;'and'spread=4.0;'

*filter2=GenFilters[distortFilterLen, 2, spread];

*filter3=GenFilters[distortFilterLen, 3, spread];

*filters={filter2, filter3};

＊′sx′is the output data record length for each test frequency pair. 

*'freqs'is one more than the number of output data frequency

*sets, e.g.for 127 sets, freqs=128. 

＊Note sx/freqs must be integer. 

＊′distortion is the distortion level in output lsb units

＊‘signal’is the two-tone amplitude in output lsb units. 

＊′noise′is the added noise standard deviation in output lsb units. 

＊Values that match the records acquired for the AD6645 are;

＊DataGen[filters, 1024, 128, 0.04, 15000.0, 2.0];

*)

DistortCal[fiters_, sx_, freqs_, distortion_, signal_, noise_]:=

Module[

{error, order, q, span, data, i, j, len, qfreq, f1, f2, phase1, phase2,

sines, samples, t, c, sum, maxOrder},

maxOrder = Lenqth[filters];

len=Length[filters[[1,1]]];

error=False;

If [! IntegerQ[sx],

Print["Analysis length is not integer."];

error=True;

];

If [! IntegerQ[freqs],

Print[″Frequency range is not integer.″];

error=True;

Page9

Slavin.txt

];

If [! IntegerQ[len],

Print[″Filter lengths are not integer.”];

error=True;

];

If[Mod[sx,4]! = 0,

print[″Analysis length is not a multiple of 4.″];

error=True;

];

If[Mod[sx, freqs]! = 0,

Print[″Analysis length is not a multiple of the number of″

"frequencies.n"];

error=True;

];

If [! OddQ[len],

Print["Filter lengths are not odd."];

error=True;

];

If[error,

Abort[];

];

c=(len+1)/2;

q=2.0＊Pi/sx;

span=sx/(2*freqs);

data = Table[0,{freqs-1}];

For[i=1, i<freqs, ++i,

qfreq=span*i;

f1=qfreq-1;

f2=qfreq+1;

phase1 = Random[Rea], {0, 2*pi}, 10];

Phase2 = Random[Rea], {0, 2*pi}, 10];

(*Length[sines]=sx+len-1*)

Sines=Table[(Cos[q*f1*t-phase1]+Cos[q*f2*t-phase2]),

{t,0,sx+len-2}];

(*Length[samples]=Length[sines]+1-len=sx*)

sum=0;

For[j=1, j<=maxOrder, ++j,

  sum + = GenDistortion[sines, filters[[j]]];

];

samples=distortion*sum;

samples+=Table[Gaussian[0, noise], {t, 1, sx}];

  samples+＝signal＊Take[sines, {c, sx+c-1}] 

(*output highest frequencies first*)

data[[Freqs-i]]=samples;

];

Return[data];

]

(＊＊)

DistortData[data_, filters_, distortion_, signal_, noise_]: =

Module[

{maxOrder, len, c, sum, j, nrSamples, samples, t},

maxOrder = Length[filters];

len=Length[filters[[1,1]]];

c = (len-1)/2;

sum=0;

For[j=1, j<=maxOrder, ++j,

  sum + = GenDistortion[data, filters[[j]]];

];

nrSamples = Length[sum];

Page10

Slavin.txt

Print[″nrSamples=″, nrSamples];

samples=distortion*sum;

If [noise! = 0,

samples+=Table[Gaussian[0, noise], {t, 1, nrSamples}];

];

samples+=signal* Take[data, {c+1, nrSamples+c}];

Return[samples];

]

(*

*'symbVars'is the list of symbolic variables;

* Flatten[{k, Table[a[j], {j, 1, order}], Table[d[j], {j, 1, order}}];

＊′ix′is the index symbol to use(passed in at a higher level). 

＊′sx′is the symbolic time scale factor.

*'dft'is a symbolic variable that is expanded to a list, used internally

＊ and also returned as the list′dfts′. 

*Returns;

＊'symbolicBins': a list of IM and harmonic frequencies as a function on of the

＊ frequency index'ix'.

＊ ′e2′: a symbolic error equation as a function of ix, sx, symbVars

* and'dfts'.

＊′g’: a list of symboli gradient equations that are derivatives of

＊ 'e2'with respect to each of the'symbVars'variable list. 

＊'dfts': a list of symbolic variables: {dft[1], dft[2],...} 

*)

DoInit[symbVars_, ix_, sx_, dft_]:=

Module[

{e2, g, j, nrAnalysisBins, symbolicBins, nrVars, dfts, dft, p},

(*

＊Create the symbolic expressions for the objective function e2

  ＊and all its partial derivatives with respect to 'vars'. 

＊Note: only really needs to be called once for multiple analysis

  ＊frequencies, as the results are purely symbolic. 

**)

{symbolicBins, p} = Terms[ix, sx, symbVars];

nrAnalysisBins = Length[symbolicBins];

dfts = Table[dft[j], {j, 1, nrAnalysisBins}];

(*e2 is the symbolic sum of the squares of the differences between the

  predicted vlues'p'and'dfts'at each analysis bin＊)

e2＝Sum[(p[[j]][[1]]-Re[dfts[[j]]])^2+

(p[[j]][[2]]-Im[dfts[[j]]])^2,{j,1,nrAnalysisBins}];

(*calculate symbolic partial derivatives*)

nrVars = Length[symbVars];

g = Table[D[e2, symbVars[[j]]], {j, 1, nrvars}];

Return[{symbolicBins, e2, g, dfts}];

]

(*

＊′e2′is symbolic error equation. 

＊′g′is symbolic gradiente equation. 

＊′dfts′is array of symbolic dft results. 

＊′symbVars′is list of veriables that BGGS is attempting to solve. 

＊′y′is list of dft results. 

＊′x′is set of start values for BFGS. 

＊′sx′is samples symbol ′s′is actual sample value. 

＊'ix'is frequency index symbol，'i'is actual index value. 

＊‘acc’ is BFGS accuracy goal.

＊'maxIterations'limits BFGS search time. 

Page11

Slavin.txt

*)

Dosolve[e2_, g_, dfts_, symbVars_, y0_, x_, sx_, s_, ix_, i_, acc_,

MaxIterations_]:＝

Module[

{e2n, gn, j, nrAnalysisBins, subst, substx, xnew, subst2, residual, len, y},

y = SetPrecision[y0, 39];

nrAnalysisBins = Length[y];

(*evaluate the symbolic expressions with the known constants*)

substx = Table[dfts[[j]]->y[[j]],{j,1,nrAnalysisBins}];

Subst = Flatten[{sx->s, jx->i, substx}];

e2n=e2/.subst;

gn=g/.subst;

xnew=BFGS[x,symbVars,gn,e2n,acc,maxIterations];

len = Length[xnew];

Subst2=Table[symbVars[[j]]->xnew[[j]],{j,1,len}];

Residul=e2n/.subst2;

Return[{xnew, residual}];

]

ExpFunc[q_]:＝If[q<-100, 0, Exp[q]]

(*takes sparse gain and delay lists, places them in a'size'length list

at "start'(zero based), for iterative filter calculations＊)

FilterDesign[gain_, dlyin_, size_, start_, requiredtaps_]:=

Module[

{l1, l2, ampl, delay, dly, n, taps, nrtaps, offset, i, j,

fftsize, deviation, avgdev, min, bestTaps},

l1=Length[gain];

l2=Length[dlyin];

If[l1! = l2,

Print["Missmatched gain and delay list lengths."];

Abort[];

];

If[EvenQ[requiredtaps],

Print[″Number of taps must be odd.″];

Abort[];

];

If[(start<=0)||((start+l1)>=size)

Print["start and size values illegal for", l1, "gains and delays."];

Abort[];

];

dly=dlyin+size; (*center the delays*)

n=size+1;

ampl = Flatten[{Table[0, {start}], gain, Table[0, {n-start-l1}]}];

delay=Flatten[{Table[size,{start}],dly,Table[size,{n-start-l2}]}];

If[Length[ampl]! = n,

Print[″Internalerror:Length[ampl]=″, Length[ampl],″,n=″,n];

Abort[];

];

If[Length[delay]! = n,

Print[″Internal error:Length[delay]=″, Length[delay],″,n=″,n];

Abort[];

];

fftsize=2*size;

offset=Floor[reqiredtaps/2];

min=Infinity;

For[j=0, j<100, ++j, (*note:Break[] for small j in loop*)

Taps = Taps[ampl, delay];

nrtaps = Length[taps];

(*set coefficients outside desired filter size to zero

Page12

Slavin.txt

For[i=0, i<fftsize, ++i,

If[(i<(size-offset))||(i>=(size+offset)),

Taps[[i+1]]=0;

];

];

{ampl, delay} = Response[taps];

(* overwrite gain and delay with original values in critical

region*)

deviation=0.0;

For[i=start, i<start+l1, ++i,

deviation+＝Abs[delay[[i+1]]-dly[[i-start+1]]];

];

Avgdev＝deviation/(l1+1);

If[avgdev<min,

min=avgdev;

bestTaps = taps;

];

Print[″Loop″,j,″:average absolute delay deviation=″,avgdev];

If[j＝＝25,

Break[];

];

For[i=start, i<start+l1, ++i,

Ampl[[i+1]]＝gain[[i-start+1]];

Delay[[i+1]]=dly[[i-start+1]];

];

    For[i=0, i<start, ++i,

Delay[[i+1]]=delay[[start+1]];

];

];

Taps = Table[bestTaps[[i+1]], {i, size-offset, size+offset}];

dcdelay = DCdelay[taps];

margin=Floor[0.45 requiredtaps];

If[dcdelay<margin,

Print[″limiting lower,″];

Taps = DCadjust[taps, margin];

{ampl, delay} = Response[taps],

(*else*)

If[dcdelay>=requiredtaps-margin,

     Print[″limiting upper,″];

Taps = DCadjust[taps, required taps-margin-1];

{ampl, delay} = Respons[taps];

];

];

Return[taps];

]

(*x and ds are 2-D column vectors, 'vars' is a list of variables in the

Gradients(derivatives of the objective function). Return the objective

function value at vars=x*)

FletcherLineSearch[x_, ds_, vars_, gradients_, objectiveFunc_]:=

Module[

{m, tau, rho, sigma, gamma, mhat, epsilon, xk, s, f0, gk, deltaf0, dk, aL, aU, fL, dfL,

a0, a0hat, a0Lhat, a0Uhat, gtemp, df0, deltaa0, alpha, div},

m=0;

sigma=0.1;

tau=Min[sigma, 0.1]: (*tau<=sigma and tau<=0.1*)

rho=sigma; (*rho<=slgma*)

gamma=0.75;

mhat=400;

Epsilon=1.0*10^-10;

xk=x; (*xk is a column vector*)

s=ds; (*'s'a column vector*)

Page13

Slavin.txt

(*step 1: evaluate fk and gk at k＝0.f0 is a scalar，g0 is a column

vector*)

f0 = ObjectiveFunc[xk, vars, objectiveFunc];

gk = GradFunc[xk, vars, gradients]; (*gk is a column vector*)

m+=2;

deltaf0=f0;

(*step2:initialize the search*)

dk=s; (*dk is a column vector*)

aL=0.0;

aU＝1.0＊＊10^99;

fL=f0;

dfL＝(Trranspose[gk].dk)[[1]][[1]];

If[Abs[dfL]>epsilon,

If[div==0.0,

Print[″FletcherLineSearch divide by dfL＝0.0″];

Abort[];

];

a0＝-2.0＊deltaf0/dfL;

,

a0=1.0;

];

If[(a0×10^-9)||(a0>1.0),

a0=1.0;

];

(*step3: ＊)

While[True,

deltak＝a0＊dk; (＊deltak is a column vector＊)

f0 = ObjectiveFunc[xk+deltak, vars, objectiveFunc];

m+=1;

If[(f0>(fL+rho＊(a0-aL)＊dfL))&&

(Abs[fL-f0]>epsilon)&&(m<mhat),

If[a0<aU,

aU=a0;

];

div=2.0*(fL-f0+(a0-aL)*dfL);

If[div==0.0,

Print[″FletcherLineSearch divide by div＝0.0″];

Abort[];

];

a0hat＝aL+(dfL＊(a0-aL)^2)/div;

a0Lhat＝aL+tau＊(aU-aL);

If[a0hat<a0Lhat,

a0hat＝a0Lhat;

];

a0Uhat＝aU-tau＊(aU-aL);

If[a0hat>a0Uhat,

a0hat=a0Uhat;

];

a0=a0hat;

(*else*)

(*gtemp is a vector*)

gtemp = GradFunc[xk+a0*dk,vars,gradents];

df0 = (Transpose[gtemp].dk)[[1]][[1]];

m+=1;

If[(df0<sigma＊dfL)&&(Abs[fL-f0]>epsilon)&&

(m<mhat)&&(dfL!=df0),

div = dfL-df0;

  deltaa0=(a0-aL)*df0/div;

If[deltaa0<=0,

A0hat=2*a0;

, (*else*)

A0hat=a0+deltaa0;

];

Page14

Slavin.txt

a0Uhat＝a0+gamma＊(aU-a0);

If[a0hat>a0Uhat,

A0hat＝a0Uhat;

];

aL=a0;

A0=a0hat;

fL=f0;

dfL=df0;

, (*else*)

Break[];

];

];

];

If[a0<1.0＊10^-5,

  alpha=1.0＊10^-5; 

, (*else*)

alpha=a0;

];

Return[alpha];

]

FTSize[len_,min_]:=

Module[

{},

Return[Max[min, 2^Ceiling[Log[2,len]]]];

]

Gaussian[mean_, sd_] :=

Module[

{result},

Return[sd＊Sqrt[2.0]*InverseErf[Random[]]＊Sign[Random[]-0.5]+mean];

]

(＊finds sample delay as a function of frequency. 

uses ALL'list'taps-no symmetry is assumed＊)

GenDelay[list_,min_:16]:=

Module[

{len, delays, endindex, nx},

len=Length[list];

`nx=FTSize[len,min];

endindex=nx/2;

Return[Take[Delays[list,nx],{1,endindex+1}]];

]

(*

＊distort the data for a particular filter order Length[taps]. 

*)

GenDistortion[data_, taps_] :=

Module[

{len, filters, flen, delay, output, i, file},

len=Length[data];

filters = Length[taps];

flen = Length[taps[[1]]];

If[EvenQ[flen],

Print["Filter lengths are not all odd."];

Abort[];

];

For[i=1, i<=taps1, ++i,

len1=Length[taps[[i]]];

Page15

Slavin.txt

If[len1! =flen,

Print["Filter lengths are not all equal."];

Abort[];

];

];

delay=(1+flen)/2;

output = Table[1.0, {len-flen+1}];

For[i=1, i<=filters, ++i,

filt = taps[[i]];

If[Length[filt]! =flen,

Print["Filter lengths are not all the same."];

Abort[];

];

output*=ListConvolve[filt, data];

];

Return[output];

]

(＊generate a set of filters for a particular order of distortion＊) 

GenFilters[len_, order_, spread_]:=

Module[

{error, taps, i, j, c},

error=False;

If [! IntegerQ[len],

Print[″Filter lengths are not integer.″];

error=True;

];

If [! OddQ[len],

Print["Filter lengths are not odd."];

error=True;

];

If[error,

Abort[];

];

c = (len+1)/2;

Taps=Table[Gaussian[0,1]*ExpFunc[-Abs[(j-c)*2/spread]],

{i, 1, order}, {j, 1, len}];

Return[taps];

]

(*

*'samples' gives nr of output sampels(plus extra), and also fundamental

＊length, so all generated sines are harmonics of the fundamental. 

＊′extra′gives the number of extra samples in the output record that may be

＊needed to compensate for convolution data loss due to filters. 

＊′freqLo′is the low normalized frequency for sine output: 

＊0<freqLow<0.5.The nearest higher harmonic of the'samples'frequency

＊is the lowest output frequency. 

＊′freqHi′is the high normalized frequency for sine output:

＊0<freqLow<0.5.The nearest lower harmonic of the'samples'frequency

＊is the highest output frequency. 

*'freqInterval'is the normalized frequency step, also a harmonic of the

＊‘samples’frequency. 

*)

GenTest[samples_, extra_, freqLo_, freqHi_, freqInterval_]: =

Module[

{kInterval, kStart, kEnd, sum, x, j, t},

kInterval=Floor[freqInterval*samples+0.5];

Page16

Slavin.txt

kStart = Ceiling[freqLo*samples];

kEnd = Floor[freqHi*samples];

sum=0;

x=2.0＊Pi/samples;

For[j=kStart, j<=kEnd, j+=kInterval,

   sum+=Table[sin[x*j*t], {t, 0, samples+extra-1}];

];

Return[sum];

]

(*′x′is a 2-D column vector,′vars′is a list of variables in the list of

′gradients′functions which are the partial derivatives of the objective

function with respect to each variable.＊)

GradFunc[x_, vars_, gradients_] :=

Module[

{g, subst, j, len},

len=Length[x];

Subst = Table[vars[[j]]->x[[j]][[1]],{j,1,len}];

g = gradients/.subst;

Return[Transpose[{g}]]; (*return a column vector as a 2-D matrix*)

]

IsSame[xx_]:＝

Module[

{nrVars, order, indices, permArray, perms, same, i, j, perm, v0, v1}

nrVars = Length[xx];

order=(nrVars-1)/2;

indices = Table[i,{i,1,order}];

permArray = Permutations[indices];

perms = Length[permArray];

same=False;

For[j=1, j<=perms, ++j,

Perm = permArray[[j]];

v0 = perm[[1]];

v1 = pem[[2]];

If[v1×v0,

Same=(Chop[xx[[v0+order+1]]-xx[[v1+order+1]]]==0);

If[same,

Break[];

];

];

];

Return[same];

]

Log10List[list_]:＝

Module[

{minv, lminv, log, i, item, L, result},

L = Length[list];

minv=Min[list];

If[minv<=10^-10,

   minv = 10^-10; 

];

  result = TaBle[0,{L}]; 

For[i=0, i<L, ++i,

     item = list[[i+1]]; 

If[item<minv,

`item=minv;`

];

   result[[i+1]]＝20.0 Log[10, item] 

Page17

Slavin.txt

];

Return[result];

]

LogFFT[list_]:＝

[Mdule[

{i, L, magfreg, logmagf, data, bins},

L = Length[list];

bins=Floor[(L+1)/2];

magfreq = Abs[Take[Fourier[list, {FourierParameters->{-1, 1}}], bins]];

  freq=Table[N[i/L], {i, 0, bins-1}];

` logmagf = Log10List[magfreq];

data = Transpose[{freq, logmagf}];

Return[data];

]

Magnitude[list_,n_]:=

Module[

{i, L, reals, magfreq},

L = Length[list];

Reals=Join[list, Table[0, {i, 1, n-L}]];

magfreq = Abs[Take[Fourier[reals, {FourierParameters->{1, 1}}], n/2+1]];

Return[magfreq];

]

(*

＊patt1 and patt2 numeric list lengths are equal.Find the permutation index

＊of patt2 that minimizes the sum of the squares of the differences of patt1

*and patt2, modulo: 'mod'.-Used for reorganizing phases. Also returns a

＊modified version of patt2 which is close to patt1, modulo 'mod'. the

＊For example:

＊MatchMod[{0.2, 0.5, 0.8}, {1.22, 0.88, 0.47}, 1] 

＊yields the indicesmodified patt2, and odd status:

* {{1, 3, 2}, {0.22, 0.47, 0.88}, False} 

*)

MatchMod[patt1_,patt2_,mod_]:=

Module[

{len1, len2, indices, permArra, perms, min, sum, perm, i, j, diff, index, x, new,

 step, steps, total, odd},

len1=Length[patt1];

len2=Length[patt2];

If[len1! =len2,

Print["Length mismatch."];

Abort[];

];

indices=Table[i, {i, 1, len1}]

permArray = Permutations[indices];

perms = Length[permArray];

min=Infinity;

step=mod*0.5;

For[j=1, j<=perms, ++j,

sum=0;

Perm = permArray[[j]];

total=0;

    For[i=1, i<=len1, ++i,

diff = patt1[[i]]-patt2[[perm[[i]]]];

  steps = Floor[diff/step+0.5];

total+=steps;

     sum+=(diff-step*steps)^2; 

Page18

Slavin.txt

];

If[sum<min,

(*found best solution so far*)

   odd＝OddQ[total; (＊record the odd-ness of the solution＊) 

index=j;

min=sum;

];

];

perm = permArray[[index]];

new=Table[x=patt2[[perm[[i]]]];

x+step*Floor[(patt1[[i]]-x)/step+0.5],

{i,1,len1}];

Return[{perm, new, odd}];

]

(*′x′is a 2-D column vector,′vars′is a list of all variables in the

Objective function.Return the objective function value at vars＝x＊)

ObjectiveFunc[x_, vars_, objectiveFunc_] :=

Module[

{subst, j, len},

len=Length[x];

Subst = Table[vars[[j]]->x[[j]][[1]],{j,1,len}];

Return[objectiveFunc/:subst]; (*return a scalar*)

]

PlotColor[list_, color_]:=

Module[

{min, max, r, g, b, data},

r = color[[1]];

g = color[[2]];

b = color[[3]];

data={0,Transpose[list][[2]]};*always include zero*)

max = Max[data];

min=Min[data];

Return[ListPlot[list, {PlotJoined->True,

GridLines->Automatic,

PlotRange->{min,max},

ImageSize->600,

PlotStyle->RGBColor[r,g,b]}]];

]

(*startIndex and endIndex indices(0-based)apply to the list*)

plotFreq[list_, startIndex_:0, endIndex_:-1]:=

Module[

{i, L, delays, freq, data, idata, endi, max, min, selected},

L = Length[list];

If[endIndex==-1,

endi=L-1,

endi＝endIndex

];

freq=Table[N[i/(2*(L-1))], {i, startIndex, endi}];

selected = Take[list, {startIndex+1, endi+1}];

data = Transpose[{freq, selected}];

idata={0, selected}; (*always include zero*)

max = Max[idata];

min=Min[idata];

Return[ListPlot[data,{PlotJoined->True,

GridLines->Automatic,

PlotRange->{{0,0.5},{min,max}},

ImageSize->600,

Page19

Slavin.txt

PlotStyle->RGBColor[1, 0, 0]}]];

]

(＊finds sample delay as a function of frequency. 

uses ALL 'list' taps -no symmetry is assumed，plot in red＊)

PlotRawDelay[delays_]:＝

Module[

{i, len, freq, data},

len=Length[delays]-1;

freq=Table[N[i/(2*len)], {i, 0, len}];

data = Transpose[{freq, delays}];

Return[PlotColor[data,{1,0,0}]];

]

PlotUnwindowedFFT[list_]:＝

Module[

{},

Return[WfmColor[LogFFT[list],{1,0,0}]];

]

PolarToComplex[ampls_, angles_]:=

Module[

{l1, l2, fl, i, len, reals, angle, ampl, real, imag, out},

l1=Length[ampls];

l2=Length[angles];

If[l1! = l2,

Print["Missmatched amplitude and angle list lengths."];

Abort[];

];

len=l1-1;

`fl=2*len;

out = Table[0,{i,1,fl}];

For[i=0, i×len, ++i,

angle = angles[[i+1]];

ampl = ampls[[i+1]];

real = ampl Cos[angle];

imag = ampl Sin[angle];

  out[[i+1]]=real+I imag;

If[(i != 0),

    out[[fl-i+1]]＝real-I imag; 

];

];

out[[len+1]]=ampls[[len+1]];

Return[out];

]

PrAmplPhase[data_, bin_]:=

Module[

{out,ampl,ph},

out = DFT[data, bin];

ampl = Abs[out];

ph = Arg[out];

Return[{ampl,ph}];

]

ProcessCalData[data_, order_, lowf_, hiF_, acc_]: =

Module[

{fi,

Page20

Slavin.txt

delta, freqs, start, end, skip, i, j, t, freqIndex, overhang,

i1, i2, ampl, ph, len, bins, sBins, w1, w2, scaling, temp, cycle,

amplRet, phRef, amplRatio, q, scale,

amplitude2, relPhase2, relDelay2, amplitude3, relPhase3, relDelay3,

corrPh, turns, rph, xbin, analyze1, reference, redo, d1, d2, x, y,

dly, dlyA, dlyB, amplA, amplB, donesomething,

descr, expected, K, logSize, logFregs, amplx,

magA, magB, partAm1, partAp1, partBm1, partBp1, transpose,

delays, base, partial, de, amplI, ii, nrBins,

filter1, filter2, frBins, index, f1, f2, fs, out1, out2, out3,

e2, g, s, nrVars, symbolicBins, symbVars,

k, a, d, ix, sx, xnew, yl, allSymBins, blen, slen, numericBins,

dfts, dft, results, maxorde, gainA, gainB, delavA, delayB,

maxIterations, substx, jj, kk, random,

residual, lowestResiduals, bestResults, startPoints, xStart, gains,

solutions, found, xx, qresults, solution, smallesttResidueIndex,

sortList, sorted, state, const, maxTries, extraTries,

symbvarsW, symbolicBinsW, e2W, gW, dftsW},

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊＊＊＊＊＊"];

Print

Print［″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊＊＊＊＊＊"];

(＊Various″magic″parameters＊) 

(＊maximum tries to find a first BFGS solution. Too large means time

may be wasted finding a solution when none can be found. Too small

and a solution may not be found.＊)

maxTries=100;

(＊extra attempts to find more BFGS solutions after the first is found. 

Too large and time may be wasted finding additional solutions when

None exist. Too small, and a better solution can be missed.The

Probability of missing a second solution is about 1/2 per try,

So the probability of missing it is(1/2)^extraTries＊)

extraTries = 25;

(＊maximum BFGS iterations allowed in which to converge to a solution. 

Too large and time may be wasted trying to find a BFGS solution,

too small and a solution may not be found.Typically, 2nd order

needs<30 and 3rd order needs<80.＊)

maxIterations = 100;

(＊maximum number of total unique solutions at a cal frequency. 

2nd order: 3^2＝9 solutions per group, 3rd order: 3^3＝27 per group. 

Too larg a value and more memory is used. Too small and there is

Insufficient storage for multiple groups. I′ve never seen more than 3

groups, or 3*27=81 solutions for any frequency.*)

solutions=200;

(＊sampling frequency(MHz)＊)

fs=100.0;

overhang=63;

logFreqs=7;

symbVars =

Flatten[{k, Table[a[j], {j, 1, order}], Table[d[j], {j, 1, order}]}];

{symbolicBins, e2, g, dfts} = DoInit[stmbvars, ix, sx, dft];

If[(lowF<0)||(lowF>0.5)||(hiF<0)||(hiF>0.5)||

(lowF>=hiF),

Print[″low or high frequency bounds are<0 or×0.5, or″];

Print["low frequency bound", lowF, "is>=hidh bound", hiF];

Abort[];

];

Print["sampling frequency:", fs, "MHz"];

logSize＝logFreqs+3; (＊added value must be>＝3＊) 

Page21

Slavin.txt

delta=1;

samples=2^logSize; (*power of 2*)

n＝(2＊overhang+1); (＊n is odd filter length＊) 

len=samples+2*(n-1); (*overhang required at each end*)

If[logFreqs>logSize-3,

Print[″logFreqs must be<=logSize-3″];

Abort[];

];

freqs＝2^logFregs; (＊frequencies to Nyguist. A power of 2＊)

start=Max[Floor[lowF*freqs*2], 1];

end=Min[Ceiling[hiF*freqs*2], freqs-1];

If[start==freqs/2,

Print[″Error: start frequency evaluates to fs/2.This cannot be″

"used because it starts on an alias notch."];

Abort[];

];

If[end==freqs/2,

Print[″Error: end frequency evaluates to fs/2.This cannot be″

"used because it starts on an alias notch."];

Abort[];

];

skip=samples/(2*freqs);

Print["Frequency indices: start = ", start, ", end = ", end];

Print["Range:0:", freqs-1, ", index skip=", skip];

(＊nr of different IM and harmonic frequency bins to consider＊) 

nrBins = Length[symbolicBins];

amplitude=Table[0, {i, 1, freqs+1}, {j, 1, nrBins}];

yl=Table[0,{i,1,freqs+1},{j,1,nrBins}];

bins=Table[0,{i,1,freqs+qs+1},{j,1,nrBins}]; 

amplRatio=Table[0.0,{i,1,freqs+1}];

Scaling=Table[0.0,{i,1,freqs+1}];

descr=Table[0.0,{i,1,freqs+1}]; 

allSymBins = AllTerms[ii, orde];

blen = Length[allSymBins];

dly=0.0;

redo=-1;

const=2*Pi/samples;

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊＊＊＊＊＊"];

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ DFT analysis＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊＊＊＊＊"];

(*

*default x[[i]]=False indicates not measurable due to aliasing,

＊True indicates measurable but failed to converge on search,

＊otherwise x[[i]]equals a list of solutions. 

*)

x = Table[False, {i, 1, freqs+1}];

For[fi=start, fi<=end, ++fi,

i=fi+1; (*offset for non-zero based array indexing*)

index = freqs-fi;

analyze1 = data[[index]];

freqIndex = fi skip;

i1 = freqIndex-delta;

i2=freqIndex+delta;

w1=2*Pi*i1/samples;

w2=2*Pi*i2/samples;

Page 22

Slavin.txt

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″];

Print["fi=", fi, ", index=", index, ", freqIndex=", freqIndex];

(＊look for alias collisions using all frequencies＊)

numericBins = allSymBins/.ii->freqIndex;

sBins = Sort[AliasF[numericBins, samples]];

Print[″sBins=", sBins];

For[j=1, j<blen, ++j,

If[sBins[[j]]==sBins[[j+1]],

Redo = i;

Break[];

];

];

If[redo==i,

Continue[]

];

x[[i]]=True; (*indicate measurable at i*)

(*symbolicBins represents only the frequencies used for analysis,

which excludes Dc and the ibput cal frequencies themselves＊)

bins[[i]]=symbolicBins/.ix->freqIndex;

Print[″bins″, order, ″=″, bins[[i]]];

(*generate timing references from signal+small distortion*)

reference=analyze1^order;

temp = Abs[DFT[analyze1, i1]];

temp+=Abs[DFT[analyze1,i2]];

scaling[[i]]=temp*0.5;

(*calculate amplitudes, phases and delays*)

For[j=1, j<=nrBins, ++j,

xbin = bins[[i, j]];

temp = DFT[analyze1, xbin];

{amplRef, phRef = PrAmplPhase[reference, xbin];

yl[[i,j]]=temp*Exp[-I*phRef];

];

f1=N[fs*i1/samples, 8];

f2=N[fs*i2/samles, 8];

Print[″f1=", fs-f1,"(aliased=", f1,"), ",

″f1=″, fs-f2, ″(aliased=″, f2,″).″];

Print[″yl[″,i,″]=″,yl[[i]]];

If[redo==(i-1),

redo=-1;

];

];

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊"];

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊"];

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊"];

search=1; (*related search radius*)

For[fi=start, fi<=end, ++fi,

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″];

i=fi+1; (*offset for non-zero based array indexing*)

If [! x[[i]];

Continue[];

];

Page23

Slavin.txt

cycle=samples/freqIndex;

(＊measurable＊) 

freqIndex = fi*skip;

Print["fi=", fi, ", freqIndex=", freqIndex,

″, cycle delay(samples) = ″, N[cycle]];

Print[″y[″, fi, ″]=″, yl[[i]]];

lowestResiduals = Table[Intinity, {solutions}];

bestResults = Table[{}, {solutions}];

startPoints = Table[{},{solutions}];

solution=0;

smallestResidueIndex = 1;

kk=maxTries;

nestStart = {};

random=False;

While[--kk>=0,

slen = Length[x[[i-1]]];

random=(slen==0)||(kk!=maxTries-1)||IsSame[x[[i-1]]];

If[random,

If[slen==0,

xStart＝Flatten[

{0.1*(Random[]-0.5),

       Table[Exp[0.1*(Random[]-0.5)],{order}], 

Table[cycle*(Random[]-0.5)*0.5,{order}]}];

, ,

xStart=x[[i-1]]

Table[Exp[0.4*(Random[]-0.5)],{slen}];

];

, ,

xStart = x[[i-1]];

];

{results, residual} = DoSolve[e2, g, dfts, symbVars,

yl[[i]], xStart, sx,

samples, ix, freqIndex, acc,

maxIterations];

If[Length[results]<=0,

Continue[];

];

If [! random,

    Print[″result=″,results]; 

Print[″residual=″, residual];

x[[i]]=results;

Break[];

];

(＊center the angles as much as possible＊)

xStart = results;

parity=0;

For[jj=1, jj<=order, ++]],

xx=Floor[[results[[order+jj+1]]*(2.0/cycle))+1/2];

xStart[[order+jj+1]]-＝(xx＊cycle＊0.5);

Parity+＝xx;

];

If[OddQ[parity],

xStart[[1]]=-xStart[[1]];

];

(＊reoptimize from more centered starting point＊)

{results, residual} = DoSolve[e2, g, dfts, symbVars,

yl[[i]], xStart, sx,

samples, ix, freqIndex, acc,

maxIterations];

If[Length[results]<=0,

Continue[];

Page24

Slavin.txt

];

sortList＝{Table[results[[jj+order+1]], {jj, 1, order}],

      Table[results[[jj+1]],{jj,1,order}], 

     Table[xStart[[jj+order+1]],{jj,1,order}], 

      Table[xStart[[jj+1]], {jj, 1, order}]}:

(＊rank solutions in order of second sub-list (angles)＊)

sorted=Transpose[Sort[Transpose[sortList]]];

results=Flatten[{results[[1]], sorted[[2]], sorted[[1]]}];

xStart = Flatten[{xStart[[1]], sorted[[4]], sorted[[3]]}];

(＊see if solution found already＊)

found=False;

For[jj=1, jj<=solutions, ++jj,

If[lowestResiduals[[jj]]==Infinity,

(＊reached end of valid solutions，nothing found＊)

Break[];

, (*else*)

If[Chop[lowestResiduals[[jj]]-residual]＝＝0.0，

found=True;

Break[]; (*solution already found*)

];

];

];

If[found,

Continue[];

];

(＊no solution found＊)

If[jj>solutions,

(*finding another solution when no more are allowed*)

Print["＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊",

"Too many solutions found.",

″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″];

Abort[];

];

(＊new solution, so reset loop limit＊)

kk = extraTries;

(＊find index into data set with lowest error residual＊)

If[jj>1,

If[residual<lowestResiduals[[smallestResidueIndex]],

The smallestResidueIndex = jj;

];

];

(*find a multiplicity of local solutions at approximately Pi

intervals on all angles＊)

lowestResiduals[[jj]]=residual;

bestResults[[jj]]=results;

startPoints[[jj]]=xStart;

++solution;

Print["Solution", solution];

Print[″Starting vector x=″, xStart];

Print["Minimum at", results];

Print[″Lowest residual error=″, residual];

++jj;

state=Table[-search,{order}];

While[True,

(*find next state*)

xStart = results;

(*Print[″state=″, state];*)

parity=0;

For[xx=1, xx<=order, ++xx,

Page25

Slavin.txt

Parity+=state[[xx]];

  xStart[[order+1+xx]]+＝(cycle＊state[[xx]]＊0.5);

];

If[OddQ[parity],

xStart[[1]]=-results[[1]];

];

{qresults, residual} = DoSolve[e2, g, dfts, symbVars,

y[[i]], xStart, sx,

samples, ix, freqIndex, acc,

maxIterations];

If[Length[qresults]>0,

If[jj>solutions,

                  

       "Too many solutions found.", 

  ″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″];

Abort[];

];

lowestResiduals[[jj]]=residual;

BestResults[[jj]]=qresults;

startPoints[[jj]]=xStart;

++jj;

];

xx=1;

While[xx<=order,

State[[xx]]+=1;

If[state[[xx]]<=search,

Break[];

];

State[[xx]]=-search;

      ++xx;

];

If[xx>order,

Break[];

];

]; (*end of While[True,]*)

]; (*end of While[--k>=0,]*)

If[random,

If[Length[bestResults[[1]]]==0,

Print["Solutions not found."];

(*leave x[[i]] as True indicates measurable but not valid*)

x[[i]]=bestResults[[smallestResidueIndex]];

];

];

];

(*

＊for each i, x[[i]] is either a legitimate list of solutions, or is

＊boolean False indicating it is not measurable due to alias frequency

＊conflicts, or is True indicating it is measurable, but that a

＊solution hasn′t been previously toound. 

*)

Print[″Filling in holes using search:″];

(＊fill in holes where no solution was found＊)

doneSomething = True;

While[doneSomething,

doneSomething = FALSE;

for[fi=start, fi<=end, ++fi,

i＝fi+1; (*offset for non-zero based array indexing*)

If[(Length[x[[i]]]==0)&& x[[i]]],

(*x[[i]] is measurable*)

    freqIndex = fi*skip; 

Page26

Slayin.txt

If[(fi>start)&&(Length[x[[i-1]]]>0),

xStart=x[[i-1]];

{results, residual} = DoSolve[e2, g, dfts, symbVars,

yl[[i]],xStart,sx,

samples, ix, freqIndex, acc,

maxIterations];

If[Length[results]>0,

x[[i]]=results;

doneSomething=True;

Print[″x[″,fi,″]=″,x[[i]]];

Print[″residual=″, residual];

];

];

If[(Length[x[[i]]]==0)&&

... (fi<end)&&(Length[x[[i+1]]]>0),

xStart=x[[i+1]];

{results, residual} = DoSolve[e2, g, dtts, symbVars,

yl[[i]],xStar,sx,

samples, ix, freqIndex,

acc, maxIterations];

If[Length[results]>0,

x[[i]]=results;

doneSomething=True;

Print[″x[″,Fi,″]=″x[[i]]];

Print[″residual=″,residual];

];

];

];

];

];

Print[″Forcibly filling all renaining holes using previous valuesls:″];

For[fi=start+1, fi<=end, ++fi,

     i＝fi+1; (*offset for non-zero based array indexing*) 

If[Lenth[x[[i]]]==0,

x[[i]]=x[[i-1]];

];

];

Print[″Scale output values:″];

(*scale as separate step to avoid problems with valid handling*)

For[fi=start, fi<=end, ++fi,

i=fi+1;

If[scaling[[i]]==0.0,

Temp = (scaling[[i-1]]+scaling[[i+1]]*0.5;

           

temp = scaling[[i]];

];

x[[i,1]]/=(temp^order);

(*Print[″output:x[″, fi,″]=", x[[i]]];*)

];

Return[x];

]

(*finds amplitudes and delay (in sample units) response of even length real

tap filter*)

Response[real]taps_]:＝

Module[

{taps, dcDelay, complex, len, l1, ampls, delays, i, ampl,

radiansPerSample, delay, prevangle, angle, diffangle, newdelay, mod,

scale, xtaps},

taps = Length[realtaps];

Page27

Slavin.txt

(*Print[″realtaps=″,realtaps];*)

(*calculate the dc delay without Arg[] aliasing.*)

dcDelay = DCdelay[realtaps];

(＊extend input to scale＊taps to use a longer FFT which

should reduce phase wrapping ambiguities＊)

scale=2;

xtaps = Flatten[{realtaps, Table[0, {(scale-1)taps}]}];

(*taps input returns l1=(taps/2)+1 items*)

complex = Fourier[xtaps, {FourierParameters->{1, 1}}];

(*Print[″complex=″, Take[complex, {1, taps*scale, scale}]];*)

len=scale*taps/2;

l1=len+1;

ampls = Table[0,{l1}];

delays = Table[0,{l1}];

delays[[1]] = dcDelay;

ampls[[1]]=Abs[complex[[1]]];

mod = Pi;

(*for all non-dc frequencies*)

For[i=1, i<=len, ++i,

ampl = Abs[complex[[i+1]]];

ampls[[i+1]]=ampl;

If[ampl<10^-12,

(*if amplitude too small for accurate angle, use previous

  delay as this is needed by the delay reconstruction＊) 

Delays[[i+1]]=delays[[i]],

(* else*)

(*reconstruct delay using previous delay which should

  work well as long as delay changes, when translated

     to angles, are not>Pi bwtween adjacent frequencies

at the new frequency＊)

  (＊conversion ratio between delays and radians＊) 

radiansPerSample=N[Pi*i/len];

delay = delays[[i]]; (*get previous delay*)

(＊convert previous delay to angle at new frequency＊)

Prevangle=delay＊radiansPerSample;

(*get angle at new frequency*)

 anale=Mod[Arg[complex[[i+1]]], Pi, prevangle-Pi/2];

(*Print[″angle=″, angle]:*)

NewDelay = angle/radiansPerSample;

Delays[[i+1]]=newDelay;

];

];

(*where magnitude is too small, angles are unreliable, so

Average adjacent delays to get delay instead＊)

For[i=1, i<len, ++i,

If[ampls[[i+1]]<10^-12,

Delays[[i+1]]=(delays[[i]]+delays[[i+2]])/2;

];

];

If[ampls[[len+1]]<10^-12,

delays[[len+1]]=2 delays[[len]]-delays[[len-1]];

];

Return[{Take[ampls, {1, len+1, scale}],

Take[delays, {1, len+1, scale}]}];

]

Page28

Slavin.txt

(*in each list, first value is DC, last is at w=Pi=fs/2*)

Taps[ampls_, delays_]:=

Module[

{l1, l2, len, taps, angles, complex, targetDelay, realtaps},

l1=Length[ampls];

l2=Length[delays];

If[l1! = l2,

Print[″Missmatched amplitude and delay listl lengths.″];

Abort[];

];

len=l1-1;

Taps=2*len;

angles=Table[delays[[i+1]]*i*Pi/len, {i, 0, len}];

(*Print[″angles=″, angles];*)

complex = PolaroComplex[ampls, angles];

(*Print[″complex=″, complex];*)

realtaps = InverseFourier[complex, {FourierParameters->{1, 1}}];

Return[realtaps];

]

(*

＊Finds symbolic IM and harmonic frequency terms using the frequency

＊index and samples symbols 'ix'and'sx'.'vars'is the list of variables

＊defining the filter product to use.It is of the form;

* Flatten[{k, Table[a[j], {j, 1, order}], Table[d[j], {j, 1, order}]}];

*)

Terms[ix_,sx_,vars_]:=

Module[

{r0, r1, r2, r3, r4, r5, t, tt, a, b, c, j, cost, sint, len, i1, i2, p, q, nrVars, order,

q0, q1, q2, q3, q3a, q4, q5, q6, q7, q8, terms, j1, j2, t1, t2},

nrVars = Length[vars];

order=(nrVars-1)/2;

If [! IntegerQ[order],

Print[″Vars list must be of from:″

  "{kx,ax0,ax1,...,axn,dx0,dx1,...,dxn}"];

Abort[];

];

q0=TrigReduce[(Cos[(ix-1)tt]+Cos[(ix+1)tt])^order];

(＊extract tt variable outside as a common factor for each cos term＊)

q1=q0/.Cos[a_]:>Cos[Apart[a]];

(＊expand into separate terms＊)

q2 = Expand[q1];

(*make into a list*)

q3 = Table[q2[[j]], {j, 1, Length[q2]}];

(*eliminate DC frequencies from list*)

q4=Table[If[NumericQ[q3[[j]]], 0, q3[[j]]/(q3[[j]]/.tt->0)],

        {j, 1, Length[q3]} 

(*eliminate the Cos[] around each term*)

q5=q4/.Cos[a_]:>Apart[a];

(*form(k*ix+m)*tt->k*ix+m for integer k, m*)

q6=q5/.a-*tt:>a;

(＊make ix term positive:if k is-ve，then its derivative is-ve＊)

q7=Table[If[D[q6[[j]],ix]<0,-q6[[j]],q6[[j]]],{j,1,Length[q6]}];

(*eliminate lower-order harmonics from term list*)

If[order>=2,

j1=j2=1;

  (＊a higher order contains all the terms in the lower order＊) 

t1=Sort[Expand[q7]];

t2 = Sort[AllTerms[ix, order-2]];

len1=Length[t1];

Page29

Slavin.txt

len2=Length[t2];

While[j2<=len1,

If[t1[[j1]]==t2[[j2]],

t1[[j1]]={};

If[j2==len2,

Break[];

];

      j2+=1; 

];

     j1+=1; 

];

terms = Flatten[t1];

, ,

terms = q7;

];

len=Length[terms];

(＊convert into form:Cos[(k＊ix+m)＊2＊Pi＊t/sx]＊)

Cost=terms/.a_:>Cos[Apart[a＊2＊Pi＊t/sx]];

(*find corresponding Sin[] terms*)

sint=cost/.Cos[a_]:>sin[a];

(＊expand general result for order＊)

r0＝product[Cos[2*Pi/sx(ix-1)(t-vars[[j+1+order]])]+

vars[[j+1]]Cos[2*Pi/sx(ix+1)(t-vars[[j+1+order]])],

{j,1,order}];

r1=TrigReduce[r0];

r2=r1/.Cos[a_]:>Cos[Collect[Together[a],t]];

r3=r2/.Cos[b_t+c_]->Cos[b*t]*Cos[c]-Sin[b*t]*Sin[c];

r4=r3/.{Cos[a_]:>Cos[Apart[a]], Sin[a_]:>Sin[Apart[a]]};

r5=Table[{Coefficient[r4, cost[[j]]], Coefficient[r4, sint[[j]]]},

{j,1,len}];

Return[{terms, vars[[1]]*r5}];

]

WfmColor[list_, color_]:=

Module[

{max, r, g, b, data},

r = color[[1]];

g = color[[2]];

b = color[[3]];

data = Transpose[list][[2]];

max = Max[data];

Return[ListPlot[list, {PlotJoined->True,

GridLines->Automatic,

PlotRange->{max-150, max+5},

ImageSize->600,

Plotstyle->RGBColor[r,g,b]}]];

]

EndPackage[]

Page30

Claims (4)

1. calibration system comprises:

First analog signal generator produces first analog signal;

Second analog signal generator produces second analog signal;

Analog adder is connected to said first analog signal generator and said second analog signal generator, be used for said first analog signal and said second analog signal mutually adduction simulation output is provided;

ADC has the input and numeral output that are connected to said simulation output;

Linearity corrector with the input that is connected to said numeral output is realized the filter product is sued for peace;

Acquisition memory is connected to said numeral output; With

Processor is used to control said first analog signal generator and said second analog signal generator, is connected to read said acquisition memory data, also is connected to said linearity corrector so that be said filter product programming filter factor;

Wherein said first analog signal generator produces the primary sinusoid, and said second analog signal generator produces second sine wave; And

The wherein said primary sinusoid and said second sine wave have the frequency corresponding to the coprime multiple of the frequency of base data record length.

2. calibration system as claimed in claim 1; The rated frequency index wherein is provided; Said rated frequency index is the even multiple corresponding to the frequency of said base data record length; The said primary sinusoid subtracts one corresponding to said rated frequency index, and said second sine wave adds one corresponding to said rated frequency index.

3. calibration system as claimed in claim 1, wherein said processor comprises the special circuit that is used to seek the target function minimum value.

4. calibration system as claimed in claim 1, wherein said processor are that runs software program is to seek the general processor of target function minimum value.

Images