Real-time battery impedance spectrum is acquired using a one-time record. Fast Summation Transformation (FST) is a parallel method of acquiring a real-time battery impedance spectrum using a one-time record that enables battery diagnostics. An excitation current to a battery is a sum of equal amplitude sine waves of frequencies that are octave harmonics spread over a range of interest. A sample frequency is also octave and harmonically related to all frequencies in the sum. The time profile of this signal has a duration that is a few periods of the lowest frequency. The voltage response of the battery, average deleted, is the impedance of the battery in the time domain. Since the excitation frequencies are known and octave and harmonically related, a simple algorithm, FST, processes the time record by rectifying relative to the sine and cosine of each frequency. Another algorithm yields real and imaginary components for each frequency.

Claim:

The invention claimed is:

1. A method for detecting function of a unit under test by measuring frequency response the method comprising; the acts of: setting measurement conditions comprisingchoosing a Root Mean Square (RMS) level of current, identifying a start frequency, and choosing a period of a lowest frequency, selecting a number of octave harmonic frequencies relative to the start frequency over which the function of the unit undertest will be tested, wherein the number of octave harmonic frequencies is one or more; assembling an excitation time record including a sum of sinusoids of the number of octave harmonic frequencies and a duration of the period of the lowest frequency; conditioning the excitation time record to be compatible with the measurement conditions; exciting the unit under test with the excitation time record and simultaneously capturing a resulting response time record with a data acquisition system having asample frequency that is octave related and harmonically related to the frequencies within the sum of sinusoids and at least four times the highest frequency in the sum of sinusoids; rectifying the response time record relative to a sine at eachfrequency of the octave harmonic frequencies with a first sine rectification sum m.sub.11, comprising: .times..times..times..times..function..times..times..DELTA..times..times.- .times..times..times..times..function..times..times..DELTA..times..times. ##EQU00025## Where: TR is a captured response time record with discrete time index n and duration of one period of a first frequency f.sub.1 N.sub.i is a number of sample points for a period of the first frequency f.sub.1 wherein for each sum of sinusoidfrequencies f.sub.i, N.sub.i is constrained as log.sub.2 (N.sub.i) and is an integer greater than 1 and wherein the frequency of interest is given as: .times..DELTA..times..times. ##EQU00026## wherein .DELTA.t is a sample period summing the sinerectification sum for each frequency of interest of the octave harmonic frequencies to create a sine-rectified response time record sum, normalizing the sine-rectified response time record sum to a number of periods known to be within the sine-rectifiedresponse time record sum for each frequency, and storing a normalized sine sum result, m.sub.1,i, where i is an i.sup.th frequency of interest; rectifying the response time record relative to a cosine at each frequency of the octave harmonic frequencieswith a first cosine rectification sum m.sub.21, comprising: .times..times..times..times..times..times..times..DELTA..times..times..ti- mes..times..times..times..times..function..times..times..DELTA..times..tim-es..times..times..times..times..times..function..times..times..DELTA..time- s..times. ##EQU00027## summing the cosine rectification sum for each frequency of interest of the octave harmonic frequencies to create a cosine-rectified response time recordsum, normalizing the cosine-rectified response time record sum to the number of periods known to be within the cosine-rectified response time record sum for each frequency, and storing a normalized cosine sum result, m.sub.2,i, where i is the i.sup.thfrequency of interest; determining a real and an imaginary form of a magnitude and a phase of each of the octave harmonic frequencies: .times..times..times..PHI..times..times..times..PHI..times..times..times.-.times..times..times..times..times..function. ##EQU00028## Where: V.sub.i is a magnitude of the voltage for the frequency f.sub.i .phi..sub.i is a phase of the voltage for the frequency f.sub.i m.sub.1,i, m.sub.2,i are rectification sums for thefrequency f.sub.i .times..function..times..pi..times..times..function..times..pi..times..ti- mes..function..times..pi..times..times..function..times..pi..times..times.- .function..times..pi..times..times..times..function..times..pi..times..tim-es..times..function..times..pi..times..times..function..times..pi..times..- times..times..function..times..pi..times..times..times..function..times..p- i..times. ##EQU00029## and, obtaining the magnitude and the phase via a standard rectangular to polarconversion: .times..times..times..times..times..PHI..times..times..times..PHI..times.- .times..PHI..function..times. ##EQU00030## assembling the magnitude and the phase at each of the octave harmonic frequencies to determine the frequency response.

2. The method of claim 1, wherein rectifying the response time record relative to the sine at each of the octave harmonic frequencies is performed for signals that pass through zero at a discrete time step, wherein N is a number of samples overthe period: .function.<<<< ##EQU00031##

3. The method of claim 1, wherein rectifying the response time record relative to the cosine at each of the octave harmonic frequencies is performed for signals that pass through zero at a discrete time step, wherein N is a number of samplesover the period: .function..ltoreq.<<<.times..times..times.<< ##EQU00032##

4. The method of claim 1, further comprising the acts of: repeating the original acts of; selecting the number of octave harmonic frequencies, assembling the excitation time record, conditioning the excitation time record, exciting the unitunder test, rectifying the response time record relative to the sine at each frequency, summing the sine rectification sum for each frequency of interest, rectifying the response time record relative to the cosine at each frequency, summing the cosinerectification sum for each frequency of interest, determining the real and the imaginary form, and assembling the magnitude and the phase, at least once, wherein the octave harmonic frequencies are shifted in the repeated acts relative to the originalacts; and interleaving the frequency response of the original acts with the frequency responses of the repeated acts to obtain a higher resolution frequency response.

5. The method of claim 1, wherein: rectifying the response time record relative to a sine at each frequency of the octave harmonic frequencies includes a second sine rectification sum m.sub.12, comprising; .times..times..function..times..times..DELTA..times..times..times..functi- on..times..times..DELTA..times..times..times..function..times..times..DELT- A..times..times..times..times..function..times..times..DELTA..times..times- . ##EQU00033####EQU00034## and rectifying the response time record relative to a cosine at each frequency of the octave harmonic frequencies includes a second cosine rectification sum m.sub.22, comprising; .times..times..function..times..times..DELTA..times..times..times..times.- .function..times..times..DELTA..times..times..times..times..function..time- s..times..DELTA..times..times..times..function..times..times..DELTA..times-..times..times..times..function..times..times..DELTA..times..times..times.- .times..times..function..times..times..DELTA..times..times..times. ##EQU00035## .times. ##EQU00035.2## where: N.sub.2 is a number of sample points in a second frequency at anoctave harmonic of the first frequency and 2N.sub.2=N.sub.1.

6. The method of claim 5, wherein: rectifying the response time record relative to a sine at each frequency of the octave harmonic frequencies includes a third sine rectification sum m.sub.13, comprising; .times..times..function..times..times..DELTA..times..times..times..functi- on..times..times..DELTA..times..times..times..function..times..times..DELT- A..times..times..times..times..times..times..times..DELTA..times..times..t-imes..times..times..function..times..times..DELTA..times..times..times..ti- mes..times..function..times..times..DELTA..times..times..times..times..tim- es..function..times..times..DELTA..times..times..times..times..times..func-tion..times..times..DELTA..times..times. ##EQU00036## ##EQU00037## and rectifying the response time record relative to a cosine at each frequency of the octave harmonic frequencies includes a third cosine rectification sum m.sub.23, comprising,comprising: .times..times..function..times..times..DELTA..times..times..times..times.- .function..times..times..DELTA..times..times..times..times..function..time- s..times..DELTA..times..times..times..function..times..times..DELTA..times-..times..times..times..function..times..times..DELTA..times..times..times.- .times..times..times..times..times..DELTA..times..times..times..times..tim- es..times..times..times..times..DELTA..times..times..times..times..times..-times..times..times..times..DELTA..times..times..times..times..times..time- s..times..times..times..DELTA..times..times..times..times..times..times..t- imes..times..DELTA..times..times..times..times..times..times..times..times-..times..DELTA..times..times..times..times..times..times..function..times.- .times..DELTA..times..times..times. ##EQU00038## .times. ##EQU00038.2## where: N.sub.3 is a number of sample points in a third frequency at an octave harmonic of the firstfrequency and 4N.sub.3=N.sub.1.

