The present invention relates to signal processing methods, apparatus and systems. The invention is useful in a diverse range of fields including, but not limited to, electrophysiology, non-destructive testing (NDT) and geoscience. More particularly, the invention relates to a signal processing technique which is analogous to the operation of a physical gradiometer device, but is of wider application and does not suffer the limitations imposed by the physical structure of such devices. The technique relates especially to the extraction of signals and other information from data subject to noise and/or background drift.

Experimental data in fields as diverse as electrophysiology [1], NDT [2], geoscience etc. often comprise one or more short signal pulses superimposed on a slow moving background drift. In many of these cases the background drift is generated by a quite different source to that which produces the signal pulses, and so it is desirable to be able to separate these two distinct signal components. A variety of signal processing techniques - using, for example, Fourier [3] or Walsh [4] Transforms, or, more recently, wavelet theory [5] - have addressed this problem. However, for most such data sets it is necessary, if not crucial, that the chosen technique should preserve the position, amplitude and shape of the signal pulses, and it is well known that existing methods can fail to preserve all three of these characteristics [3,6].

By way of background to the present invention, the operation of physical gradiometer devices will now be discussed. However, it will be understood that the invention itself does not relate to physical devices of this type, but to more general data processing techniques which operate in a manner analogous to such devices.

Gradiometers are devices which detect rates of change of physical parameters, most commonly of electrical, magnetic or gravitational fields. Typically, they are required to discriminate between field components deriving from large, distant sources and smaller local sources.

For example, magnetic field gradiometers coupled to Superconducting Quantum Interference Devices (SQUIDs) are widely used for the detection of small fields from local sources in the presence of much larger fields from more distant sources. This situation arises, for example, in NDT and in biomagnetism (the study of magnetic fields from biological sources such as electrical activity in the brain). In the case of biomagnetism, field magnitudes are of the order of 0.2 pT for stimulus-invoked brain activity. Ambient noise levels vary widely depending on location and whether magnetic shielding is used, but white noise of 1 nT Hz^-1/2 with discrete frequency interference of up to 50 µT is typical. A bandwidth of at least 100 Hz is required for most medically useful measurements.

The field from magnetic sources generally decays with distance r from the source as r^-n, where 1 < n < 3, depending on the nature of the source. Typical signal sources in NDT are magnetic dipoles (n = 3) and those in biomagnetism are current dipoles (n = 2). The sources are normally external to the sensor structure. Discrimination against distant sources can be achieved by a device geometry responsive primarily to the gradient or higher spatial derivatives of the field. In practice, this is approximated by a form of spatial filtering. An N-th order gradiometer is an arrangement of inductive pick-up loops designed so that the first N terms of a Taylor expansion of the field threading the gradiometer couple zero net magnetic flux into it.

Early gradiometers, as illustrated in Fig. 8(a), were fabricated by wire wound on a machined transformer, with complex mechanisms to compensate for geometrical imperfections. In the drawing, "S" indicates a typical location of the source of the field of interest. Such devices were designed to sense field components in the x (or radial) direction, for two reasons. Firstly, this component of the field has a simpler intuitive interpretation and, secondly, early simplistic models of biomagnetic sources indicated no y or z (tangential) components. In fact, more realistic source models and measurements using tangential field sensors show the tangential and radial field components to be comparable in magnitude. In the case of NDT all field components are significant.

More recent planar gradiometers sense tangential z components of the field and its differentials along the x axis (see Fig. 8(b) and Fig. 9). They are fabricated with precise geometrical definition using photolithographic thin-film patterning in a plane, obviating the need for mechanical trimming. A simple gradiometer of this type is illustrated in Fig. 9. It consists of a series of adjacent, rectangular pick-up loops, made of super-conducting thin films, wired in series with an insulated cross-over between each pair of loops. In the example shown, insulated cross-overs exist at x₁ and x₂; thus, the only continuous current carrying path is ABGFCDEFCBGHA. The structure is fabricated with wires or thin-film metal tracks, and any electric current induced by magnetic flux is monitored. The dimensions are such that a desired signal source will normally lie in the interval 1 < x < 2 (i.e. the source distance from the gradiometer should be comparable to the size of the gradiometer), and the source sensitivity is greatest in the plane of the gradiometer.

The positions and number of the cross-overs determine not only the gradiometer order N, but also the detailed response to nearby sources. In particular, the number of cross-overs m, where m ≥ N, should be minimised, subject to further considerations, to minimise the gradiometer self-inductance. This, due to details of SQUID operation, maximises sensitivity. The cross-over positions are selected so as to achieve a given gradiometric order N in the case where one or more cross-over positions have been pre-selected to achieve a given near-source response. This is required to give two specific properties. Firstly, it is highly desirable to locate the SQUID near the gradiometer on the same substrate, typically in the interval -2 < x < -1 (Fig. 9). However, the SQUID itself, as a superconductor, constitutes a nearby magnetic anomaly; i.e. it causes local distortion of the quasi-uniform interference fields and thus behaves as a local noise source. If the position of the cross-over nearest the SQUID (i.e. x₂ in Fig. 9) can be adjusted, a null-response to this source can be designed in. Secondly, the position of the cross-over nearest the signal source (x₁ in Fig. 9) can be adjusted to give a peak in gradiometer sensitivity as a function of distance at x > 1. This can be used either for depth-profiling (by moving the gradiometer along the x-axis) or for discrimination against an overlying source (e.g. detection of a nerve underlying surface muscle).

