There are no related applications.
The methods and systems relate to identifying system parameters, and more particularly to methods and systems for identifying parameters in Finite Impulse Response (FIR) systems.
Mathematical models of systems can be used in a wide variety of disciplines. These mathematical models can be used to predict how a system will react in different environments and can be used to improve the system. One approach to mathematical modeling may be known as “white box” modeling. In this approach, the laws governing a system are represented mathematically. For example, for a physical system, the law of gravity may be represented in the model. The outputs of the system may be determined by applying the mathematical model to the inputs. Sometimes, however, system behavior is too complex to be easily modeled using the “white box” approach.
Another approach to mathematical modeling may be known as system identification (Sysld). System identification may include mathematical modeling of systems from experimental data and may be applied in many engineering areas, e.g., communications, mechanical engineering, or geophysical engineering for purposes such as spectral analysis, fault detection, adaptive filtering, linear prediction, and pattern recognition. Sysld techniques may also find application in fields such as biology, environmental sciences and econometrics. In Sysld, a system can be treated as a black box, i.e., inputs can be applied to an unknown system, and a corresponding (noise corrupted) output signal from that system can be observed. These input and output signals may then be used, together with a mathematical model of the system, to estimate the model parameters that best characterize the system. When the system has M input channels, traditionally one may repeat the experiment M times with a different input channel selected to be active for each repetition. Additionally, traditional methods may involve inputting one or more sets of input probe signals and observing the respective output signals. Auto-correlations among the probe signals and cross correlations between the input and output signals may then be computed and used to solve a linear system of equations to determine the system coefficients, or parameters.
The need to solve a linear system of equations can be expensive, in terms of computations. As such, there is a need for computationally efficient methods and systems that eliminate the need to solve a linear system of equations.
Accordingly, the present invention is directed to methods, systems and/or computer programs that obviate the above and other disadvantages of the prior art. In general, methods, systems and/or computer programs disposed on a computer readable medium can determine characteristics of a Finite Impulse Response (FIR) system having M≧1 inputs and a response duration T≧1, by: determining a snapshot having a duration S at least equal to a multiple D of the number of inputs, the multiple D being at least equal to the response duration of the system; determining a number of identical sets of probe signals having probe components taken from a sequence of length S; applying at least one of the sets to the system to obtain output signals; obtaining transformed companion components related to the probe components; obtaining a correlation estimate based on components of the output signals and the transformed companion components; and selecting the characteristics from the correlation estimate.
In one embodiment, the output signals can be averaged over the number of sets of probe signals applied so as to obtain averaged output signals and the correlation estimates can be obtained based on the averaged output signals. In another embodiment, the transformed companion components can be obtained as an inverse of the components of a discrete Fourier transform (DFT) of the sequence of length S. The correlation estimate can be obtained by multiplying corresponding components of a DFT of components of the output signals and the transformed companion components to obtain product components and computing the inverse DFT of the product components.
A preamble sequence of a size equal to the multiple D can be attached to a first of the sets of probe signals or can be attached to each set of probe signals. The preamble sequence can be defined such that the probe components of the preamble correspond to the D last components of one of the sets of probe signals. In one embodiment, a designated probe component m of an nth sample of a designated set of probe signals can be set equal to the (n−mD)th component of the sequence of length S, with (n−mD) being evaluated modulo S.
In yet another embodiment, the sequence of length S can be determined by selecting a partial set of transform components having arbitrary, non-zero, real number first and last partial set transform components and arbitrary, non-zero, complex number components for other than the first and last partial set transform components, adding the non-zero, complex number components to the partial set in reverse order and in complex conjugate to obtain a full set of S transform components and obtaining an inverse DFT of the full set of transform components. The transformed companion components can be obtained by taking inverses of the full set of transform components.
The correlation estimate can be obtained by obtaining a DFT of components of the output signals, multiplying corresponding components of the DFT of components of the output signals and the transformed companion components to obtain product components and computing the inverse DFT of the product components. The characteristic can be selected by setting a characteristic at a designated tap q for a designated system output k from a designated system input p to a correlation estimate at correlation lag (pD+q) corresponding to the kth component of the output response to the applied probe signals. The DFT's and inverse DFT's can be obtained by the use of a Fast Fourier Transform algorithm and an inverse Fast Fourier Transform algorithm, respectively.
In further embodiments, methods and computer programs for determining characteristics of a FIR system can include applying identical sets of probe signals to the FIR system to obtain output signals, where the probe signals can have probe components taken from a circulant matrix of a size at least equal to a result of a multiplication of a number of FIR system inputs and a response duration of the FIR system. The method can obtain product components based on multiplication of corresponding components of a DFT of components of the output signals and a DFT of companion components derived from a circulant inverse of the circulant matrix and the characteristics can be selected based on computing an inverse DFT of the product components.
In other embodiments, methods and computer programs can determine system characteristics of a FIR system by: selecting a partial set of transform components having arbitrary, non-zero, real number first and last partial set transform components and arbitrary, non-zero, complex number components for other than the first and last partial set transform components; adding the non-zero, complex number components to the partial set in reverse order and in complex conjugate to obtain a full set of transform components, the full set having a number of transform components at least equal to a multiple D of the number of inputs, the multiple D being at least equal to a response duration of the system. The method can obtain an inverse DFT of the full set of transform components to obtain a set of probe components and can apply a number of identical ones of the set of probe components to the FIR system to obtain output signals. The method can further obtain transformed companion components from inverses of the full set of transform components, obtain product components based on multiplication of corresponding components of a DFT of components of the output signals and the transformed companion components and select the characteristics based on computing an inverse DFT of the product components.
In still another embodiment, systems for determining characteristics of a FIR system can include means for inputting identical sets of probe signals to the FIR system, means for measuring output signals in response to the probe signals, means for determining transformed companion components related to components of the probe signals, means for obtaining a DFT of components of the output signals, means for multiplying corresponding components of the DFT of components of the output signals and the transformed companion components to obtain product components and means for selecting characteristics of the FIR system based on an inverse DFT of the product components.
In a further embodiment, a system can include means for selecting a partial set of transform components having arbitrary, non-zero, real number first and last partial set transform components and arbitrary, non-zero, complex number components for other than the first and last partial set transform components and means for adding the non-zero, complex number components to the partial set in reverse order and in complex conjugate to obtain a full set of transform components, the full set having a number of transform components at least equal to a multiple D of the number of inputs, the multiple D being at least equal to a response duration of the system. The system can include means for obtaining a set of probe components by taking an inverse DFT of the full set of transform components, means for obtaining output signals by applying a number of identical ones of the set of probe components to the FIR system and means for obtaining transformed companion components from inverses of the full set of transform components. The system can further include means for obtaining product components based on multiplication of corresponding components of a DFT of components of the output signals and the transformed companion components and means for selecting the characteristics based on computing an inverse DFT of the product components.
The following figures depict certain illustrative embodiments in which like reference numerals refer to like elements. These depicted embodiments are to be understood as illustrative and not as limiting in any way.
By analyzing response 110 of system 112 to a given set of inputs 100, a relationship can be established between a given input 100 and a given output 110. One method to determine such relationships may be to apply a zero input signal 100 to all but one of the actuators 104, i.e., apply only one scalar input signal 100 to one of the actuators 104, and measure the responses 110 of the sensors 106 to the single non-zero input signal 100. This process can be repeated until the different input signals 100 have been applied. An approach, or method can be described herein that applies non-zero input signals 100 to multiple input channels 104 simultaneously. Such input signals 100 can be considered vector probe signals 100. The method can provide a low-cost method of Sysld, in terms of computational requirements.
In greater detail, response 110 of system 112 to the inputs 100 can be modeled as the noise-corrupted output of a Finite Impulse Response (FIR) system. The input-output relationship of such a system can be compactly represented according to:
where Bl can represent the matrix of unknown system coefficients for the lth tap; en can represent the nth time sample of a vector-valued zero-mean noise process; and the input-output signals un and Vn, respectively, can be observed in K≧1 equal-duration snapshots Wk≡{[un,vn(k)];n=0,1, . . . ,S−1 for k=0,1, . . . ,K−1. The input signal, un can be the same for each snapshot. Moreover, for some D≧T, the input signal may be defined and applied for the entire time indicated by the index range −D≦n≦S−1. The output signals may be averaged over the snapshots
and Bl can be estimated using [un, {overscore (v)}n] for n=0,1, . . . , S−1. It can be understood that the averaging of the output signals may be used to improve the signal to noise ratio (SNR) of the output signal and the resulting fidelity of the estimates. However, if one snapshot can be applied at a sufficiently high level, the fidelity of the resulting estimate may be adequate, so that no additional snapshots may need to be applied and averaged, i.e., K=1.
The choice of input signals un may include the use of a circulant matrix of probe components
and the companion components derived from circulant inverse C−1={overscore (C)}, where
the element in the mth row and nth column of C and {overscore (C)} can be expressed as Cm,n=cn⊖m and {overscore (C)}m,n={overscore (c)}n⊖m, respectively, where n⊖m≡n−m mod S for n,m=0,1, . . . ,S−1. Thus,
where δm,τ=1 for m=τ and δmτ=0 for m≠τ.
If un can be taken to have M components and S≧MD, then, for n=0,1, . . . ,S−1 and m=0,1, . . . ,M−1, let
un;m≡cn⊖mD, [6]
where un;m can denote the mth component of un. If {overscore (v)}n may have N components, the input-output relationship of Equation [1] can be rewritten as:
for n=0,1, . . . ,S−1 and k=0,1, . . . ,N−1, where Bl;k,m may be the element in the kth row and mth column of Bl, while {overscore (v)}n;k can denote the kth component of {overscore (v)}n. Additionally, the use of a preamble sequence, which may be referred to as cyclically padding the probe signal, can place the unknown system in a known state at the beginning of each snapshot. The preamble sequence can be defined for n=−D,−D+1, . . . ,−1 by letting un=uS+n for n=−D,−D+1, . . . ,−1. For the exemplary embodiment described herein, the preamble sequence may cyclically pad the first of the sets of probe signals. However, it can be noted that the preamble sequence may be used to cyclically pad each of the sets of probe signals. Using the cyclically-padded circulant probe signal in Equation [7] can yield
since un−l;m=c(n−l)⊖mD=cn⊖(mD+l).
Correlations can be considered such that, for τ=0,1, . . . , S−1 and using Equation [5],
If τ=pD+q for some pε{0,1, . . . ,M−1} and some qε{0,1, . . . ,T−1}, then Equation [9] can become
{circumflex over (R)}(k, pD+q)=ηpD+q,k+Bq;k,p, [11]
otherwise {circumflex over (R)}(k,τ)={overscore (n)}τ,k. Thus, it can be seen to be reasonable to use the parameter estimates
{circumflex over (B)}q;k,p≡{circumflex over (R)}(k,pD+q) [12]
for pε{0,1, . . . ,M−1}, qε{0,1, . . . ,T−1}and kε{0,1, . . . ,N−1}.
It can be understood from the above that one can obtain estimates of the system parameters by computing the circular correlation {circumflex over (R)}(k,τ), as defined in Equation [9]. As previously noted, traditional methods of determining the system parameters can involve solving a linear system of equations using auto-correlations of the input signals and cross correlations between the input and output signals. For the method as described above, the parameter estimates can be obtained without the need to explicitly solve a system of equations.
The circular correlation {circumflex over (R)}(k,τ) can be computed using one of a multiple of techniques that may be known in the art. In a preferred embodiment, one can compute the parameter estimates {circumflex over (B)}q;k,p, by considering the following transforms:
and
for m=0,1, . . . ,S−1. Referring to Equation [9], it can be verified that
for τ=0,1, . . . ,S−1. Thus, for each kε{0,1, . . . , N−1}, {overscore ({tilde over (v)}m;k may be computed in accordance with Equation [13] using the known Fast Fourier Transform (FFT) algorithm, and {overscore ({tilde over (c)}m{overscore ({tilde over (v)}m;k can be computed for each m=0,1, . . . , S−1. Then, {circumflex over (R)}(k,τ) can be computed in accordance with Equation [15] using the inverse FFT algorithm. The estimates can then be selected from these results, in accordance with Equation [12].
The input signals un can then be applied to the system for the K equal duration snapshots at 206, and the output signals vn can be observed at 208. It can be noted that the input signal un can be the same for each snapshot and, for some D≧T, may be defined and applied for the entire time indicated by the index range −D≦n≦S−1. Using Equation [2], the output signals can be averaged, at 210, over the number of snapshots K to obtain {overscore (v)}n;k, where {overscore (v)}n;k can denote the kth component of averaged output signals {overscore (v)}n. It can be understood, as previously noted, that only a single snapshot (K=1) may be required. Thus, 210 may be shown in phantom in
Using Equation [13], one can obtain {overscore ({tilde over (v)}m;k for m=0,1, . . . ,S−1, i.e., the DFT of {overscore (v)}n;k, by using the FFT algorithm at 212, for each kε{0,1, . . . ,N−1}. The inverse of the circulant matrix C, i.e., C−1, can be obtained at 214, with the companion component sequence {overscore (c)}0,{overscore (c)}1, . . . ,{overscore (c)}S−1 taken from C−1. Using Equation [14], and with {overscore (c)}n denoting the companion component sequence {overscore (c)}0,{overscore (c)}1, . . . ,{overscore (c)}S−1, one may obtain {overscore ({tilde over (c)}n for m=0,1, . . . , S−1, i.e., the DFT of {overscore (c)}n, again by using the FFT algorithm at 216. Components {overscore ({tilde over (v)}m;k of the transform of {overscore (v)}n;k can be multiplied at 218 with corresponding components {overscore ({tilde over (c)}m of the transform of {overscore (c)}n to obtain product components {overscore ({tilde over (c)}m{overscore ({tilde over (v)}m;k and Equation [15] can be used to obtain, at 220, the inverse DFT of these product components, also using the FFT algorithm. The parameter estimates, {circumflex over (B)}q;k,p of Equation [7], can then be selected at 222 from the result of 220 using Equation [12].
Alternatively, those of skill in the art may recognize that the correlation estimate {circumflex over (R)}(k,τ) can be obtained directly from Equation [9] using known techniques other than the discrete Fourier transforms described herein. For example, if one chooses not to use the discrete Fourier transforms, as determined at 224, the companion components can be determined at 226 as described previously with relation to obtaining the inverse circulant matrix C−1. Using the averaged output signals {overscore (v)}n;k obtained at 210 and the companion components {overscore (c)}n obtained at 226, the correlation estimate can be obtained from Equation [9], as shown in 228. The parameter estimates can be selected (222) based on the correlation estimate so obtained.
It may also be understood that in defining the input signals at 202, the circulant matrix C need not be explicitly formed and that the inverse circulant matrix C−1 need not be explicitly computed at 214. In practice, one can first select transform components {tilde over (c)}1,{tilde over (c)}2, . . . ,{tilde over (c)}(S/2)−1 as arbitrary (non-zero) complex numbers and transform components {tilde over (c)}0,{tilde over (c)}S/2 as arbitrary (non-zero) real numbers (202a). Then one can define the transform components {tilde over (c)}n={tilde over (ċ)}S−n, where {tilde over (ċ)}S−n denotes the complex conjugate of {tilde over (c)}S−n for n=1+S/2,2+S/2, . . . ,S−1 (202b), so that the inverse DFT of {tilde over (c)}0,{tilde over (c)}1, . . . ,{tilde over (c)}S−1 (202c) can be the real sequence of probe components c0,c1, . . . ,cS−1. Furthermore, the transformed companion sequence obtained at 216 can be obtained directly from {overscore ({tilde over (c)}0,{overscore ({tilde over (c)}1, . . . ,{overscore ({tilde over (c)}S−1=1/{tilde over (c)}0,1/{tilde over (c)}1, . . . , 1{tilde over (c)}S−1, without explicitly computing the companion sequence {overscore (c)}0,{overscore (c)}1, . . . ,{overscore (c)}S−1, though, if desired, the companion sequence {overscore (c)}0,{overscore (c)}1, . . . ,{overscore (c)}S−1 (obtained at 214) can be computed as the inverse DFT of the transformed companion sequence {overscore ({tilde over (c)}0,{overscore ({tilde over (c)}1, . . . ,{overscore ({tilde over (c)}S−1.
The process 200 of
The memory 312 may also include instructions 310. The instructions 310 can be transferred, in the course of operation, from the non-volatile memory 312 to the volatile memory 308 and processor 306 for execution. The instructions 310 may output system input values, as indicated at 302, and receive system output values, as indicated at 304, during system testing. The device 300 may communicate with a user via a monitor or other input/output device 316 such as a keyboard, mouse, microphone, and so forth. Additionally, the device 300 may feature a network connection 318.
The techniques described herein may not be limited to a particular hardware or software configuration; they may find applicability in many computing or processing environments. The techniques may be implemented in hardware or software, or a combination thereof. Preferably, the techniques are implemented in computer programs.
The computer programs may be implemented in high level procedural language (e.g., MATLAB®, The MathWorks, Inc., Natick, Mass.) or an object oriented programming language to communicate with a computer system. However, the programs can be implemented in assembly or machine language, if desired. In any case the language may be compiled or interpreted language.
The computer programs can be stored on a storage media or device(s) (e.g., CD-ROM, hard disk, or magnetic disk) that are readable by a general or special purpose programmable computer for configuring and operating the computer when the storage medium or device is read by the computer to perform the procedures described herein. The system may also be considered to be implemented as a computer-readable storage medium, configured with a computer program, where the storage medium so configured causes a computer to operate in a specific and predefined manner.
Additionally, while described in conjunction with an example of a shaken physical structure, the approaches described may be applied to a wide variety of systems such as biological, signal processing, and so forth. Additionally, the spectrum of the circulant probe signal sequence c0,c1, . . . ,cS−1 can be arbitrarily specified, subject only to the mild condition that the corresponding circulant matrix C must be non-singular, so that the probe energy can have any desired spectral distribution. This may be desirable if the observation noise may have a spectral distribution and/or certain spectral portions of the system 112 may be identified more accurately than others. Other embodiments are within the scope of the following claims.
The U.S. Government has a paid-up license in this invention and the right in limited circumstances to require the patent owner to license to others on reasonable terms as provided for by the terms of Contract No. N00024-91-C-4355 awarded by the Naval Sea Systems Command.
Number | Name | Date | Kind |
---|---|---|---|
5347586 | Hill et al. | Sep 1994 | A |
5796849 | Coleman et al. | Aug 1998 | A |
5991348 | Alelyunas et al. | Nov 1999 | A |
6701335 | Pupalaikis | Mar 2004 | B2 |
20020126803 | Jones et al. | Sep 2002 | A1 |