7. The method of claim 6, comprising: rectifying the response time record relative to a sine at each frequency of the octave harmonic frequencies includes at least one additional sine rectification sum m.sub.1n; and rectifying the responsetime record relative to a cosine at each frequency of the octave harmonic frequencies includes at least one additional cosine rectification sum m.sub.2n; wherein a last frequency includes 2.sup.n-1 periods within a period of the first frequency.

8. A method of calibrating phase response of a captured time record for a unit under test by measuring frequency response of a reference standard the method comprising the acts of: setting measurement conditions comprising choosing a Root MeanSquare (RMS) level of current, identifying a start frequency, and choosing a period of a lowest frequency; determining a theoretical frequency response for the reference standard; selecting a number of frequencies and a frequency spread over which thereference standard will be tested, wherein the frequency spread is by octaves; assembling a phase-shifted excitation time record including a sum of sinusoids of the number of frequencies with a pre-selected phase shift, wherein a duration of theexcitation time record is greater than or equal to one period of a lowest of the number of frequencies; conditioning the phase-shifted excitation time record to be compatible with the measurement conditions; exciting the reference standard with thephase-shifted excitation time record and simultaneously capturing a response time record with a data acquisition system having an appropriate sample frequency; processing the response time record to obtain estimated magnitude and estimated phase at eachof the number of frequencies to determine the frequency response; wherein, phase calibration is achieved through a selection of multiple known phase-shifts until enough points are obtained to determine phase calibration constants for each of the numberof frequencies using a regression analysis.

9. The method of claim 8, wherein the act of processing the response time record includes: rectifying the response time record relative to a sine at each of the number of frequencies, with a first sine rectification sum m.sub.11, comprising:.times..times..times..times..function..times..times..DELTA..times..times.- .times..times..times..times..function..times..times..DELTA..times..times. ##EQU00039## Where: TR is a captured voltage time record with discrete time index n and duration of oneperiod of a first frequency f.sub.1 N.sub.1 is a number of sample points for a first frequency f.sub.1 wherein for each sum of sinusoid frequencies f.sub.i, N.sub.i is constrained as log.sub.2 (N.sub.i) and is an integer greater than 1 and wherein thefrequency of interest is given as: .times..DELTA..times..times. ##EQU00040## wherein .DELTA.t is a sample period summing the sine rectification sum for each frequency of interest of the number of frequencies to create a sine-rectified response timerecord sum, normalizing the sine-rectified response time record sum to a number of periods known to be within the sine-rectified response time record sum for each frequency, and storing a normalized sine sum result, m.sub.1,i, where i is the i.sup.thfrequency of interest; rectifying the response time record relative to the cosine at each of the number of frequencies, with a first cosine rectification sum m.sub.21, comprising:.times..times..times..times..times..times..times..DELTA..times..times..ti- mes..times..times..times..times..function..times..times..DELTA..times..tim- es..times..times..times..times..times..function..times..times..DELTA..time- s..times. ##EQU00041##summing the cosine rectification sum for each of the number of frequencies to create a cosine-rectified response time record sum, normalizing the cosine-rectified response time record sum to the number of periods known to be within the cosine-rectifiedtime record sum for each frequency, and storing a normalized cosine sum result, m.sub.2,i, where i is an i.sup.th frequency of interest; determining the real and imaginary form of the magnitude and phase of each octave harmonic frequency:.times..times..times..PHI..times..times..times..PHI..times..times..times.- .times..times..times..times..times..function. ##EQU00042## Where: V.sub.i is a magnitude of the voltage for the frequency f.sub.i .phi..sub.i is a phase of the voltage for thefrequency f.sub.i m.sub.1,i, m.sub.2,i are rectification sums for the frequency f.sub.i .times..function..times..pi..times..times..function..times..pi..times..ti- mes..function..times..pi..times..times..function..times..pi..times..times.-.function..times..pi..times..times..times..function..times..pi..times..tim- es..times..function..times..pi..times. ##EQU00043## .times..times..function..times..times..pi..times..times..times..function.-.times..times..pi..times..times..times..function..times..times..pi..times. ##EQU00044## and, obtaining the magnitude and the phase via a standard rectangular to polar conversion: .times..times..times..times..times..times..PHI..times..times..times..time-s..PHI..times..times..times..times..times..PHI..function..times..times. ##EQU00045##

10. The method of claim 9, wherein: rectifying the response time record relative to a sine at each of the number of frequencies includes a second sine rectification sum m.sub.12, comprising:.times..times..function..times..times..DELTA..times..times..times..functi- on..times..times..DELTA..times..times..times..times..times..times..DELTA..- times..times..times..times..function..times..times..DELTA..times..times. ##EQU00046## ##EQU00047## andrectifying the response time record relative to a cosine at each of the number of frequencies includes a second cosine rectification sum m.sub.22, comprising: .times..times..function..times..times..DELTA..times..times..times..times.-.function..times..times..DELTA..times..times..times..times..function..time- s..times..DELTA..times..times..times..function..times..times..DELTA..times- ..times..times..times..function..times..times..DELTA..times..times..times.-.times..times..function..times..times..DELTA..times..times. ##EQU00048## ##EQU00049## where: N.sub.2 is a number of sample points in a second frequency at an octave harmonic of the first frequency and 2N.sub.2=N.sub.1.

11. The method of claim 10, wherein: rectifying the response time record relative to a sine at each of the number of frequencies includes a third sine rectification sum m.sub.13, comprising:.times..times..function..times..times..DELTA..times..times..times..functi- on..times..times..DELTA..times..times..times..function..times..times..DELT- A..times..times..times..times..times..times..times..DELTA..times..times..t-imes..times..times..function..times..times..DELTA..times..times..times..ti- mes..times..function..times..times..DELTA..times..times..times..times..tim- es..function..times..times..DELTA..times..times..times..times..times..func-tion..times..times..DELTA..times..times. ##EQU00050## ##EQU00051## and rectifying the response time record relative to a cosine at each of the number of frequencies includes a third cosine rectification sum m.sub.23, comprising, comprising:.times..times..function..times..times..DELTA..times..times..times..times.- .function..times..times..DELTA..times..times..times..times..function..time- s..times..DELTA..times..times..times..function..times..times..DELTA..times-..times..times..times..function..times..times..DELTA..times..times..times.- .times..times..times..times..times..DELTA..times..times..times..times..tim- es..times..times..times..DELTA..times..times..times..times..times..times..-times..times..times..DELTA..times..times..times..times..times..times..time- s..times..times..DELTA..times..times..times..times..times..times..times..t- imes..DELTA..times..times..times..times..times..times..times..times..times-..DELTA..times..times..times..times..times..times..function..times..times.- .DELTA..times..times..times. ##EQU00052## .times. ##EQU00052.2## where: N.sub.3 is a number of sample points in a third frequency at an octave harmonic of the first frequencyand 4N.sub.3=N.sub.1.