Gradiometers of this general type have the effect of filtering out low-order spatial components of the parameter being monitored, such as constant or uniformly varying background fields which would otherwise obscure a more rapidly-varying, local signal source. Theoretically, the operation of such devices would be improved by the inclusion of more "stages", corresponding to higher-order components of the parameter being monitored. However, there are practical difficulties in manufacturing such devices. Also, particular devices have to be designed for and limited to quite specific functions.

The present invention relates to a signal processing technique which, in principle, is somewhat similar to the principle of operation of physical gradiometers. For this reason, the inventors have coined the term "Gradiometer Transform" to describe the basic method which underlies the invention. The technique in accordance with the invention operates on sampled, numerical data, so that it can be extended to arbitrarily high orders, and may be applied to any type of data.

The invention embraces: a general data or signal processing method, which may be applied to stored data or to real-time or "live" signals and which may be implemented in software, firmware or hardware; apparatus incorporating means for applying the method to signals or data generated by any kind of transducer or other signal or data source; and apparatus adapted to apply the method for specific purposes.

In accordance with a first aspect of the invention, there is provided a data processing method comprising generating a plurality of components of a gradiometer transform by correlating a normalised data-set with a plurality of gradiometer step functions of differing orders.

Said step functions are preferably selected to successively eliminate components of a Taylor expansion of said data-set.

Preferably, said data-set comprises sampled data and said correlation of the data-set with said step functions is performed using a modified data set in which each data sample is replaced by a piecewise continuous rectangular step function.

The components of the gradiometer transform can be used to obtain useful information regarding the data-set, such as information regarding the background drift of a time-varying signal.

The method preferably includes a subsequent step of calculating a plurality of coefficients of at least one polynomial representing said data-set. Preferably, coefficients of a plurality of polynomials of differing orders are calculated.

Preferably, said coefficients are calculated using a table of values specific to said step functions which can be used to operate on said transform components to obtain said coefficients. Said table of values is preferably calculated in advance and stored as a look-up table for use in the execution of the method.

Preferably also, said polynomials are fitted to said data-set and the fitted polynomials are subtracted from the data-set.

It is particularly preferred that the data-set comprises a signal and the method is employed to remove background drift from said signal.

In accordance with a second aspect of the invention there is provided data processing apparatus comprising data processing means for applying the method in accordance with the first aspect of the invention to a data-set from a data source.

Preferably, said apparatus includes output means for generating an output data-set resulting from the application of said method.

In certain embodiments, the apparatus includes input means for generating an input signal in response to variations in a parameter to be monitored by the apparatus. Preferably also, the apparatus includes sampling means for dividing said signal into discrete samples.

The data-processing means may comprise a suitably programmed computer or a dedicated unit in which the method is implemented by means of a program encoded in hardware or firmware.

The data-source may be a sensor for monitoring a physical parameter or a store of previously obtained or synthesised data.

The underlying theory, examples and embodiments of the invention will now be described, by way of example only, with reference to the accompanying Tables and drawings in which:

• Table I shows the first fifteen rows and columns of the generic matrix J^-1.
• Table II shows the locations of the maxima and minima of the function illustrated in Fig.4(a) compared with the locations of the maxima/minima of the gradiometer transform fits to the data given in Figs.4(b)-(d).

Figure 1 illustrates the first six gradiometer step functions. The step function h_m(t) has crossovers at ${\mathit{\text{t}}}_{\mathit{\text{s}}}\text{= cos(}\mathit{\text{s}}\text{π/(}\mathit{\text{m}}\text{+1))}$, s=1,2,...,m.

Figure 2 illustrates: (a) the background drift $\mathit{\text{d(t) =}}\text{5-3(}\mathit{\text{t}}\text{+1)}{\mathit{\text{e}}}^{\mathit{\text{-2t}}}$; (b) the pulse signal; and (c) the test function comprising the sum of the pulse signal in (b) and the drift in (a); and (d)-(f): the gradiometer transform components F_m for each of the functions in (a),(b),(c).

Figure 3 illustrates: (a) the test function f(t); (b)-(g) the function f(t) after using Gradiometer differencing to remove 0^th, 1^st, 2^nd, 3^rd, 4^th and 5^th order polynomials; and (h) the original signal pulses after subtraction of their mean value over [-1,1]. As Gradiometer differencing will always return functions whose mean value is zero, the mean has been removed from the pulses to aid comparison.

Figure 4. illustrates random noise (in the form a rand(rand-0.5) where rand is a random number between 0 and 1) added to the signal in (a); (b)-(d) show the generated data (a=0.5,1,2 respectively); and (e)-(g) show the reconstructed signal obtained by applying the Gradiometer Transform to the generated data.

Figure 5. illustrates: (a) a foetal magnetocardiogram in which magnetic field density B (arbitrary units) is plotted against time t (scaled so that -1≤t≤t1); (b) the gradiometer transform components F_m of the signal in (a); and (c) the same data after 30^th order gradiometer differencing.

