Real-time digital signal recovery for a low-pass transfer function system with multiple complex poles

In order to solve the problems of waveform distortion and signal delay by many physical and electrical systems with linear low-pass transfer characteristics with multiple complex poles, a general digital-signal-processing (DSP)-based method of real-time recovery of the original source waveform from the distorted output waveform is proposed. From the convolution kernel representation of a multiple-pole low-pass transfer function with an arbitrary denominator polynomial with real valued coefficients, it is shown that the source waveform can be accurately recovered in real time using a particular moving average algorithm with real-valued DSP computations only, even though some or all of the poles are complex. The proposed digital signal recovery method is DC-accurate and unaffected by initial conditions, transient signals, and resonant amplitude enhancement. The noise characteristics of the data recovery shows inverse of the low-pass filter characteristics. This method can be applied to most sensors and amplifiers operating close to their frequency response limits or around their resonance frequencies to accurately deconvolute the multiple-pole characteristics and to improve the overall performances of data acquisition systems and digital feedback control systems.


Introduction
Many sensors and amplifiers suffer from delayed and distorted responses with single-or multi-pole low-pass characteristics when operated close to their frequency response limits or resonance frequencies. Notable examples of multi-pole systems include a spring-mass system driven by external force or displacement ( Fig. 1(a)), or an L-R-C-based circuit driven by external voltage or current ( Fig. 1(b)). Unlike a first-order-response system, a second or higher order system may have complex-conjugate-paired poles even though all the coefficients of the denominator polynomial are real. This condition may induce resonant underdamped transfer characteristics and requires a carefully thought-out digital signal recovery algorithm in a conventional DSP system. Here I propose a real-time numerical waveform recovery method that is suitable for the case of arbitrary order denominator polynomials with multiple complex poles and discuss their detailed implementation method, accuracy and noise characteristics.

Generalization to high-order real systems with complex poles
Let's consider the two-pole low-pass system as shown in Fig. 1(a) driven by an external displacement to one end of the spring while the output of this system is the displacement of the mass. In case of small enough damping parameter b, the two poles can be complex and a complex-conjugate pair. This condition makes the previous method with real-valued poles [1] becomes inappropriate due to the requirement of complex-valued calculations (i.e. = − are complex for complex ) in the real-time core of the DSP/FPGA signal processor.
Here I will introduce an appropriate digital signal recovery method for a general two-pole linear low-pass system by direct mathematical expansion starting from the convolution kernel expression of the two pole transfer characteristics and show that it is equivalent to the nested multi-pole solution shown in Fig. 4 of Ref. [1]. The results will then be generalized into higher order low-pass system whose Laplace transfer function has an n-th order real-coefficient polynomial as its denominator.
The Laplace representation of the first order low pass characteristics with an initial condition (Ref. [2]) with time constant 1 = 1/ 1 can be converted to the time domain expression of ( ) = ∫ We can recover the original signal as derived in Ref.
since the contribution from the initial condition term ( ) = ( = 0) − 1 Θ( ) vanishes as for arbitrary > . Now let's consider the general second order low-pass characteristics with initial conditions This can be converted to the time domain solution of ( ) = 0 ( ) + ( ) which is the sum of the nested time domain convolution and the transient solution dependent on the initial conditions Since which is graphically illustrated in Fig. 2(d)-(h).
In the limit of small , the average value of over the square region of integration with parameter range [ − 2 , ] can be approximated by ( − ) and the left-hand side expression is simplified as (16) Therefore we have Since for both = 1, 2 we have the contribution from the initial-condition-dependent transient solution ( ) vanishes for arbitrary > 2 when inserted in the places of the 0 ( ): Therefore the signal recovery formula for ( ) maintains the same form as the Eq. (17): We can also show that the above Eq. (20) can be converted to a nested form The above result can be generalized to an arbitrary high order . The general -th order low-pass characteristics with initial conditions is given by where is a linear combination of the initial conditions (0,1,2,…, −1) ( = 0) [2].
The time-domain solution is given in the form and The signal recovery formula (corresponding to Eqs. (12) and (20)) is generalized to (29) where the coefficients for ( − ) in the numerator are simply the coefficients for the − term in the generating polynomial (30) We can again show that the transient solution ( ) gives no contribution to the right-hand side of the Eq. (29) since, for the arbitrary -th term of the Eq. (28) proportional to − , we have (32) and the right-hand side of the Eq. (29) vanishes when ( ) is replaced by ( ).
This provides a general proof that the signal recovery method shown in the Eq. (29) produces output signal completely independent of the initial conditions and the transient signals. Fig. 3 shows three equivalent representations for physical implementation of high order signal recovery in the DSP/FPGA device. They are mathematically equivalent but in case of complex poles, the first implementation ( Fig. 3(a)) requires complex computation while the second ( Fig. 3(b)) and the third (Fig. 3(c)) require only real-valued computation which has a significant advantage in the real-time process core of the DSP/FPGA. We note that the complex roots of any real-coefficient polynomial always occur in complex-conjugated pairs. Therefore in the device implementation of Fig. 3(b) (formed by combining every pair of register loops of Fig. 3a for a complex conjugate pair (such as 1 and 2 with = 2 − 4 < 0) into a single register loop), it is sufficient to show that only real-valued computations are needed for the evaluation of the combined second order Eq. (20).

Device Implementation
Let's assume that 1 = + and 2 = − where and are real and > 0. Then we have the following values all real: and all the numbers used in the evaluation of the Eq. (20) become real The device implementation of Fig. 3(c) (formed by combining all the register loops of Fig. 3a into one) also has the property of real-value-only computations. The proof that all the coefficients in the numerator of the Eq. (29) and its denominator are realvalued is easily provided noting that if 1 and 2 are a complex-conjugate pair, − 1 and − 2 are also. As a result, all the coefficients of the polynomial expansion of the Eq. (30) are real. Therefore, all the coefficients appearing in the numerator of Eq. (29) are real and the denominator (1) is real.