12. The method of claim 11, wherein: rectifying the response time record relative to a sine at each of the number of frequencies includes at least one additional sine rectification sum m.sub.1n; and rectifying the response time record relativeto a cosine at each of the number of frequencies includes at least one additional sine rectification sum m.sub.2n; wherein a last frequency includes 2.sup.n-1 periods within a period of the first frequency.

Description:

BACKGROUND OF THE INVENTION

Electrochemical Impedance Measurement Systems use the Bode analysis technique to characterize an impedance of an electrochemical process. It is a well-established and proven technique. A battery being evaluated is excited with a current thatis a single frequency and its response is measured. The process is repeated over a range of frequencies of interest until the spectrum of the impedance is obtained. The method is effective but time consuming, as the process is serial. A parallelapproach using band width limited noise as an excitation current can obtain the same information in less time. The system response to the noise is processed via correlation and Fast Fourier Transform (FFT) algorithms and many such responses areaveraged. The result is the spectrum of response over the desired frequency range. The averaging of many responses also makes this process somewhat serial. Another technique assembles the current noise waveform from a sum of sinusoids, each at adifferent frequency. The system response as a time record is acquired and processed with the FFT algorithm. To reduce noise, multiple time records of waveforms are processed and their resultant spectra averaged. This process is also serial.

There remains a need for real-time acquisition of battery impedance for control and diagnostics over a limited frequency range. This method of acquisition should be a true parallel approach that uses a single time record of battery responsewith a duration compatible with a real-time control process.

BRIEF SUMMARY OF THE INVENTION

The invention involves using a parallel approach to analyze battery impedance or other system functions. A number of frequencies are selected over which the battery is to be tested. These frequencies are assembled into an Excitation TimeRecord (ETR) that is the Sum of the Sinusoids (SOS) of the frequencies and the length of such periods of the lowest of the frequencies. The ETR is conditioned to be compatible with the battery. The battery is then excited with the ETR and a ResponseTime Record (RTR) is captured. The RTR is then synchronized to the ETR and processed by a series of equations to obtain frequency response.

In one preferred embodiment, the RTR is processed to obtain estimated frequency components of magnitude and phase for one of the selected frequencies. Processing is repeated to obtain estimated frequency components for each selected frequency. Frequency components are reassembled to obtain an Estimated Time Record (ETR). The ETR is subtracted from the captured RTR to get an error. The error is minimized to achieve the frequency response estimate. Error is minimized using CompensatedSynchronous Detection (CSD) using a CSD algorithm, which can be implemented by a neural network.

In another preferred embodiment, all excitation frequencies of the SOS are harmonics by powers of two. The sample period likewise is a power of two with all the SOS frequencies. The RTR is rectified relative to a square wave and a 90-degreeshifted square wave of one of the SOS frequencies. Integrating the processed RTR results in an "in phase" and "quadrature" sum that is easily processed to yield the magnitude and phase shift of the desired frequency components. Frequency components areassembled to obtain frequency response.

The subject method allows a parallel implementation for swept frequency measurements to be made utilizing a composite signal of a single time record that greatly reduces testing time without a significant loss of accuracy.

FIG. 30 is a flow chart showing another preferred method of the subject invention.

DETAILED DESCRIPTION OF THE INVENTION

The method of the subject invention allows for real time estimation of a battery's impedance spectrum. The shift of a battery's impedance spectrum strongly correlates to the health of the battery. Therefore, the subject method provides in-situdiagnostics for state-of-health estimation of the battery, which is critical for enhancing the overall application's reliability. The subject method measures a frequency response of a unit under test, for example, a battery. The battery under test isexcited by the sum of sinusoids of a number of test frequencies. A response time record is captured, then processed, to obtain estimates of frequency components for each of the number of frequencies. Estimated frequency components are assembled toachieve a frequency response.

A system function can be identified over a limited number of specific frequencies at step 10 as shown in the flow chart of FIG. 23. The desired frequencies are assembled as an excitation time record 12 that is a sum of those sinusoids and havea length of several periods of the lowest frequency. The time step selected must be compatible with Shannon's sampling constraints for the highest frequency component. The individual waveforms should be sine waves of equal amplitude but alternatingsigns for a phase shift of 180 degrees between the components. Alternating the 180-degree phase shift will minimize a start-up transient. The Root Mean Square (RMS) and the rogue wave peak (sum of the absolute values of all component peaks) of theassembled time record must be compatible with the system being excited and the Data Acquisition System (DAS) that will capture the response.

The Excitation Time Record (ETR), as a current for impedance identification or voltage for system function identification, is signal conditioned 14 to be compatible with the Unit Under Test (UUT). As part of the signal conditioning,anti-aliasing filters ensure that all frequencies generated by the digital-to-analog conversion process other than the intended frequencies are suppressed. The UUT is excited 16 by the ETR and a time record of the UUT response is captured by the DAS. The UUT Response Time Record (URTR) is synchronized to, and the same length as, the ETR.

A preferred embodiment of the processing of the response time record is described below. In order to fit steady-state sinusoidal response assumptions, a preselected number of data points R at the beginning of the URTR must be discarded. Ingeneral, the sum of those data points total to a time that is larger than the transient response time of the UUT at the front end of the ETR. The UUT corrected time response is referred to as the CTR.

In order to fit steady-state sinusoidal response assumptions, a preselected number of data points, R, at the beginning of the URTR must be discarded. In general, the sum of those data points total to a time that is larger than the transientresponse time of the UUT at the front end of the ETR. The UUT Corrected Time Response is referred to as the CTR.

The first estimate of components, magnitude, and phase of the frequency response is made by processing 18 the URTR via Equations 3 through 10. It is important that a zero mean be established prior to processing with Equations 3 to 10 (seebelow).

The core of this whole concept is that an estimate of the UUT Corrected Time Response, the CTR, is made by reassembling 20 the CTR using the estimates of the individual frequency components with the same time step and then discarding the first Rtime steps to become the ECTR. The difference between the CTR and the ECTR is an error 22 and minimizing 24 this error will increase the accuracy of the frequency response estimates.

A first approach to minimizing the error between the CTR and the ECTR is Compensated Synchronous Detection (CSD). The CSD algorithm synthesizes a residual time record of the original time record using the magnitudes of the in-phase andquadrature components for each frequency, except the one to be detected. This synthesized residual is then subtracted from the original time record. The resulting compensated time record is processed with synchronous detection and a new compensatedestimate of the response at the detection frequency is obtained. Since all of the other components in this compensated time record are suppressed, the error from leakage at those other frequencies is less. This process is repeated for each of thefrequencies. Assembling the residual time record and generating the compensated time record are illustrated by Equations 11 and 12 (see below).

Another approach to minimize the error between the CTR and the ECTR is to use a neural network. A first estimate of the component magnitudes and phases is made as described for the CSD technique. Those values are stored and the ECTR iscalculated. This signal is then subtracted from the CTR to produce a response residual. The synchronous detection is then performed upon this residual and the component magnitudes and phases are again stored. These components are then used toreconstruct an estimate of the residual signal. This estimated residual signal is subtracted from the initial residual signal to produce a residual, residual signal. This is then synchronously detected and the loop starts again. This is repeated asmany times as desired, each time with the resultant components being stored. The assumption is that there is a functional relationship between these resultant components and a true system response. A neural network is then used to determine thisrelationship. The previously stored results become the test dataset for the neural network. The network has been trained previously on a similar unit and a known ideal response (e.g., battery impedance measured using electrochemical impedancespectroscopy). The output of the network is the estimate of the response.