Figure 6 illustrates: (a) a somatosensory evoked field (100 averages) in which magnetic field density is plotted against time t (scaled so that -1≤t≤t1); (b) the gradiometer transform components F_m of the signal in (a); and (c) the 50^th order fit to the data in (a) which is produced by the gradiometer transform.

Figure 7 illustrates: (a) a magnetic profile over a steel tube in which magnetic field density B is plotted against displacement s (scaled so that -1≤s≤1); (b) gradiometer transform components F_m of the signal in (a); and (c)&(d) the signal after gradiometer filtering (see text for details).

Figures 8(a) and 8(b) illustrate gradiometer geometries of a conventional wirewound axial field sensor and a planar tangential field sensor, respectively, by way of background information to the present invention.

Figure 9 shows an example of planar gradiometer geometry, by way of background information to the present invention.

Fig. 10 illustrates a typical implementation of the invention using a general purpose computer.

Fig. 11 illustrates a typical implementation of the invention in a dedicated hardware unit.

The recently developed technique of Gradiometer filtering [7] (discussed in greater detail above) has been used to successfully separate signals of the type described into two separate components, namely the signal pulses and the background drift, while maintaining the position, amplitude and shape of the signal waveform. Briefly, the Gradiometer Transform of the present invention works by firstly providing an estimate of the order of drift present and secondly enabling a polynomial function of arbitrary order to be rapidly and accurately fitted to the experimental data. As the background drift can in general be represented by a low order polynomial, while rapidly varying underlying signals cannot, a differencing or filtering method can be used to separate the two signal components. In Gradiometer differencing, a polynomial of sufficiently low order to represent the drift alone is fitted to, and then subtracted from, the signal; in Gradiometer filtering such a low order polynomial is subtracted from a high order polynomial fit which captures the complete signal - the advantage of filtering is that it will simultaneously eliminate noise. The present invention is based on a novel underlying theory of Gradiometer Transforms, which will facilitate rapid and accurate computation of both the transform components and the regenerated signal, together with an evaluation of its performance using both synthesised and collected data.

The key definitions of the Gradiometer Transform are as follow.

In the present example, an m^th order gradiometer step function takes values ±1 over the interval [-1,1], with crossovers at ${\mathit{\text{t}}}_{\text{s}}\text{= cos(}\mathit{\text{s}}\text{π/(}\mathit{\text{m}}\text{+1))}$, s = 1,2,...,m; the first eight of these step functions are shown in Fig.1. The component F_m (m = 0,1,2,...) of the Gradiometer Transform is found by correlating a given signal f(t), or other data-set, (scaled, or "normalised", so that - 1≤t≤1) with the m^th order gradiometer step function h_m(t) :$${\mathit{\text{F}}}_{\mathit{\text{m}}}\mathit{\text{=}}\int_{\text{-1}}^{\text{1}} \mathit{\text{f}}\text{(}\mathit{\text{t}}\text{)}{\mathit{\text{h}}}_{\mathit{\text{m}}}\text{(}\mathit{\text{t}}\text{)}\mathit{\text{dt}}\text{.}\,d$$

Evaluating N of these transform components generates a transform vector which can be used to obtain an (N-1)^th order polynomial fit to the signal, namely f(t) ≃ c₀ + c₁t + ... + c_N-1t^N-1 [9]. The vector of coefficients c = [c₀ c₁ ... c_N-1]^T is obtained from the transform vector F using the matrix equation where J^-1 is the inverse of an N×N matrix J = [J_m,n] whose elements are given by$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\int_{\text{-1}}^{\text{1}} {\mathit{\text{t}}}^{\mathit{\text{n}}\text{-1}}{\mathit{\text{h}}}_{\mathit{\text{m}}\text{-1}}\text{(}\mathit{\text{t}}\text{)}\mathit{\text{dt}}\text{.}\,d$$ (see [7]).

The Matrices J and J^-1 will now be described.

Here it is shown that the matrix J can be decomposed into the product of four matrices, three of which have well defined inverses and the fourth of which, being a sparse matrix whose entries are 0 and 1, has an inverse which can be easily calculated. Decomposing the matrix J in this way eliminates the error which can occur in using (4) to calculate the matrix elements numerically, and also the resulting error (naturally compounded) in calculating the inverse matrix J^-1. As in previous work with gradiometer step functions [8], it is convenient initially to consider separately those cases where the gradiometer step function has an even or an odd number of crossovers.

When m is even, the gradiometer step function h_m-1(t) will have an odd number of crossovers, located at $\mathit{\text{t}}\text{= cos(}\mathit{\text{s}}\text{π/}\mathit{\text{m}}\text{)}$, s = 1,2,...,m-1. It follows from the (4) above that the element in the m^th row and n^th column of the matrix J is given by$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{\text{1}}{\mathit{\text{n}}}\left[{\text{-(-1)}}^{\mathit{\text{n}}}\text{- 1 + 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\mathit{\text{m}}\text{/2}} {\text{cos}}^{\mathit{\text{n}}}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}}{\mathit{\text{m}}}\text{- 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\mathit{\text{m}}\text{/2-1}} {\text{cos}}^{\mathit{\text{n}}}\frac{\text{2}\mathit{\text{k}}\text{π}}{\mathit{\text{m}}}\right]$$