Noise consideration
In order to understand the noise characteristics of the multi-pole recovery method quantitatively, let's assume without too much loss of generality, that the high-order-convoluted analog output signal � ( ) contains a slow-varying (compared to 2 ) raw signal ( ) plus a pseudo-random noise ( ) whose correlation time is shorter than . First, let's look at the nontrivial second-order case of Eq. (35) where the two roots 1 and 2 are a complex-conjugated pair. Then the numerical recovery operation applied to � ( ) = ( ) + ( ) can be divided into two terms This is formally identical to the result when the two poles are real Generalization to an arbitrarily high order case of Eq. (29) leads to As mentioned in Ref. [1], if we sample � ( ) (≥ 2) times over the short time intervals within [ , − ) and use their averages in place of � ( ), we may further increase the S/N by up to a factor given by a fraction of � . The factor can approach � in case the correlation time of the noise is still shorter than the sampling periods of the data points. On the other hand, when it is possible to perform multiple measurements over a repeated input signal, we can increase the S/N to an arbitrary level by choosing the averaging by

Simulated demonstration
I performed simulations for a two-pole spring-mass-damper system shown in Fig. 1(a) whose results are shown in Fig. 4-6 and for a three-pole passive electrical circuit system (shown in Fig. 1(b)) whose results are shown in Fig. 7-9. In each case, we first numerically solved the (second or third order) differential equations for three different kinds of input waveforms and processed it with the real-time data recovery scheme shown in Fig. 3(c) with optimal and non-optimal choices of parameters. The signal recovery scheme of Fig. 3(b) produced numerically identical results as the scheme of Fig. 3(c) and therefore not separately reproduced here.
As can be seen in all Figs. 4-9, the recovered waveform closely matches with the input waveform only when the parameter used in evaluating all the coefficients of the Eq. (20) or (29) matches with the actual time difference between samplings, leading to the optimal values of all the coefficients and hence the optimal overall compensation. Also it should be noted that the signal recovery is DC-accurate and independent of the initial conditions, the transient waveforms and the waveforms significantly amplified near the resonance frequency.

Conclusion
A relatively simple digital-signal-processing-based method of real-time signal recovery is proposed, which can compensate for the waveform distortion and propagation delay due to single-pole or multiple-complex-pole low-pass transfer characteristics in many physical and electronic systems. In case the transfer function has a real-coefficient polynomial as its denominator, we can use signal processing based on real-valued computations only, even though some of the poles are complex. The overall method is also initial-value-independent and will be especially useful in exactly deconvoluting the multi-pole transfer characteristics, in improving the performances of data acquisition systems and in stabilizing high speed feedback control systems with sensors and amplifiers operated close to their frequency response limits or around their resonance frequencies, utilizing modern low-cost high-speed DSPs and FPGAs [3][4][5][6][7][8][9].   (a) In case all the are real and positive, the nested multiple register scheme as illustrated in Ref. [1] can be used. However, in case some of the are complex, the scheme is not efficiently realized in real-valued digital signal processors. Two alternative solutions are suggested: (b) Combining two registers for every complex conjugated pair of the (e.g. = 2 − 4 < 0) into one, while leaving the registers for real left nested as before, all the digital signal processing is done with real parameters and real-valued digital signal processing only. (c) Combining all the registers of (a) into one register also results in the digital signal processing done with real parameters and real-valued digital signal processing only.     Fig. 3(c) for a system with three-pole transfer characteristics shown in Fig. 1(b), tested with a frequency-varying rectangular input waveform. (a) Input waveform ( ). (b) Response of the system output ( ) calculated by numerically solving the differential equation with the two initial conditions (0) and ̇(0) arbitrarily chosen. (c) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters intentionally offset from their optimal values by replacing every by 0.9 . (d) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters set at their optimal values. (e) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters intentionally offset from their optimal values by replacing every by 1.1 . Note that the strong transient waveforms after the sharp transitions and the strong resonant waveform near 4000 < < 7000 are exactly cancelled out in (d).  Fig. 3(c) for a system with three-pole transfer characteristics shown in Fig. 1(b), tested with a frequency-varying cosine-cubed input waveform. (a) Input waveform ( ). (b) Response of the system output ( ) calculated by numerically solving the differential equation with the two initial conditions (0) and ̇(0) arbitrarily chosen. (c) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters intentionally offset from their optimal values by replacing every by 0.9 . (d) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters set at their optimal values. (e) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters intentionally offset from their optimal values by replacing every by 1.1 . Note that the strong transient waveform near 0 < < 800 in (b) and the strong resonant waveform near 3500 < < 5500 are exactly cancelled out in (d). Fig. 9. Numerical simulation of the real-time waveform recovery in the scheme of Fig. 3(c) for a system with three-pole transfer characteristics shown in Fig. 1(b), tested with an aperiodic input waveform. (a) Input waveform ( ). (b) Response of the system output ( ) calculated by numerically solving the differential equation with the two initial conditions (0) and ̇(0) arbitrarily chosen. (c) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters intentionally offset from their optimal values by replacing every by 0.9 . (d) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters set at their optimal values. (e) Compensated waveform ( ) (overlaid with the input waveform ( )) calculated from the waveform in (b) with all the parameters intentionally offset from their optimal values by replacing every by 1.1 . Note that the strong transient waveforms after the sharp transitions are exactly cancelled out in (d) and the recovery is DC-accurate even with the parameters slightly offset from their optimal values as shown in (c)-(e).