A further step of the subject method is to shift the complete set of desired frequencies and to repeat the whole process. This step could be repeated many times with different shifts to develop a high-resolution frequency response that forbattery impedance could be comparable to that provided by electrochemical impedance spectroscopy. Thus, the subject method can provide capability for both limited frequency response in real-time or high-resolution frequency response not in real-time forperiodic system in-depth diagnostics.

The system of the subject invention is based on the following theoretical design. The Unit Under Test (UUT) is excited with a limited sum of sinusoids, each at a different frequency that is spread over the range of interest. The magnitude,frequency and phase of each sinusoid making up the sum are known. If a total response of the system is measured via a sample data system at an acceptable sample rate and an adequate duration time record is acquired, then a simple algorithm that uses theknown magnitude, frequency and phase of each individual sinusoid will process the single time record. This analysis will obtain the true Bode response at the selected frequencies spread over the range of interest all in parallel. The followingsynchronous detection analysis is the basis of this simple algorithm. The reference waveform is chosen as a sine, as at time zero everything will be at zero.

Each component magnitude and phase of the system response at all the excitation frequencies can be obtained via the following synchronous detection analysis. This analysis quantifies the response at the k.sup.th radian frequency .omega..sub.kwith the "in phase" and "quadrature" response. The analysis incorporates a feature of discarding a preselected number of points R at the beginning of the system response in order to meet the assumption of steady-state sinusoidal response. Additionally,for most applications, prior to processing the data, the mean of the acquired time record should be computed and deleted. The presence of a non-zero mean could corrupt an estimate of the lowest frequency component.

Again, the summation over j for infinite time record averages everything to zero except the final result of Equation 8. Equations 5 and 8 can be combined to give magnitude and phase for the k.sup.th frequency response.

Equations 9 and 10 give the final results for synchronous detection. This process is repeated for all M of the excitation sinusoids. As stated, the results depend on the time record being infinite in duration. This is not the case and, thus,the results are in error by leakage from the other components. This error can be reduced without increasing the time record using the following algorithm for Compensated Synchronous Detection (CSD).

The CSD algorithm synthesizes a residual time record of the original time record using the magnitudes of the in-phase and quadrature components for each frequency except the one to be detected. This synthesized residual is then subtracted fromthe original time record. The resulting compensated time record is then processed with synchronous detection, and a new compensated estimate of the response at the detection frequency is obtained. Since all of the other components in this compensatedtime record are suppressed, the error from leakage at those other frequencies will be less. This process is repeated for each of the frequencies. Assembling a residual time record and generating the compensated time record are illustrated by Equations11 and 12.

F.sub.P is the estimated in phase amplitude response at the p.sup.th frequency

F.sub.Qp is the estimated quadrature amplitude response at the p.sup.th frequency

.omega..sub.p is the radian frequency of the p.sup.th sinusoid

.DELTA.t is the time step of the data system

N is the number of points of the output time record

M is the number of different sinusoids of the excitation function

R is the number of points of the output time record that are discarded

The synchronous detection algorithm described by Equations 1 through 8 is applied to the compensated time record of Equation 12 and better estimates of F.sub.P and F.sub.Qp are obtained. This process can be repeated again until the totaldifference between a completely synthesized time record response and the original time record is minimized.

The following examples are offered to illustrate the method of the subject invention and should not be construed as limiting.

Example 1

Analytical Testing on a Sum of Sines

The CSD algorithm was evaluated using a simple signal that was assembled from a finite sum of equal amplitude sine waves (Sum of Sines (SOS)) with frequencies distributed logarithmically over a limited range. The objective of the analysis wasto assess how well the CSD algorithm could pick out an amplitude for each component.

To check out the concept analytically, a MATLAB.RTM. matrix calculation computer software code was written that was a logarithmic mix of 5 equal unity amplitude frequencies (5.sup.0.5, 5.sup.1, 5.sup.1.5, 5.sup.2, 5.sup.2.5 Hz). The acquiredtime record was set to 10 periods of the lowest frequency and the time step was set to 1/10 of the period of the highest frequency. As per Equation 9, error-free detection should estimate the amplitude of 0.5 for each component. Table 1 gives anestimate for a first pass or simple synchronous detection and a second pass, the CSD algorithm. The MATLAB.RTM. matrix calculation computer software code for this analysis is disclosed by Morrison et al., 2005, Algorithms as MATLAB Code for Real timeEstimation of Battery Impedance.

As seen in Table 1 analytically, the compensation technique does appear to work as the error for the second pass is much reduced. This initial check of the concept was applied to detect amplitude only and not phase detection. The signal is asum of equal amplitude sine waves being decomposed into the individual components by the algorithm.

Example 2

Analytical Testing of a Low-Pass Filter

A recursive model of a second order low-pass function was excited with an SOS input signal. The CSD algorithm was then used to estimate the frequency response at each of the specific frequencies making up the SOS.

A spread of 13 specific frequencies was chosen that were spaced in 1/4 decade steps starting from 0.1 Hz up to 100 Hz. Using these frequencies, a mix of equal unity amplitude sine waves was created. This range of frequencies was picked asresearch performed at the Idaho National Laboratory with batteries is typically over this frequency spread. The signal was discretized with a time step that was 10% of the period of the highest frequency. The length of the time record was set at 10periods of the lowest frequency.

A recursive model of a second order Butterworth low-pass function was developed. A center frequency was set at the middle of the SOS frequency spread. A filter response to the SOS input time profile was computed.

Finally, the CSD code was developed to estimate the filter frequency response from the time profile generated by the recursive model code. That code is the implementation of Equations 1 through 12. Additionally, the implementation has theability to discard a number of user-selected points at the beginning of the time profile such that the remaining data better fits the assumption of "steady-state" response. The following 6 figures are MATLAB.RTM. matrix calculation computer softwarediscrete plots of an ideal frequency response, an uncompensated frequency response, and finally, a compensated frequency response. All 3 have a magnitude plotted against a Log of frequency. FIGS. 1 and 2 are an ideal magnitude and phase response. FIG.3 is an uncompensated magnitude response, while FIG. 4 is a compensated magnitude response. The improvement seen by comparing the last two figures is clear FIG. 5 is the compensated phase response. Error at the higher frequencies for the compensatedresult is likely due to processing a small signal at the high frequencies relative to larger signals at the lower frequencies. The net result is that the one-time record yields a limited number of points for both magnitude and phase. This techniqueshows promise for real-time applications. The MATLAB.RTM. matrix calculation computer software code for this analysis is disclosed by Morrison et al., 2005, Algorithms as MATLAB Code for Real time Estimation of Battery Impedance. The next approachwill apply the concept analytically in an attempt to identify the impedance of a computer model of a battery.

Example 3

Analytical Testing of CSD with a Battery Model

The CSD algorithm was evaluated analytically via a computer simulation of the detection of the impedance of the Lumped Parameter Model of a battery (LPM) that was developed by the Idaho National Laboratory (INL) (see, FreedomCAR Battery TestManual, 2003). A computer model for the LPM that will simulate battery voltage response to an arbitrary battery current was also developed at INL by Fenton et al., 2005. The voltage response of the model normalized to the current in the frequencydomain will be the battery impedance. The equivalent circuit for the LPM with parameter identification is shown in FIG. 6.