If we now expand the powers of $\text{cos}\mathit{\text{x}}$ in terms of cosines of multiple angles$${\text{cos}}^{\mathit{\text{n}}}\mathit{\text{x}}\text{=}\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}\text{cos}\mathit{\text{rx}}$$ and use the identities$$\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}\text{= 1 and}\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}{\text{(-1)}}^{\mathit{\text{r}}}{\text{= (-1)}}^{\mathit{\text{n}}}$$ (obtained by substituting x=0 and x=π respectively in the above), then

While the values of b_r,n, the coefficients of $\text{cos}\mathit{\text{rx}}$ in the expansion of ${\text{cos}}^{\mathit{\text{n}}}\mathit{\text{x}}$, are well documented [9], their precise values are not required here. Following Jones & Bain [8], we note that as m is even, the roots of the equation$${\mathit{\text{z}}}^{\mathit{\text{m}}}\text{+ 1 = 0}$$ are$$\mathit{\text{z}}\text{= cos}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}}{\mathit{\text{m}}}\text{±}\mathit{\text{i}}\text{sin}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}}{\mathit{\text{m}}}\text{(}\mathit{\text{k}}\text{= 1,2,...,}\frac{\mathit{\text{m}}}{\text{2}}\text{),}$$ and the roots of the equation$${\mathit{\text{z}}}^{\mathit{\text{m}}}\text{- 1 = 0}$$ are$$\mathit{\text{z}}\text{= ±1 and}\mathit{\text{z}}\text{= cos}\frac{\text{2}\mathit{\text{k}}\text{π}}{\mathit{\text{m}}}\text{±}\mathit{\text{i}}\text{sin}\frac{\text{2}\mathit{\text{k}}\text{π}}{\mathit{\text{m}}}\text{(}\mathit{\text{k}}\text{= 1,2,...,}\frac{\mathit{\text{m}}}{\text{2}}\text{-1),}$$ from which it follows that the sum of the r^th powers of the roots of (8), S_m,r, is$${\mathit{\text{S}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{= 2}\sum_{=\mathit{\text{r}}\text{=1}}^{\mathit{\text{m}}\text{/2}} \text{cos}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}\mathit{\text{r}}}{\mathit{\text{m}}}\text{,}$$ and the sum of the r^th powers of the roots of (9), T_m,r, is$${\mathit{\text{T}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}{\text{= 1 + (-1)}}^{\mathit{\text{r}}}\text{+ 2}\sum_{=\mathit{\text{r}}\text{=1}}^{\mathit{\text{m}}\text{/2-1}} \text{cos}\frac{\text{2}\mathit{\text{k}}\text{π}\mathit{\text{r}}}{\mathit{\text{m}}}\text{.}$$ It follows from Eq.(7) that J_m,n may be written as$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{\text{1}}{\mathit{\text{n}}}\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}\text{[}{\mathit{\text{S}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{-}{\mathit{\text{T}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{]}\text{.}$$

When m is odd, the gradiometer step function h_m-1(t) has an even number of crossovers, located at $\mathit{\text{t}}\text{= cos(}\mathit{\text{s}}\text{π/}\mathit{\text{m}}\text{)}$, s = 1,2,...,m-1, and the element J_m,n is given by$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{\text{1}}{\mathit{\text{n}}}\left[{\text{-(-1)}}^{\mathit{\text{n}}}\text{+ 1 - 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\text{(}\mathit{\text{m}}\text{-1)/2}} {\text{cos}}^{\mathit{\text{n}}}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}}{\mathit{\text{m}}}\text{+ 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\text{(}\mathit{\text{m}}\text{-1)/2}} {\text{cos}}^{\mathit{\text{n}}}\frac{\text{2}\mathit{\text{k}}\text{π}}{\mathit{\text{m}}}\right]$$ which, using the identities (5) and (6), may be written as$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{\text{1}}{\mathit{\text{n}}}\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}\left[\text{-2}\sum_{=\mathit{\text{k}}\text{=1}}^{\text{(}\mathit{\text{m}}\text{-1)/2}} \text{cos}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}\mathit{\text{r}}}{\mathit{\text{m}}}\text{+ 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\text{(}\mathit{\text{m}}\text{-1)/2}} \text{cos}\frac{\text{2}\mathit{\text{k}}\text{π}\mathit{\text{r}}}{\mathit{\text{m}}}{\text{+ 1 - (-1)}}^{\mathit{\text{r}}}\right]\text{.}$$

When m is odd, the m roots of Eq.(8) are$$\mathit{\text{z}}\text{= -1,}\mathit{\text{z}}\text{= cos}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}}{\mathit{\text{m}}}\text{±}\mathit{\text{i}}\text{sin}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}}{\mathit{\text{m}}}\text{(}\mathit{\text{k}}\text{= 1,2,...,}\frac{\mathit{\text{m}}\text{-1}}{\text{2}}\text{) ,}$$ and the m roots of Eq.(9) are$$\mathit{\text{z}}\text{= +1 and}\mathit{\text{z}}\text{= cos}\frac{\text{2}\mathit{\text{k}}\text{π}}{\mathit{\text{m}}}\text{±}\mathit{\text{i}}\text{sin}\frac{\text{2}\mathit{\text{k}}\text{π}}{\mathit{\text{m}}}\text{(}\mathit{\text{k}}\text{= 1,2,...,}\frac{\mathit{\text{m}}\text{-1}}{\text{2}}\text{) .}$$

In this case, the sums of the r^th powers of the roots of Eqs.(8) and (9) are given by$${\mathit{\text{S}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}{\text{= (-1)}}^{\mathit{\text{r}}}\text{+ 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\mathit{\text{m}}\text{-1)/2}} \text{cos}\frac{\text{(2}\mathit{\text{k}}\text{-1)π}\mathit{\text{r}}}{\mathit{\text{m}}}\text{and}{\mathit{\text{T}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{= 1 + 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\text{(}\mathit{\text{m}}\text{-1)/2}} \text{cos}\frac{\text{2}\mathit{\text{k}}\text{π}\mathit{\text{r}}}{\mathit{\text{m}}}\text{,}$$ respectively, so that J_m,n may be written as$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{\text{1}}{\mathit{\text{n}}}\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}\text{[}{\mathit{\text{T}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{-}{\mathit{\text{S}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{].}$$

It follows from Eqs.(10) and (11) that in all cases, the element J_m,n of the matrix J is$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{{\text{(-1)}}^{\mathit{\text{m}}}}{\mathit{\text{n}}}\sum_{=\mathit{\text{r}}\text{=0}}^{\mathit{\text{n}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}{\text{δ}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}$$ where$${\text{δ}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{=}{\mathit{\text{S}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{-}{\mathit{\text{T}}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{,}$$ and S_m,n and T_m,r are the sums of the r^th powers of the roots of Eqs.(8) and (9) respectively. Newton's Formulae [8,10] state that the sums S_r of the r^th powers of the roots of the polynomial equation$${\mathit{\text{z}}}^{\mathit{\text{n}}}\text{+}{\mathit{\text{a}}}_{\mathit{\text{1}}}{\mathit{\text{z}}}^{\mathit{\text{n-1}}}\text{+}{\mathit{\text{a}}}_{\mathit{\text{2}}}{\mathit{\text{z}}}^{\mathit{\text{n-1}}}\text{+...+}{\mathit{\text{a}}}_{\mathit{\text{n-1}}}\mathit{\text{z}}\text{+}{\mathit{\text{a}}}_{\mathit{\text{n}}}\text{= 0}$$ satisfy the relationships

Using these relationships, it can be seen that S_m,r, the sum of the r^th powers of the roots of ${\mathit{\text{z}}}^{\mathit{\text{m}}}{\text{}}^{\text{+1}}\text{=0}$, is given by and the sum T_m,r of the r^th powers of the roots of ${\mathit{\text{z}}}^{\mathit{\text{m}}}\text{-1=0}$ is so that δ_m,r is given by

We note that b_r,n - defined as the coefficient of $\text{cos}\mathit{\text{rx}}$ in the expansion of ${\text{cos}}^{\mathit{\text{n}}}\mathit{\text{x}}$ - will be zero whenever r>n, and furthermore as ${\text{δ}}_{\mathit{\text{m}}}{\text{}}_{\text{,0}}\text{=0 ∀}\mathit{\text{m}}$, we may without prejudice alter the limits of summation in Equ.(12) to run from r=1 to r=N (n≤N):-$${\mathit{\text{J}}}_{\mathit{\text{m}}\text{,}\mathit{\text{n}}}\text{=}\frac{{\text{(-1)}}^{\mathit{\text{m}}}}{\mathit{\text{n}}}\sum_{=\mathit{\text{r}}\text{=1}}^{\underline{\mathit{\text{N}}}} {\mathit{\text{b}}}_{\mathit{\text{r}}\text{,}\mathit{\text{n}}}{\text{δ}}_{\mathit{\text{m}}\text{,}\mathit{\text{r}}}\text{,}$$ which may be written in matrix form as

It is now clear that the matrix J may be written as or where D₁ and D₂ are the diagonal matrices D₁ = diag[-1,1,-1,1,...,(-1)^N] and D₂ = diag[1,½,$\frac{\text{1}}{\text{3}}$,¼,...,¹/_N], Δ is the matrix and B is the matrix Using Equ.(13), the matrix Δ can be written in the form$$\text{Δ =}\mathit{\text{diag}}\text{[-2,-4-6,...,-2}\mathit{\text{N}}{\text{]Δ}}_{\mathit{\text{s}}}$$ where is a sparse matrix such that the elements in the m^th row are 1 in those columns whose column number is an odd multiple of m, and zero otherwise. By combining Eqs.(15) and (16) we obtain where D₃ = diag[2,-4,6,-8,...,2(-1)^N-1N].

It follows from Equ.(18) that the inverse matrix J^-1 is given by Calculating D₂^-1 and D₃^-1 is trivial : D₂^-1 = diag[1,2,3,4,...,N] and D₃^-1 = diag[1/2,-1/4,1/6,-1/8,...,

(-1)^N-1/2N]. As stated above, the elements of B are simply the coefficients of $\text{cos}\mathit{\text{rx}}$ in the expansion of ${\text{cos}}^{\text{n}}\mathit{\text{x}}$. It is easy to show that the elements of the inverse matrix B^-1 are the coefficients of ${\text{cos}}^{\text{r}}\mathit{\text{x}}$ in the expansion of $\text{cos}\mathit{\text{nx}}$ - that is, these elements are the coefficients of the Chebyshev polynomials of the first kind, which are well documented [11]. The element in the i^th row and j^th column of B^-1, t_i,j, is the coefficient of ${\text{cos}}^{\mathit{\text{j}}}\mathit{\text{x}}$ in the expansion of $\text{cos}\mathit{\text{ix}}$, and is given by$${\mathit{\text{t}}}_{\mathit{\text{i}}\text{,}\mathit{\text{j}}}\text{=}\frac{{\text{2}}^{\mathit{\text{i}}\text{-1}}\mathit{\text{j}}\text{(-1)}{\text{}}^{{\text{}}^{\mathit{\text{j}}\text{-}\mathit{\text{i}}}{\text{/}}_{\text{2}}}\left(\frac{\mathit{\text{i}}\text{+}\mathit{\text{j}}}{\text{2}}\text{-1}\right)\text{!}}{\left(\frac{\mathit{\text{j}}\text{-}\mathit{\text{i}}}{\text{2}}\right)\text{!}\mathit{\text{i}}\text{!}}$$ when i≤j and i and j are either both even or both odd, and zero otherwise.

While the inverse matrix Δ_s^-1 is not so easily obtained, we do note that as the elements of Δ_s are all either 0 or 1, the inverse may be easily calculated. Further, as Δ_s^-1 will be an upper triangular matrix, all of whose entries are integers,any numerical inaccuracy in the calculation of the inverse matrix will rapidly become apparent. In summary, we may write the inverse matrix J^-1 in the form where and Δ_s^-1 is the inverse of the matrix given in (17). As all four of the above matrices are rational and upper triangular, it is clear that the matrix J^-1 is also a rational, upper triangular matrix.

We now state the main property of the matrix J^-1. So far, we have considered the inverse of the finite N×N matrix J defined by (4). While the order of the matrix N will in all practical cases be finite, it may, nevertheless be arbitrarily large (and there is no reason why the theory may not be extended to the case where J and J^-1 are infinite matrices). Because J is an upper triangular matrix, the submatrix J_n which consists of the first n rows and columns of J will have an inverse which is the corresponding submatrix - the first n rows and columns - of J^-1. This property means that calculating J^-1 need only be done once, provided of course that it is done for a sufficiently large matrix. All lower orders of the Gradiometer transform can then be dealt with by reading the appropriate submatrix from a 'look-up' table. We provide, in Table I, the first 15 rows and columns of the matrix J^-1.

Examples will now be described of the uses of the Gradiometer Transform.

To evaluate the ability of the Gradiometer transform to extract signal pulses from a slow moving background drift, we consider a known function f(t) which is both piecewise continuous and integrable over the interval [-1,1]. As such, the transform components may be evaluated directly using (1), in which case the component F_m is given by$${\mathit{\text{F}}}_{\mathit{\text{m}}}{\text{= (-1)}}^{\mathit{\text{m}}}\mathit{\text{F}}\text{(1) - 2}\sum_{=\mathit{\text{k}}\text{=1}}^{\mathit{\text{m}}} {\text{(-1)}}^{\mathit{\text{m}}}{\text{}}^{\text{+}\mathit{\text{k}}}\mathit{\text{F}}\left(\text{cos}\frac{\mathit{\text{k}}\text{π}}{\mathit{\text{m}}\text{+1}}\right)\text{,}\mathit{\text{m}}\text{= 0,1,2,...}$$ where$$\mathit{\text{F}}\text{(}\mathit{\text{t}}\text{) =}\int_{\text{-1}}^{\mathit{\text{t}}} \mathit{\text{f}}\text{(}\mathit{\text{t'}}\text{)}\mathit{\text{dt'}}\,d \text{.}$$

The function we have selected is the sum of a background drift and a signal containing three distinct pulses and is shown, along with the corresponding transform components in Fig.2. As expected [7] only the very low order transform components of the drift are significant - implying that the drift can be well approximated by a low order polynomial, whereas the presence of non-zero components over the whole range shows that the signal pulses cannot be successfully fitted by a polynomial of order less that 100.

In Fig.3 we illustrate the effect of removing successively higher order fitted polynomials from our chosen function. In each case the coefficients of a fitted polynomial are obtained from (3) using the transform vector of the combined signal f(t), which is then simply subtracted from f(t). It is apparent that this method of Gradiometer differencing rapidly and effectively eliminates the background drift while leaving each of the signal pulses intact.

In almost all practical applications the signal of interest will consist of sampled data, in which case the transform coefficients are obtained as follows. For sampled data given by$$\text{{ (}{\mathit{\text{t}}}_{\mathit{\text{i}}}\text{,}{\mathit{\text{f}}}_{\mathit{\text{i}}}\text{) :}\mathit{\text{i}}\text{= 1,2,...,}\mathit{\text{N}}\text{} ,}$$ scaled in such a way that -1<t₁<t₂< ... < t_N < 1, we replace the discrete sample with a piecewise continuous rectangular step function given by$${\mathit{\text{f}}}_{\mathit{\text{rect}}}\text{(}\mathit{\text{t}}\text{) =}{\mathit{\text{f}}}_{\mathit{\text{i}}}\text{,}{\mathit{\text{s}}}_{\mathit{\text{i}}}{\text{}}_{\text{-1}}\text{≤}\mathit{\text{t}}\text{<}{\mathit{\text{s}}}_{\mathit{\text{i}}}\text{(}\mathit{\text{i}}\text{= 1,2,...,}\mathit{\text{N}}\text{)}$$ where ${\mathit{\text{s}}}_{\mathit{\text{i}}}\text{= ½(}{\mathit{\text{t}}}_{\mathit{\text{i}}}\text{+}{\mathit{\text{t}}}_{\mathit{\text{i+1}}}\text{)}$ (i = 1,2,...,N-1)

and s₀ = -1, s_N = 1;

we note the data need not be equispaced. The transform component F_m of the sampled data is obtained by replacing f(t) in (1) with f_rect(t). Correlating f_rect(t) with the m^th order gradiometer step function is computationally efficient as it involves simply a change of sign in f_rect at the appropriate gradiometer crossover positions and summing the signed areas of the resulting step function. We provide, in an Appendix, algorithms for calculating F_m for sampled data.

To demonstrate the ability of the Gradiometer Transform to extract signal pulses from noisy data we have sampled a known function - chosen to mimic a variety of typical signal pulses - and added noise. In Fig.4 we illustrate signals at three different noise levels along with the result of using the Gradiometer Transform to fit a 50^th order polynomial to the generated data. Table II gives the true and estimated location of each of the signal peaks marked. It is clear that the reconstructed signals retain both the features and the locations of the original signal pulses, even at moderately high noise levels.

To explore the effectiveness of using the Gradiometer Transform to eliminate noise and/or background drift in experimental data, we have chosen data from the fields of biomagnetism and non-destructive testing. In Fig.5a we show a typical foetal magnetocardiogram [12]. The rapidly varying waveforms in the signal are due to magnetic fields associated with electrical activity in the foetal heart; these are superimposed on a slow moving quasi-periodic drift which is mainly induced by maternal respiration. Because the interval between successive pulses is not strictly periodic, it contains diagnostic information (as do, to a lesser extent, the amplitude and shape of the pulses) [1,12], and so such information must be preserved by any signal processing technique. Calculating the transform components of this signal (Fig.5b) shows an abundance of both low and high order components. While the presence of high order components indicates that high order polynomial fitting is inappropriate here (as may be apparent from the original signal), the dominance of lower order components suggests that low order differencing may be applied in order to eliminate the background drift. Fig.5c shows the effect of removing a 30^th order polynomial from the data; the foetal heartbeat is now clearly visible. The resultant signal is now suitable for standard cardiographic analysis.

In Fig.6a we show the magnetic activity induced at one point on the human skull in response to a nerve stimulus applied to the wrist at time t=0. It is diagnostically important to obtain accurate estimates of the location of the maxima and minima which occur approximately 20ms, 30ms, 60ms and 120ms after the application of the stimulus [1]. While these signal peaks are often visible, determining their exact location is hampered by the noisy nature of the signal and by the possible presence of an underlying magnetic drift, and it is well known that these locations are normally estimated by eye with concomitant experimental bias. It is clearly desirable to quantify this process. Using the Gradiometer Transform we can obtain a polynomial fit to the data in which case the positions of the maxima/minima can be easily determined. In Fig.6b we show the transform components of the signal in 6a; while it is clear that the data may be fitted by a polynomial of order ∼x⁵⁰, the components also show that there is no significant underlying drift to be eliminated here. In Fig.6c we illustrate the 50^th order polynomial fit to the data obtained by using the Gradiometer Transform. A simple Newton-Raphson may be used to obtain the roots of the derivative of this polynomial which will give the location of the peaks in question - in this case the peaks are located at 27.0ms, 33.0ms, 59.1ms and 109.5ms respectively.

Our third data set comes from the field of non-destructive testing. Fig.7a shows the magnetic profile over a steel tube [13]. Each of the pulses corresponds to a defect in the tube, and here it is clear that the signal contains a background magnetic drift (generally present in experimental situations). As the amplitude and shape of each pulse contains information on the shape (i.e. depth and width) of the defect, while the position of the pulse determines its location, it is clear that all three characteristics should be preserved in extracting the pulses from the signal. The Gradiometer transform components of these data are shown in Fig.7b; the dominance of the first two components reflects the clear linear trend in the data. Although not immediately apparent, the data are noisy; the effective absence of high order transform components (|F_m/F₀|<0.0012 ∀ m ≥ 40) shows that in this case it is appropriate to fit a high order polynomial to the data. Figs.7c,d show the result of fitting a 40^th order polynomial to the data and filtering, respectively, linear and 6^th order polynomials; the defects are now clearly visible and as such are amenable to standard analysis.

We have presented in detail the underlying theory of the Gradiometer Transform - a new transform designed to extract signal pulses from data subject to noise and background drift. While calculating the transform components of any given signal (the forward transform) is straightforward, using these to obtain the coefficients of a fitting polynomial (the back-transform) requires knowledge of the generic matrix J^-1. We have presented the underlying theory which allows this matrix to be broken down in such a way that it may be easily and accurately computed, and furthermore shown that its elements are independent of the data set under consideration.

The Gradiometer Transform may be used in a variety of ways to aid signal extraction. Although designed primarily as a tool for extracting pulse signals from background drift [7], the Gradiometer Transform has proved surprisingly effective at extracting signal pulses, by high order polynomial fitting, from noisy data in which there is no drift, while preserving both the general shape and location of the pulses. The transform may also be used simply to eliminate underlying background drift by removing a low order polynomial fit from the original data. In both these cases the Gradiometer Transform will be ∼ 95% computationally more efficient than standard polynomial fitting techniques [7]. We have shown that the Gradiometer Transform is even more efficient when used to simultaneously fit a polynomial to data and remove background drift, as this can be done by filtering in transform space.

The transform components themselves contain useful information about data under consideration, in that these components can be used to indicate the order of fitting and/or filtering that can be reasonably applied, although we note that in all the cases considered here, the precise order of fit or filter is non-critical. We have shown, using data from the fields of biomagnetism and non-destructive testing, hoe the transform components and back-transform may be used to effectively extract the desired signal.

The signal processing method described herein and referred to as the Gradiometer Transform may be applied for practical purposes in a variety of ways, all of which are embraced by the present invention, including the following.

The method may be implemented as a software program running on a general purpose computer. Such a program may operate on "live" or stored data from any applicable source, such as a sensor. An example of this is illustrated in Figure 10, which shows a computer 10, such as a personal computer (PC), for running a program to perform the Gradiometer Transform, connected to a data source 12, which may be a sensor adapted to generate a signal in response to variations in a physical parameter of interest or a source of stored data.

Where the data source 12 is a sensor or the like which outputs an analogue signal, or an analogue data store, then the complete signal processing unit will include analogue to digital converter (ADC) means 14 for sampling the analogue signal or data prior to processing by the computer 10. It will be understood that the ADC means might comprise a separate unit connected between the data source 12 and the computer 10, as illustrated, or may be incorporated into either the data source 12 or the computer 10.

It will also be understood that, where the data source comprises a store of digital data, this store might be external to the computer 10 or may comprise the computer's own data store (typically the hard disk of the computer). The digital data may originate from any applicable source and may be transferred to the data store by any suitable means.

The computer 10 may be adapted to generate an external output 16 from the Gradiometer Transform program which may be employed for any required purpose.

Alternatively, the method may be implemented by a program stored in hardware or firmware in a dedicated data processing unit 16, as is illustrated schematically in Figure 11. Again, the program operates on "live" or stored data from an internal or external data source 20 (sensor, data store or the like, as above), the complete system including ADC means 22 as required. The data processing unit 18 generates an output 24 which may be employed for any required purpose.

Where the system is intended for a specific purpose, then the components of the system may be integrated into a single a single unit comprising, for example, a sensor, ADC means and data processing means. A general purpose unit might comprise the data processing means alone or the ADC and data processing means, with provision for connection to other components as required.

In any of the cases mentioned above, the program implementing the Gradiometer Transform, and other system components, may be dedicated to a particular application or may be user-configurable to suit specific requirements.

The detailed design of software, firmware and hardware components of data processing systems to implement the Gradiometer Transform, as discussed above, is a matter of routine engineering and will vary widely according to the requirements of particular applications. Accordingly, these aspects are not described in detail herein.

Improvements and modifications may be incorporated without departing from the scope of the invention.

References

• 1. Regan and Regan, Fundamentals of electrophysiology, Elsevier (1980).
• 2. H.Weinstock, A review of SQUID magnetometry applied to non-destructive evaluation, IEEE Trans. Magnetics, MAG 27, 3231 (1991)
• 3. T.R.Clune, Intelligent Instruments and Computers 7, 103 (1989)
• 4. K.G.Beauchamp, Walsh Transforms and their applications, Academic Press, London, (1975)
• 5. Wavelets.(J.M.Combes, A.Grossmann, Ph.Tchamitchian Eds.), Springer-Verlag, 1989.
• 6. P.J.Maccabee, N.F.Hassan, Muscle & Nerve, 15, 865 (1992)
• 7. R.J.P.Bain, A.E.Jones, Gradiometer transforms for elimination of polynomial-like background drift. submitted to IEEE Trans. on Signal Processing.
• 8. A.E.Jones, R.J.P.Bain, A generalisation of planar magnetic gradiometer design via orthogonal polynomials, J.Comp.Phys 118,191-199 (1995).
• 9. I.S.Gradshteyn, I.M.Ryzhik (1980), Table of integrals, series and products, Academic Press, New York.
• 10. C.V.Durrell, A.Robson, Advanced Algebra Vol.2, (G.Bell and Sons, London, 1958) p.302.
• 11. M.Abramowitz, I.A.Stegun, Handbook of Mathematical Functions, (Dover, New York, 1970) p.775.
• 12. A.Quinn, A.Weir, U.Shahani, R.J.P.Bain, P.Maas, G.Donaldson, Brit. J. Obstet. Gynaecol 101 886 (1994)
• 13. A.Cochran, G.B.Donaldson, S.Evanson, R.J.P.Bain, IEEE Proc. A, 140 113 (1993)

True location

Estimated location of maxima and minima

a = 0.5

a = 1

a = 2

-0.6000

-0.5978

-0.5974

-0.5965

-0.4925

-0.4928

-0.4905

-0.4835

-0.3091

-0.3061

-0.2888

-0.3356

-0.1000

-0.0829

-0.0878

-0.0916

0.2000

0.2015

0.1968

0.1880

0.5500

0.5611

0.5614

0.5620