The LPM was excited with a current source I.sub.IN that was an SOS, and the CSD algorithm was used to identify the impedance seen looking into the LPM over a limited range of discrete frequencies. It should be noted that the polarity of thevoltage response was defined as negative because the SOS excitation current was a discharge (negative relative to impedance). The CSD algorithm was used to obtain the frequency response of the LPM. The CSD response magnitude and phase were compared tothe ideal response. Initially, the algorithm failed to match ideal impedance. The response of a battery terminal voltage to a discharge SOS current signal will contain a DC term caused by the DC battery voltage. Synchronous detection of any specificfrequency in the response will cause a noise frequency in the resultant spectrum at the frequency being detected, at an amplitude of the product of the DC battery voltage, and an amplitude of the detection signal. That noise amplitude will be largerelative to the signal being detected. Averaging enough cycles in the resultant time record will reject this noise. However, for a real-time application, the length of the time record needs to be short and not long. The problem was resolved when themean was deleted from the prediction response of the LPM. The number of frequency lines was set to 13 and were logarithmically spread from 0.01 Hz to 1.0 Hz. The time record was set to 10 periods of the lowest frequency. In the CSD algorithm, no datapoints were discarded. Table 2 gives analysis specifics with LPM data that is typical for some lithium batteries that INL had recently tested. INL performed the testing per methods in the FreedomCAR Battery Test Manual, 2003.

The following 7 figures are the plots of the analysis results. FIG. 7 is the time record of the SOS current signal. FIG. 8 is the time record of the LPM voltage response. FIGS. 9 and 10 are the LPM ideal impedance magnitude and phase. FIGS.11 and 12 are the uncompensated and the compensated magnitude response. FIG. 13 is the compensated phase response.

The number of frequency lines was increased to 25 and everything else was left the same. The following 7 figures illustrate the results. FIG. 14 is the time record of the SOS current signal. FIG. 15 is the time record of the LPM voltageresponse. FIGS. 16 and 17 are the LPM ideal impedance magnitude and phase. FIGS. 18 and 19 are the uncompensated and the compensated magnitude response. FIG. 20 is the compensated phase response.

Additional cases run showed that as the length of the acquired time record in the number of periods of the lowest frequency gets cut back, the number of frequency lines that can be resolved without a big increase in error must be cut back. Forexample, 5 periods with 25 lines gave terrible results, but 5 periods with 13 lines was fine. These positive results are only analytical. Nevertheless, they offer promise of positive expected performance when applied to a physical system. All theMATLAB.RTM. matrix calculation computer codes for this analysis are given by Morrison, 2005, Algorithms as MATLAB.RTM. Code for Real time Estimation of Battery Impedance.

Example 4

Neural Network Enhanced Synchronous Detection

In order to improve the accuracy of the CSD, studies were conducted upon neural network enhancement of the detection of the individual frequency components of the response of a linear system to an SOS input signal. The concept is very similarto the CSD, with some slight changes and a neural network output layer. For a second order low-pass filter, ordinary synchronous detection performed on the filter response showed a mean squared error (MSE) of 2.6.times.10.sup.-3. The compensatedsynchronous detection technique displayed an MSE of 1.6.times.10.sup.-3. The neural network enhanced synchronous detection showed an MSE of 0.2.times.10.sup.-3. Results for the lumped parameter model of a lithium ion battery were similar.

The theory of neural network enhanced synchronous detection is based on a classical synchronous detection and an inherent error associated with time records of finite length.

Given an input signal comprised of a sum of N sinusoids:

.function..times..function..omega..times. ##EQU00007## the output of the system would be:

.function..times..function..omega..times. ##EQU00008##

Being that the input and output are sinusoids, they are assumed to have started at t=-.infin. and continues to t=.infin.. In reality, however, this is not the case. First, the time record of the signal is finite in length; second, it issampled. Given a sampling frequency of at least the Nyquist frequency (twice the highest frequency in the signal), the signal can be reconstructed without error. This does not, however, rectify the finite time length of the signal. Because of this,errors enter into the synchronously detected frequency components.

Synchronous detection involves multiplying the acquired signal by sines and cosines of the desired frequencies and summing the results, as shown below:

The magnitude of a given frequency is obtained as follows: M.sub..omega..sub.n=2 {square root over (.alpha..sub..omega..sub.n.sup.2+.beta..sub..omega..sub.n.sup.2)} and the phase may also be obtained as follows:

.omega..function..alpha..omega..beta..omega. ##EQU00010## where:

y[j]=the sampled signal

.alpha.=the sine component of the response for the frequency .omega..sub.n

.beta.=the cosine component of the response for the frequency .omega..sub.n

.DELTA.t=the sample time step

Notice that the summations are infinite in length. In application, the summation would be from 0 to the length of the recorded signal. If the time record was infinite, then any errors would cancel out. Since our time record is not infinite,errors remain in the calculated response.

It is important to note that an estimate of the original signal can be reconstructed by summing the sine and cosine signals of the different frequencies, where each sine signal is multiplied by its sine response component and each cosine signalis multiplied by its cosine response component. The resulting estimate would be exact if the time record was infinite, but since it is of finite length, the reconstructed signal is not exactly correct. This may be viewed as containing noise.

In the compensated synchronous detection approach, this difference between the actual signal and the reconstructed signal is exploited to filter out all but one of the frequencies and to increase the accuracy of the measurement (see above). InNeural Network Enhanced Synchronous Detection (NNESD), the approach is slightly different. The premise is that the residual left from subtracting the reconstructed signal from the actual signal still contains some useful information. r[i]=y[i]-{tildeover (y)}[i]

For each frequency of interest, the sine and cosine components are detected on the original signal. These are then used to reconstruct the signal and then stored. The reconstructed signal is then subtracted from the original signal to obtainthe residual signal. Synchronous detection is then performed upon this residual signal; once again the sine and cosine components are used to reconstruct the residual signal and are also stored. This second residual signal is then subtracted from thefirst residual signal. This loop continues until a sufficient number (M) of sine and cosine responses are obtained for each frequency. Now, the assumption that is made is that there is some functional relationship between these M responses and the trueresponses. .alpha..sub..omega..sub.n.sub.truth=f(.alpha..sub..omega..sub.n.sub.1,.be- ta..sub..omega..sub.n.sub.1,.alpha..sub..omega..sub.n.sub.2,.beta..sub..om- ega..sub.n.sub.2, . . . ,.alpha..sub..omega..sub.n.sub.M,.beta..sub..omega..sub.n.sub.M).beta..sub..omega..sub.n.sub.truth=f(.alpha..sub..omega..sub.n.sub.1,.bet- a..sub..omega..sub.n.sub.1,.alpha..sub..omega..sub.n.sub.2,.beta..sub..ome- ga..sub.n.sub.2, . . . ,.alpha..sub..omega..sub.n.sub.M,.beta..sub..omega..sub.n.sub.M)

This functional relationship most likely varies from system to system and also is based upon system operating conditions. To deduce what this relationship is for any system would be time consuming. Instead, a generalized regression neuralnetwork is implemented to predict the true response.

A Generalized Regression Neural Network (GRNN) is a form of a radial basis neural network suitable for problems involving regression, function estimation, or prediction of a continuous value, Wasserman, 1993, Alpaydin, 2004. This is in contrastto a very similar network, the Probabilistic Neural Network, which is used for classification applications. Unlike multi-layer perceptrons, which require training, a GRNN is constructed from a training set of example inputs and the correspondingoutputs. The spread of the radial basis function is used to compute the bias. The spread, in effect, controls the spread of the radial basis function. The example inputs are in the form of an m.times.n matrix where each of the m rows represents anobservation and each of the n columns represents a feature, or parameter. The corresponding example output is an m.times.1 vector. The input is said to be n-dimensional. The network is divided into 2 layers. The first, or input layer, consists of theexample inputs. The output layer consists of the example outputs. The network generates its output in the following manner. The geometric distance is calculated between the newly presented input and each of the example inputs in the input layer:

.times..times..times. ##EQU00012## This produces a vector of length m. Each element of the vector is then multiplied by the bias, which is 0.8326/Spread. The vector then goes through the radial basis function. y=e.sup.-x.sup.2 The radialbasis function produces an output that gets closer to 1 as the input gets closer to 0. The resulting vector is a ranking of how close each of the example cases are to the new input. The normalized dot product is then calculated between the vector andthe example output vector. This is the network output.

.times..times..times..times..times..times. ##EQU00013##

The NNESD is based upon the GRNN approach. Analytical testing will provide a preliminary validation of the concept.

The NNESD concept was analytically tested on a 2.sup.nd order Butterworth filter using the MATLAB.RTM. matrix calculating computer software BUTTER function and the INL LPM FreedomCAR Battery Test Manual, 2003, Fenton et al., 2005, for thelithium ion battery. An SOS input signal was used for both cases. For the filter, the component frequencies were varied for each run in order to build up a training set with 10 component lines to be detected. The filter was run 100 times and theoutput and target response was calculated. The data was randomized and half of the data was used to construct the GRNN and the other half was used to verify the network. The mean squared error (MSE) of the predicted value was calculated and the processwas repeated 20 times, shuffling the data each time. The results are shown in FIG. 21. The MSE of both the standard synchronous detection and the CSD are also shown for comparison. The MATLAB.RTM. matrix calculation computer software code used forthis analysis is disclosed by Morrison, 2005, Intelligent Self-Evolving Prognostic Fusion.

The technique was then tested on the LPM FreedomCAR Battery Test Manual, 2003, Fenton et al., 2005. The input parameters to the model were nominally set as per Table 2 and randomly varied by up to 5% each run for 100 runs to generate thedataset. The number of lines, frequency spread and time step were held as per Table 2. The same training and testing scheme that was outlined above was used. The results are shown in FIG. 22. The MATLAB.RTM. matrix calculation computer software codeused for this analysis is given in Morrison, 2005, Intelligent Self-Evolving Prognostic Fusion.

This limited analytical validation has shown that the NNESD concept will significantly reduce the error in the estimate of the frequency components of a given system response signal. This concept allows a parallel implementation for sweptfrequency measurements to be made by utilizing a composite signal of a single time record that greatly reduces testing time without a significant loss of accuracy.

Conclusions

The physical validation of the CSD or NNESD concept will rely heavily on work performed by W. Albrecht in his thesis research: "Battery Complex Impedance Identification with Random Signal Techniques," Albrecht, 2005. In this approach, aNational Instruments data acquisition and processing system was used along with a custom analog conditioning system. The CSD or NNESD algorithm could be installed directly on the custom analog conditioning system. The system software will be CSD orNNESD rather than Noise Identification of Battery Impedance (NIBI). The NIBI approach would acquire about 100 time records of the battery response to noise current. Clearly, a time record would have to be of a length of multiple periods of the lowestfrequency of interest. The CSD or NNESD approach will acquire one time record of a length of multiple periods (exact number is still to be determined) of the lowest frequency of interest. The excitation signal would be generated analytically withsoftware, as in the NIBI system. The analytical signal would be preconditioned with a digital low-pass filter as in the NIBI system. The CSD or NNESD system may require an analog filter prior to the current driver and after the D/A. This analog filterat the current driver could serve as the prime anti-aliasing filter. This system will use the same bias compensation approach to remove most of the DC battery voltage from the acquired signal. Improved noise rejection and increased sensitivity could beachieved if the voltage sensing were upgraded to full differential via a 4-wire system rather than the 2-wire single-ended system of the NIBI. An increase in sensitivity will enable a reduction in the level of the excitation signal required. It isanticipated that the sampled voltage will be processed directly with the CSD or NNESD algorithm. It is also anticipated that a system calibration would be done exactly as the NIBI system by measuring the impedance of the test leads to determine anysystem measurement offset and phase shift.

Currently, calibrating the CSD system for magnitude response is done using current shunts in place of a test battery. Multiple shunts encompassing range of response are measured and the resulting data are used to generate a linear regressioncalibration curve. This technique is also essentially a one-point phase calibration with the shunt, by definition, having zero degrees of phase shift over the frequency range of interest. Since the entire measurement system is a series process, phaseshift in the excitation signal, the test object, the detection amplifier or the processing algorithm all look like phase shift. Thus, during shunt calibration, additional multi-point phase calibrations can be obtained by introducing a calibrated phaseshift into each frequency of the SOS excitation signal and the response processed in the CSD algorithm by using the non-phase shifted SOS signal. A multipoint linear regression magnitude calibration, as well as a multi-point linear regression phasecalibration, is obtained.

The principal attribute of the CSD method is the ability to obtain a limited resolution frequency spectrum of a system function in real time. Resolution of just a dozen frequency lines is traded off to obtain the very desirable feature ofreal-time response. The CSD system is able to reduce cross-talk between adjacent frequency lines over and above what would occur with only synchronous detection. Nevertheless, the number of frequency lines in the single time record must be limited. Ina particularly preferred embodiment, a user can choose to increase resolution at the expense of increased response time with either operation. The limited number of frequency lines in the SOS is adhered to in a single time record. However, for eitheroperation additional time records are run with each frequency shifted by a shift factor and the additional spectra acquired are interleaved to obtain an overall higher resolution spectrum. Thus, a user can make a custom trade-off of spectrum resolutionagainst response time for an application.

The CSD system processes the system response of the SOS excitation via repeated synchronous detection acting on a residual signal whereby cross-talk contributors have been subtracted out. An alternative embodiment of the method of detectingsystem function of the subject invention, as shown in the flow chart of FIG. 30, uses a fast summation algorithm to process the system response of the SOS excitation at step 32. In the Fast Summation Transformation (FST) embodiment, all the frequenciesof the SOS are harmonics by powers of two at step 26. Additionally, the sample period is also a power of two with all the SOS frequencies at step 32. Instead of multiplying the acquired time record by the sine and the cosine of each frequency, the SOSis simply rectified relative to the square wave and the 90-degree shifted square wave of the desired frequency. When the samples of that processed time record are summed, all the octave-related harmonics, other than the frequency of interest will alwayssum to zero. The resulting "In Phase" and "Quadrature" sums can be easily processed to yield the magnitude and phase shift of the desired frequency component at step 34.

The FST algorithm can be used to estimate a battery impedance spectrum at step 36. The excitation Sum of Sines (SOS) signal formation at steps 26, 28, 30, the methodology to apply it as a current signal to a test battery and the capture of thebattery voltage response time record at step 32 is almost identical to the CSD approach. However, the difference is, if the SOS excitation time record contains only sine waves whose frequencies are all related by octave and harmonics, and the sampletime is also octave and harmonically related, then if that captured time record is "rectified" relative to one of the SOS frequencies at step 34, when that transformed time record is summed, it will contain only battery response information relative tothat frequency. To identify a specific battery frequency response in the SOS signal, the response time record is square wave rectified with a phase relationship relative to a sine wave of that frequency at step 36. Then, all the points in the timerecord of that transformed signal are simply totaled up and normalized to the number of periods of that frequency present in the SOS signal. This result becomes a numerical parameter m.sub.1. The process is performed again except the rectificationsquare wave phase is relative to the cos of the frequency of interest. This other result becomes another numerical parameter m.sub.2. For all the other frequencies, except the one of interest, all the samples of those sine waves of the transformedrecord will always total to zero. The amplitude and phase response of that frequency are obtained as per the following relationships. Equation 1 represents the sampled signal component at a specific frequency that is to be detected. The amplitude,V.sub.P and phase, .phi. are the desired information.

In Equation 4, the signal has been rectified relative to the cosine wave of that frequency and again the samples are totaled up. Observe that rectification simply involves changing the sign of the sample values relative to the sine wave orcosine wave timing.

.times..times..function..times..pi..times..PHI..times..times..times..func- tion..times..pi..times..PHI..times..times..times..function..times..pi..tim- es..PHI..times..times..function..times..pi..times..times..times..times..PH-I..function..times..pi..times..times..times..times..PHI..times..times..fun- ction..times..pi..times..times..times..times..PHI..function..times..pi..ti- mes..times..times..times..PHI..times..times..times..times..times..PHI..tim-es..times..function..times..pi..times..times..times..times..pi..times. .times..times..times..PHI..times..times..function..times..pi..times..time- s..function..times..pi..times. .times..times..function..times..pi..times..times..times..times..PHI..func-tion..times..pi..times..times..times..times..PHI..times..times..times..fun- ction..times..pi..times..times..times..times..PHI..times..times..pi..times- ..times..times..times..PHI..times..times..times..function..times..pi..time-s..times..times..times..PHI..times..times..pi..times..times..times..times.- .PHI..times..times..times..times..times..PHI..times..times..function..time- s..pi..times..times..times..function..times..pi..times..times..times..func- tion..times..pi..times. .times..times..times..PHI..times..times..function..times..pi..times..time- s..times..function..times..pi..times..times..times..function..times..pi..t- imes. ##EQU00017## Note that the parameters K.sub.1, K.sub.2, K.sub.3, K.sub.4 are known for eachfrequency and m.sub.1, m.sub.2 are the numerical result of the rectifying algorithm for each frequency. Then, the magnitude and phase at each frequency can be obtained as follows:

For all other signal frequencies present in the rectification process, their values must total to zero whenever the rectification is not their specific frequency. The FST system must preserve the property of orthogonality between the differentfrequencies. Additionally, any noise present, even after rectification, will be mitigated by the summation process. Thus, for this technique, mathematically, cross-talk from adjacent frequencies will be zero and the SOS time record length can be asshort as one period of the lowest frequency. However, the FST method only works if the octave harmonic relationship holds for all frequencies in the SOS including the sample frequency. This ensures that over a period of any frequency present in the SOSthere will always be an even number of samples. The implementation of the rectification functions with discrete signals in a manner that preserves orthogonality between different frequencies will be discussed in the next section.

The algorithm shown above can be realized as computer code. A rectification function is simply a unity amplitude square wave. FST uses two forms at each frequency, one with the phase relationship of a sine wave and the other with the phaserelationship of a cos wave. In a continuous time domain, a perfect unit rectification function makes an instant transition from -1 to +1. In a discrete rectification function, there are two means for transitioning from +1 to -1 and the reverse: passthrough zero at a discrete time step or pass through zero midway between discrete time steps. In the discrete time domain, the implementation of a rectification function must preserve FST orthogonality between different frequencies (i.e., nocross-talk).

The sine waveform of the non-zero rectification function is given as:

.function..ltoreq.<.ltoreq.< ##EQU00019##

Where: N is the number of samples over the period.

The cos waveform of the rectification function is given as:

.function..ltoreq.<.ltoreq.<.times..times..ltoreq.< ##EQU00020## Where: N is the number of samples over the period.

The sine waveform of the rectification function with zeros is given as:

.function.<<<< ##EQU00021## Where: N is the number of samples over the period.

The cos waveform of the rectification function with zeros is given as:

.function..ltoreq.<<<.times..times..times.<< ##EQU00022## Where: N is the number of samples over the period.

There are two means to implement the non-zero rectification function (Equations 9 and 10), but only one means to implement the rectification with zeros (Equations 11 and 12). All means give nearly identical results. For the non-zerorectification function, the more complicated means that will work is to take the average between consecutive sample pairs to obtain the samples to be summed. Specifically, if one averages between points of the sinusoid going around a circle, with thelast point of the period averaged with the first point. The sign of the non-zero rectification function is applied to these averaged samples. As an example, consider an eight point unity amplitude discrete sine and cos signal.

Observe the existence of orthogonality with the sine rectified by the cos, or the cos rectified by the sine. This is an interesting but unnecessary result that is not needed by the FST algorithm.

The second means of rectification for the non-zero rectification function, and the preferred, is the simplest for the FST technique. For rectification, one just changes the signs of the captured time record as per Equations 9 and 10. Considerthe same eight sample discrete sine and cos waves rectified by the second means. First, the sine rectified by the sine, Equation 9: rectification: :FSTss=[0+0.707+1+0.707+0-(-0.707)-(-1)-(-0.707)]=4(0.707)+2

Next, the cos rectified by the sine, Equation 9: FSTcs=[1+0.707+0-0.707-(-1)-(-0.707)-(0)-(0.707)]=2

Now, the sine rectified by the cos, Equation 10 FSTsc=[0+0.707-(1)-0.707-(0)-(-0.707)+(-1)+(-0.707)]=-2

Finally, the cos rectified by the cos, Equation 10: FSTcc=[1+0.707-(0)-(-0.707)-(-1)-(-0.707)+0+0.707]=4(0.707)+2

For the rectification with zeros, Equations 11 and 12, consider the same eight sample discrete sine and cos waves rectified by a second mean. First, the sine rectified by the sine, Equation 11: rectification::FSTss=[0+0.707+1+(0.707)+0-(-0.707)-(-1)-(-0.707)]=4(0.707)+2

Next, the cos rectified by the sine, Equation 11: FSTsc=[0+0.707+0-(0.707)+0-(-0.707)+0-0.707]=0

Now, the sine rectified by the cos, Equation 12 FSTsc=[0+0.707+0-0.707-0-(-0.707)-(0)-(0.707)]=0

Finally, the cos rectified by the cos, Equation 12: FSTcc=[1+0.707+0-(-0.707)-(-1)-(-0.707)+0+0.707]=4(0.707)+2

Again, observe the existence of orthogonality with the sine rectified by the cos or the cos rectified by the sine. This is an interesting but unnecessary result that is not needed by the FST algorithm.

The non-zero sine and cos rectification functions in the preferred non-averaging means, Equations 9 and 10, respectively, when implemented into the FST Equations 1 through 8, preserve orthogonality between frequencies. This is shown via aMATLAB.RTM. matrix calculation computer software code that builds a discrete time record of an SOS where the frequencies are octave harmonics including the sample frequency. All the sinusoids in the SOS start at time zero, have zero phase shift andunity amplitude. The process starts with all waves being sine waves and being rectified with the sine rectification function, Equation 9, for each frequency and then repeat with rectification by the cos rectification function, Equation 10. Then, theSOS is all cos waves and is rectified with the sine rectification function, Equation 9. Cross-talk is always expected to be zero and, thus, there is orthogonality between frequencies. For this evaluation, the SOS is picked to consist of 10 distinctfrequencies, all of which are octave and harmonically related with the highest frequency having 32 samples per period. The time record length is 1 period of the lowest frequency. For Case 1, the SOS is made from discrete sine waves and for Case 2, theSOS is made from discrete cos waves. Orthogonality is proven by deleting an arbitrary frequency from the SOS or SOC and then trying to detect it with the FST Equations 1 through 8. If the result of detection is zero, then there is no cross-talk fromany of the other frequencies and orthogonality is established.

Case 1

SOS is unity amplitude sine waves, 10 octave frequencies with 32 samples per period of the highest frequency. The SOS time record length is 1 period of the lowest frequency. Orthogonality is tested between frequencies by deleting from the SOSan arbitrary frequency and then using FST to detect it. Anything that is detected will be cross-talk corruption. The 3rd frequency was deleted. The results for sine and cos rectification FST are:

As expected, there is no response at the 3rd frequency. In fact, whichever frequency is deleted, there will be no response at that frequency; thus, orthogonality between frequencies is considered valid relative to an SOS.

Case 2

The parameters are the same as in Case 1, except SOS is cosine waves (SOC). Now, if an arbitrary frequency is deleted from the SOC, detection of the deleted frequency measures cross-talk. If there is orthogonality between frequencies, theresult should be zero. The 5th frequency was deleted and the cos and sine results, respectively, are:

Again, as expected, there is no response at the 5th frequency. In fact, whichever frequency is deleted, there will be no response at that frequency; thus, orthogonality between frequencies is considered valid relative to a SOC.

Example 6

Non-Zero Rectification FST Applied to a Battery Model

The method of the subject invention using the FST algorithm with non-zero rectification was applied to a battery model. Specifically, by recursive implementation of the INL Lumped Parameter Model (LPM) (FreedomCAR Battery Test Manual, 2003). LPM parameters were the same as those used in Example 3. Several cases were analyzed and results plotted as FST compared to ideal LPM impedance (classical j.omega. circuit analysis). The plots were Bode (magnitude and phase) and Nyquist. All theseplots are shown in continuous format rather than the discrete format of Example 3. Nevertheless, the individual points as computed are clearly indicated.

Case 1

In the first case analyzed, the SOS started at 0.01 Hz, had 10 frequencies, the time record length was 1 period of the 0.01 Hz, and the highest frequency was set to 4 samples as this is considered a worst-case lower limit. FIGS. 24A-24C showthe impedance spectrum results, respectively, as magnitude vs. frequency, phase vs. frequency and the Nyquist plot of imaginary component vs. real component.

Case 2

It is suspected that the observed error between FST and the ideal response is due to the transient response of the LPM. Thus, this case is the same as Case 1 except the time record was increased to 3 periods of the 0.01 Hz frequency. FIGS.25A-25C show the results in the same format as Case 1.

Case 3

In Case 2, it was observed that with the time record of Case 1 expanded to 3 periods of the lowest frequency, the error was greatly reduced but still present. Case 3 is the same as Case 1 except the time record is set at 2 periods of the lowestfrequency and then the first period is ignored and only the second period is processed. The objective was to delete from the time record the LPM-corrupting transient response. As seen in the following plots, FIGS. 26A-26C in the same format as Case 1,the approach worked as the FST and ideal results overlay each other. This observation is very strong evidence that the LPM transient response is corrupting the non-zero rectification FST algorithm.

Example 7

Analytical Testing of the Rectification Function with Zeros

The second option of a discrete rectification function for transitioning from +1 to -1, and the reverse, and passing through zero at a discrete time step is examined. Equations 11 and 12, when implemented into the FST algorithm and Equations1-8, are shown to preserve orthogonality between frequencies. The sine and cosine waveform of such a function are given, respectively, as Equations 11 and 12. These sine and cosine rectification functions with zeros preserve orthogonality betweenfrequencies shown via a MATLAB.RTM. matrix calculator computer software code that builds a discrete time record of an SOS where the frequencies are octave harmonics including the sample frequency. All the sinusoids in the SOS start at time zero, havezero phase shift and unity amplitude. All sine waves are used initially and are rectified with the sine and cosine rectification functions for each frequency, as shown in Equations 11 and 12. Then, the SOS is all cosine waves, SOC and rectified withthe sine and cosine rectification functions, as shown in Equations 11 and 12. The cross-talk results are expected to always be zero. For this evaluation, the SOS is picked to consist of 10 distinct frequencies, all of which are octave and harmonicallyrelated with the highest frequency having 32 samples per period. The time record length is 1 period of the lowest frequency. For Case 1, the SOS is made from discrete sine waves and, for Case 2, the SOS is made from discrete cos waves. Orthogonalityis proven by deleting an arbitrary frequency from the SOS or SOC and then trying to detect it with the FST Equations 1 through 8. If the result of detection is zero, then there is no cross-talk from any of the other frequencies and orthogonality isestablished.

Case 1

Now, if an arbitrary frequency is deleted from the SOS, cross-talk is measured when trying to detect the arbitrary frequency. If there is orthogonality between frequencies, the result should be zero. The 2nd frequency was deleted and the sineand cos results respectively are:

As expected, there was no response at the 2nd frequency or any other frequency that might be deleted from the SOS, thus, there is orthogonality between frequencies. Additionally, there are 13 orders of magnitude difference between sine and cosrectification results; thus, orthogonality also exists between sine and cosine for an SOS. This feature is of interest but not necessary for the subject method employing the FST algorithm.

Case 2

Parameters are the same as in Case 1, except SOS is cosine waves (SOC). Deleting an arbitrary frequency from the SOC and then trying to detect it measures cross-talk. If there is orthogonality between frequencies, the result should be zero. The 7th frequency was deleted and the cos and sine results, respectively, are:

As expected, there is no response at the 7th frequency or any other frequency that might be deleted from the SOC. Thus, orthogonality also holds for the SOC. Additionally, there are 13 orders of magnitude difference between cos and sinerectification results; thus, orthogonality between sine and cos for an SOC also exists, an interesting but unnecessary feature.

Example 8

Rectification with Zeros FST Applied to a Battery Model

The rectification with zeros version of the FST embodiment was evaluated with the recursive implementation of the INL LPM exactly as in Example 6.

Case 1

In the first case analyzed, the SOS started at 0.01 Hz, had 10 frequencies, the time record length was 1 period of the 0.01 Hz, and the highest frequency had 16 samples. Results shown in FIGS. 27A-27C show the impedance spectrum results,respectively, as magnitude vs. frequency, phase vs. frequency and the Nyquist plot of imaginary component vs. real component.

Case 2

Case 2 was likewise run with the parameters of Example 6 (the time record is 3 periods of the lowest frequency). Results are shown in FIGS. 28A-28C in the same format as Case 1.

Case 3

Again, Case 2 for rectification with zeros is virtually the same as non-zero rectification. Corruption by LPM transient response is suspected. Case 3 was run exactly as Case 3 for non-zero rectification (a time record of 2 periods of lowestfrequency and process with only the second period via zero-rectification FST). The results with the FST algorithm that uses rectification functions with zeros is virtually identical to the non-zero rectification-based FST in that eliminating the firsthalf of the 2-period time record totally eliminates the corruption of results by the LPM transient response. Results are shown in FIGS. 29A-29C in the same format as Case 1.

Conclusions

Zero and non-zero rectification for the subject method using the FST system have both been validated analytically. The rectification without zeros using the preferred means of non-averaging is probably easier to implement for FST as it involvespreprocessing the time record by just sign changes, and then just summing all the samples.

It is to be understood that the foregoing examples are merely illustrative of the present invention. Certain modifications of the articles and/or methods employed may be made and still achieve the objectives of the invention. Suchmodifications are contemplated as within the scope of the claimed invention.