WO2006064264A1

WO2006064264A1 - Method of and apparatus for nqr testing

Info

Publication number: WO2006064264A1
Application number: PCT/GB2005/004884
Authority: WO
Inventors: John Alec Sydney Smith; Samuel Somasundaram; Andreas Jakobsson; Magnus Mossberg; Michael David Rowe
Original assignee: King's College London
Priority date: 2004-12-16
Filing date: 2005-12-16
Publication date: 2006-06-22

Abstract

A method of testing comprises irradiating a sample, receiving a response signal and analysing the response signal by combining a plurality of resonance parameters as a function of a variable environmental parameter.

Description

METHOD OF AND APPARATUS FOR NQR TESTING

The present invention relates to a method and apparatus for NQR testing, in particular a method and apparatus for processing of signals received from a sample to analyse the sample or to determine whether a particular compound or substance is present in a sample.

In existing NQR techniques a sample is irradiated by a radio frequency (RF) signal and a radio frequency receiver then "listens" for a response signal. This is typically done either by irradiating the sample using a single pulse and listening for a signal from the sample or by irradiating the sample using a series of pulses and listening for a signal in between them (often referred to as an "echo" technique). The characteristic frequency or frequencies for a large number of compounds containing quadrupolar nuclei have been characterised. A major difficulty, however, arises from the fact that the signals emitted by samples are very weak - in many cases almost indistinguishable from the noise. Another difficulty arises from the fact that the characteristic frequency or frequencies vary with temperature.

The few current publicly available approaches to processing of the NQR signal are mainly based on linear filtering via the Fast Fourier Transform or on matched filtering assuming a reliable estimate of the temperature of the target. These methods are limited due to phase and intensity uncertainties in the NQR signal as well as the difficulty in accurately measuring the temperature of the sample (under ground in the case of mine detection). At present in mine clearance, it can take around half an hour to secure one square meter of ground. Given an accurate temperature estimate, one may combine the dominant frequency responses to a single response with higher signal-to-noise ratio (SNR). See, for example, US patent 6208136. Unfortunately, an imprecise temperature estimate causes the peaks in the frequency domain to be combined sub- optimally. This will potentially result in a failure to detect the mine underground (or a bomb hidden in luggage) with obvious consequences.

According to a first aspect of the present invention, there is provided a method of testing comprising irradiating a sample, receiving a radio frequency response signal and analysing the response signal by combining a plurality of NQR parameters as a function of an variable environmental parameter (or a plurality of variable environmental parameters or a sample- dependent parameter). As a minimum, this may comprise a pair of frequencies or a single frequency and a damping coefficient. Note that there is thus no need to detect the resonance frequencies. In many applications a detection threshold will then be applied to the result to make a "present'V'not present" determination for a compound of interest.

An example of a sample-dependent parameter is the proportion of polymorphic forms in a sample, such as a sample of TNT, or relative intensities of resonance lines at particular resonance frequencies.

The step of combining the plurality of NQR parameters may comprise estimating the value of the NQR parameters and/or estimating the value of the variable environmental parameter.

According to a second aspect of the present invention, there is provided apparatus for testing a sample, the apparatus comprising a radio frequency source, a radio frequency receiver and a parameter combiner operable over a range of a variable environmental parameter.

The environmental parameter may be temperature.

In embodiments of the present invention the search is conducted over at least one other parameter, for example damping coefficient.

In one embodiment of the invention, the technique comprises utilising a nonlinear least squares (NLS) approach, exploiting the fact that the shifts of the spectral lines depend in a known way on temperature; by matching the measured data to the data model formed over a range of possible temperatures, the (unknown) temperature yielding the best match is found. The combined response for this temperature is then used as a detection variable. The NLS method is evaluated using both simulated data, and real NQR data obtained from measurements on a TNT sample. Both these evaluations indicates a strong gain for the proposed method as compared to current state of the art Fourier-based techniques.

In general, in prior art systems, RF interference (RFI) can be a major concern, especially in the detection of TNT where the NQR signal is relatively weak and lies in the AM radio band, therefore being significantly affected by the there present radio transmissions. There are two main prior art approaches to RFI mitigation, passive and active. Passive methods use specially designed antennas, called gradiometers to cancel the far field, the disadvantage being some loss in signal to noise ratio (SNR) of the NQR signal compared to using a simple coil. The active methods employ adaptive noise cancellation techniques which require reference antennas to measure the non- stationary background, the main disadvantage being the cost of extra antennas and electronics needed to measure the background RFI.

Herein, as an illustration, we discuss the detection of TNT. Detecting the masses of TNT typically found in anti-personnel mines poses a great challenge due to low signal-to-noise (SNR) and long data acquisition repeat times. The few publicly available papers discussing the detection of TNT using NQR focus on the use of relatively simple approaches, not fully exploiting the known data model of the NQR echo train.

In another embodiment an Approximative Maximum Likelihood (AML) technique is used which exhibits a significant detection gain over the demodulation method.

In a still further embodiment, a Frequency Selective AML (FSAML) technique is used which has proved to be around three times faster than the AML technique, albeit at some sacrifice in performance.

The underlying idea of the algorithms is to include the unknown temperature (T) (or other environmental parameter or sample dependent parameter) as a parameter to be estimated. We model the data as containing a (known) number of damped sinusoidal components, say D, with frequencies accurately determined by the unknown temperature. The benefit of incorporating a search over T is that we can form a model for the expected data structure for that particular T and then determine how well this particular model fits the data. How well the model fits is given as an energy measure for that given T. For the temperature equal to the unknown true temperature, the data model must fit the best which gives the maximum energy value. Thus, by evaluating how well the data model fits the measured data as a function of temperature, we can determine the target temperature as the value with maximal energy.

To generalize this somewhat further, we also note that the data will depend on the D damping constants of the sinusoidal component. Thus, we can also incorporate a search over each of these constants - and again find the best fit for the damping constants closest to the true values.

In some cases we make an approximation that all the constants have the same value, but this is only a simplification to speed up the processing; the algorithm as such is not restricted to this.

Furthermore, we note that the typical NQR data, with or without the TNT response, will have some particular statistical properties. To incorporate this knowledge, we make a simplified (low-order) model of the signal not containing the TNT response (the background signal), and use this model as a weighting in the estimation. Such a weighting will improve the estimates further. The difference between the AML and the NLS lies only in this weighting as the NLS assumes that the background signal is a completely random signal without any structure.

In a preferred arrangement, the technique is to fit, in a weighted sense, the measured data to the data model using a search over both temperature and damping constants.

The AML and FSAML embodiments may be used in analysis of echo signals.

Both the AML and FSAML detectors when used in analysis of echo signals may exploit the fine structure of the data model within an echo, allowing it to be well modelled as a sum of sinusoids which expand then decay and whose frequencies depend, in a known way, on the temperature of the examined sample.

Those detectors do not allow for loss in magnetisation which, over an echo might be negligible, but significantly affects the echo train when examined in its entirety. The detectors will offer a significant detection gain as compared to current state-of-the-art techniques not exploiting the rich data structure. In a shielded environment, the FSAML detector will offer approximately the same performance as the AML detector, but is computationally cheaper and significantly more robust to likely residual interference. The aforementioned detectors usually pre-process the echo train data by summing sequential echoes to produce a single summed echo with high SNR. However, for cases where long steady-state echo trains cannot be attained (e.g., for TNT), the echoes will decay over the echo train with a rate described by the spin-echo decay time, T_2e (for each resonant frequency) which is in part due to the loss in longitudinal magnetisation, governed by the spin-lattice relaxation time T₁, and in part due to the normal echo decay time, T₂. Consequently, echoes later in the train have a lower SNR then those earlier. Typically, the spin- echo decay time will depend on temperature in a known way for a given experimental set-up. To reflect this dependency on temperature, we use the notation T_2e(τ) where τ denotes the (unknown) temperature of the observed sample.

In further embodiments we incorporate this model knowledge, extending the AML and

FSAML detectors to function on a full echo train. We denote the so-obtained detectors the echo train AML (ETAML) detector and the frequency selective echo train AML (FETAML) detector. There are two particular advantages to modelling the echo train in this way; firstly, the detectors essentially trust earlier echoes more than later ones, allowing for a better fit between the model and the data, and, secondly, knowledge concerning the loss of longitudinal magnetisation is incorporated into the model.

The AML-based detectors may exploit an estimated model of the corrupting noise process. The underlying thermal (Johnson) noise of the RF antenna may be modelled as a white noise process. However, this noise process is shaped by the bandwidth of the receiver, due to the

Quality (Q) factor of the probe and the settings (impulse response) of the anti-aliasing filter, and so in one embodiment we use an approximative low-order autoregressive (AR) model, derived from real noise data.

The detector may be frequency selective (for instance the FSAML detector or the FETAML detector). In that case interference sources, such as, for instance, radio broadcasts, can be excluded from the processing by excluding the frequency grid points of such carriers. This can be done even if these regions are part of the expected frequency regions of interest.

Preferably, inherent dependencies in the observed response signal (for instance, temperature dependencies of resonance frequencies and/or the echo damping constants) are exploited. Different, possibly not a priori known, damping and/or echo damping constants related to one or more resonance frequencies may be incorporated into the model.

Further aspects of the invention are now described, and previously described aspects of the invention are described again or in more detail.

In a further aspect of the invention, very closely related to the first aspect of the invention mentioned above, there is provided a method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by combining a plurality of resonance parameters preferably as a function of a variable environmental parameter.

The step of combining the plurality of resonance parameters may comprise estimating the value of the NQR parameters and/or estimating the value of the variable environmental parameter.

By analysing the response signal in such a way, any relationship between one or more of the resonance parameters and the variable environmental parameter can be exploited, enabling the maximum information to be extracted from the response signal. That is particularly important when signal to noise ratios are low.

The variable environmental parameter may be one of temperature, pressure and magnetic field. Preferably the variable environmental parameter represents an environmental condition to which the sample is subject.

By a variable environmental parameter is meant an environmental parameter which may have any value in a range or set of values. It is possible that the value of the environmental parameter 'would not in fact vary during any particular performance of the method but would take one particular value.

By way of one practical example, in the case where the variable environmental parameter is temperature, the method might be performed, say, at ambient temperature, at, say, an airport. The temperature might then reasonably be considered to have a value somewhere between, say, - 10°C and 35°C, the actual value depending on the ambient conditions at the time the method was performed (and depending on the conditions to which the sample had been subject prior to performance of the method - for instance depending on whether the sample had just been brought out of an aircraft hold following a flight), but the particular value of temperature of the sample might well not change during performance of the method.

In contrast, a non-variable environmental parameter might be considered to be, say, the force of gravity which can be presumed to have the same value in all reasonable circumstances.

Advantageously, this aspect of the invention does not require the value of the environmental parameter to be known or estimated a priori. Rather, the analysis of the response signal may take into account the dependence of one or more of the resonance parameters on the variable environmental parameter. Prior art methods in contrast typically either take no account of the variation of resonance parameters with variable environmental parameter or require the measurement or estimation of the value of the variable environmental parameter. That requirement is both time consuming and potentially disadvantageous as an environmental parameter such as temperature is very difficult to measure or estimate accurately in many practical circumstances, particularly if the sample is concealed or not readily accessible. The response signal is typically a time-dependent signal. The response signal most usually comprises a radio-frequency response signal and the step of irradiating the sample usually comprises applying radio-frequency excitation to the sample.

The response signal may typically be a time varying signal generated in a probe, such as a coil, by receipt of radio-frequency electromagnetic radiation at the probe. Alternatively, the response signal may be generated or derived from such time-varying signal. A number of time varying signals generated in the probe may be accumulated and either summed or signal averaged to form the response signal.

The excitation may comprise pulsed excitation and preferably the excitation comprises a sequence of pulses.

The response signal typically comprises a resonance response signal, and the resonance response signal may be one of a nuclear quadrupole resonance (NQR) response signal, a nuclear magnetic resonance (NMR) response signal or an electron spin resonance (ESR) response signal.

Preferably the resonance parameters comprise at least one of frequency and relaxation time, such as spin-lattice relaxation time or spin-spin relaxation time.

The resonance parameters may comprise a plurality of resonance frequencies, and each of the resonance frequencies may correspond to a respective resonance arising from the same substance. Thus, this aspect of the invention may be used particularly advantageously in testing substances, or testing for the presence of substances, by exciting a plurality of resonances of the substance each at a respective resonance frequency.

The method may be a method for detecting the presence of a substance containing a given species of quadrupolar nucleus. The substance is an explosive or a narcotic and preferably the substance is TNT or RDX.

The method may suitably be used for the detection of a buried or concealed sample, preferably for the detection of a sample concealed in baggage, such as airline baggage. In such circumstances the value of the temperature, or other environmental parameter, is particularly difficult to measure or estimate a priori and so prior art methods may be particularly disadvantageous. The plurality of resonance parameters may be combined as a function of a plurality of variable environmental parameters.

The step of analysing the response signal may comprise selecting a value of the variable environmental parameter, preferably in dependence upon the response signal.

The selected value may be selected as being a value which is consistent with the response signal and/or the values of the other parameters. In that case, the selected value of the variable environmental parameter is not necessarily the actual value of the variable environmental parameter (for instance in analysing the response signal a value of temperature may be selected which is not necessarily the actual temperature of the sample), but if the analysis of the response signal is accurate it should at least be close to that actual value of the variable environmental parameter.

Preferably the step of analysing the response signal comprises, for at least one of the resonance parameters, selecting a value of that resonance parameter, preferably in dependence upon the response signal.

The step of analysing the response signal may comprise using a model of the response signal, the model combining the plurality of resonance parameters, and preferably at least one of the resonance parameters is a function of a variable environmental parameter.

The model may be a model of response signal as a function of time. Preferably the model is a model of the amplitude of the response signal as a function of time, or a model of the real and/or imaginary part of the response signal as a function of time.

Typically the form of the model of the response signal as a function of time would be dependent upon the values selected for the resonance parameters and the value selected for the variable environmental parameter.

Preferably the step of using the model comprises fitting the model to the response signal, preferably by selecting a value of the variable environmental parameter and/or a selecting a value of at least one of the resonance parameters. That feature is particularly important and so in a further aspect there is provided a method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by fitting a model of the response signal to the response signal, the model comprising a plurality of resonance parameters each being a function of a variable environmental parameter, and preferably the step of fitting the model of the response signal to the response signal comprises selecting a value of each of the resonance parameters and selecting a value of the variable environmental parameter.

Preferably the step of selecting a value of the at least one resonance parameter and/or the value of the variable environmental parameter is carried out in dependence on a pre-determined relationship between one or more of the resonance parameters and the variable environmental parameter and/or a pre-determined relationship between the plurality of resonance parameters.

Thus, advantageously, prior knowledge of the likely or possible behaviour of the sample can be used in analysing the response signal enabling information to be obtained more efficiently and accurately from the response signal.

Preferably the step of selecting a value of the at least one resonance parameter is carried out in dependence on a pre-determined relationship between one or more of the resonance parameters and the excitation applied to the sample. For instance, if one of the resonance parameters is a relaxation time, and the excitation comprises a sequence of pulses, the relaxation time may, in certain circumstances, be dependent on the spacing in time of the pulses in the sequence of pulses.

Preferably the pre-determined relationship is a relationship between resonance frequency and temperature, and is preferably a linear relationship.

In the case of NQR resonance frequencies the relationship between resonance frequency and temperature has been discussed in the prior art (see for instance Advances in Nuclear Quadrupole Resonance, VoIA, pi -69, 1980 which is herein incorporated by reference. Reference is made in particular to pages 2 to 13 of that document, including equation (9) on page 6). In a certain temperature range, NQR resonance frequency varies linearly with temperature. Preferably the method further comprises selecting a temperature range and using a pre-determined relationship between at least one of the resonance parameters, preferably resonance frequency, and temperature suitable for the selected temperature range. Alternatively or additionally the pre-determined relationship may be a relationship between decay constant and temperature.

The method ^'may comprise selecting the value of at least one further parameter of the model.

The method may comprise selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter'in dependence upon how well the response signal fits the model.

The step of selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter may be carried out using an iterative procedure.

The iterative procedure may comprise varying the value of the variable environmental parameter and/or the value of the or each resonance parameters and/or the value of the or each further parameter and comparing the model to the response signal until the model is determined to match the response signal and/or until a pre-determined number of iterations have been performed.

Preferably the step of comparing the model to the response signal comprises generating a measure of how well the model matches the response signal.

The model may be determined to match the response signal when the measure of how well the model matches the response signal is within a pre-determined limit.

The method may comprise selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter using a least squares technique, preferably a non-linear least squares technique.

The method may comprise selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter using a Maximum Likelihood Estimation technique. The Maximum Likelihood Estimation technique may be a frequency selective Maximum

Likelihood Estimation technique or a time selective Maximum Likelihood Estimation technique.

Preferably the step of selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter comprises selecting the value of one of the parameters, fixing the value of that parameter at the selected value and then selecting the value of another of the parameters.

The step of selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter may comprise selecting the value of at least two of the parameters simultaneously.

The method may comprise selecting, for at least one of the resonance parameters and/or for the environmental parameter and/or for the at least one further parameter, a range or set of values and selecting the value of that parameter from that range or set of values.

The method may comprise selecting a portion of the response signal and performing the analysis step independently of any analysis of any other portion of the or a response signal, and preferably the method comprises performing the analysis step only on the selected portion of the response signal.

Preferably the step of selecting a portion of the response signal comprises selecting the portion of the response signal being within a selected range or set of frequencies.

The step of selecting a portion of the response signal may comprise selecting the portion of the response signal being within a selected range or set of times.

The step of using the model may comprise selecting the components of the model.

The model may comprise at least one component representing a free induction decay and/or at least one component representing an echo. The model may comprise at least one component representing an echo decay.

The model may comprise at least one component representing a train of echoes. In a further aspect there is provided a method of, testing comprising irradiating a sample, receiving a response signal and analysing the response signal by identifying a component of the response signal representing a train of echoes.

Preferably the component representing the train of echoes comprises a component representing the decay of the train of echoes. The component representing the decay of the train of echoes may comprise a component representing the decay of peak echo amplitude.

Preferably in analysing the response signal more significance is given to earlier echoes in the train of echoes than to later echoes in the train of echoes..

Preferably the' model comprises at least one component representing a steady state signal associated with a train of echoes.

The model may comprise a component which represents the response signal as at least one decaying sinusoid, and preferably the or each decaying sinusoid corresponds to a respective resonance. .

At least one of the resonance parameters may comprise a decay constant of the or at least one of the decaying sinusoids.

Preferably the model comprises at least one component representing an undesired signal. The model may comprise at least one component representing radio-frequency interference. The model may comprise at least one component representing a noise signal.

The noise signal may comprise a non-white noise signal.

In a further aspect there is provide a method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by identifying a component of the response signal representing a resonance response and by identifying a component of the response signal representing a non-white noise signal.

The method may comprise determining at least one characteristic of the apparatus used to perform the method and selecting the component representing the noise signal in dependence upon the characteristic. Preferably the method comprises performing at least one measurement using the apparatus and determining the at least one characteristic in dependence upon the at least one measurement.

The method may comprise comparing the response signal to a threshold and preferably generating an alarm signal in dependence upon the comparison.

The step of comparing the response signal to a threshold may comprise generating an output in dependence upon the response signal and comparing the output to the threshold.

Alternatively or additionally, the method may comprise comparing the model to a threshold and preferably generating an alarm signal in dependence upon the comparison.

The step of comparing the model to a threshold may comprise generating an output in dependence upon the model and comparing the output to the threshold.

The output may comprise a test statistic. The test statistic may represent the likelihood of the response signal including a resonance response signal.

The method may comprise generating the test statistic according to a generalised likelihood ratio test.

The step of comparing the model to a threshold may comprise comparing a component of the model representing a resonance response to the threshold.

The step of comparing the model to a threshold may comprise comparing the value of at least one of the resonance parameters to a threshold and preferably generating an alarm signal in dependence upon the comparison.

The method may further comprise analysing the composition of the sample in dependence upon the analysis step and preferably the sample comprises a pharmaceutical. The method may further comprise outputting the value of one or more of the parameters and preferably using the outputted values in a further measurement or in a further analysis procedure.

In a further aspect of the invention there is provided apparatus for testing a sample, comprising means (for instance a probe) for irradiating the sample, means (for instance a receiver) for receiving a response signal (for instance a receiver) and means (for instance a processor) for analysing the response signal by combining a plurality of resonance parameters preferably as a function of a variable environmental parameter.

Other means, for instance as mentioned below, incorporated in or associated with the analysis means (for instance a processor) may be in the form of suitable software programs or modules for execution on the analysis means (for instance a processor).

The means for analysing the response signal may be implemented in hardware or software, and may comprise a processor. In certain embodiments, the processor may be included in a control computer used for controlling the testing.

The means for analysing the response signal may be adapted to select a value of the variable environmental parameter, preferably in dependence upon the response signal.

The means for analysing the response signal may be adapted to select, for at least one of the resonance parameters, a value of that resonance parameter, preferably in dependence upon the response signal.

Preferably the means for analysing the response signal comprises means for using a model of the response signal, the model combining the plurality of resonance parameters, and preferably at least one of the resonance parameters is a function of a variable environmental parameter.

Preferably the means for using the model is adapted to fit the model to the response signal, preferably by selecting a value of the variable environmental parameter and/or by selecting a value of at least one of the resonance parameters. In a further aspect, there is provided apparatus for testing, comprising means (for instance a probe) for irradiating a sample, means (for instance a receiver) for receiving a response signal and means (for instance a processor) for analysing the response signal by fitting a model of the response signal to the response signal, the model comprising a plurality of resonance parameters each being a function of a variable environmental parameter, and preferably the means for analysing the response signal is adapted to fit the model of the response signal to the response signal by selecting a value of each of the resonance parameters and selecting a value of the variable environmental parameter.

Preferably the means for analysing the response signal is adapted to select a value of the at least one resonance parameter and/or the value of the variable environmental parameter in dependence on a pre-determined relationship between the resonance parameter and the variable environmental parameter and/or a pre-determined relationship between the plurality of resonance parameters.

The pre-determined relationship may be a relationship between resonance frequency and temperature, and is preferably a linear relationship.

The pre-determined relationship may be a relationship between decay constant and temperature.

The means for analysing the response signal may be adapted to select the value of at least one further parameter.

The means for analysing the response signal may be adapted to select the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter in dependence upon how well the response signal fits the model.

The means for analysing the response signal may be adapted to select the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter using an iterative procedure.

The means for analysing the response signal may be adapted to compare the model to the response signal by generating a measure of how well the model matches the response signal.

The means for analysing the response signal may be adapted to select the value of the variable environmental parameter and/or .the value of the or each resonance parameter and/or the value of the or each further parameter using a least squares technique, preferably a non-linear least squares technique.

The means for analysing the response signal may be adapted to select the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter using a Maximum Likelihood Estimation technique.

The Maximum Likelihood Estimation technique may be a frequency selective Maximum

Preferably the means for analysing the response signal is adapted to select the value of one of the variable environmental parameter, the or each resonance parameter and the or each further parameter, fix the value of that parameter at the selected value and then select the value of another of the parameters.

Preferably the means for analysing the response signal is adapted to select simultaneously the value of at least two of the variable environmental parameter, the or each resonance parameter, and the or each further parameter.

Preferably the means for analysing the response signal is adapted to select, for at least one of the resonance parameters and/or for the environmental parameter and/or for the at least one further parameter, a range or set of values and to select the value of that parameter from that range or set of values. Preferably the apparatus comprises means for selecting a portion of the response signal and preferably the means for analysing the response signal is adapted to analyse the selected portion of the response signal independently of any analysis of any other portion of the or a response signal, and preferably is adapted to analyse only the selected portion of the response signal.

Preferably the analysis means is adapted to select a portion of the response signal being within a selected range or set of frequencies or being within a selected range or set of times.

Preferably the apparatus comprises means for selecting the components of the model.

Preferably the model comprises at least one component representing a free induction decay.

Preferably the model comprises at least one component representing an echo. The model may comprise at least one component representing an echo decay, and the model may comprise at least one component representing a train of echoes.

In a further aspect of the invention there is provided apparatus for testing a sample, comprising means (for instance a probe) for irradiating the sample, means (for instance a receiver) for receiving a response signal and means (for instance a processor) for analysing the response signal by identifying a component of the response signal representing a train of echoes.

The component representing the train of echoes may comprise a component representing the decay of the train of echoes. The component representing the decay of the train of echoes may comprise a component representing the decay of peak echo amplitude.

The model may comprise at least one component representing a steady state signal associated with a train of echoes.

The model may comprise a component which represents the response signal as at least one decaying sinusoid, and preferably the or each decaying sinusoid corresponds to a respective resonance. At least one of the resonance parameters may comprise a decay constant of the or at least one of the decaying sinusoids.

The model may comprise at least one component representing an undesired signal.

The model may comprise at least one component representing radio-frequency interference.

The model may comprise at least one component representing a noise signal.

The noise signal may comprise a non-white noise signal.

In a further aspect of the invention there is provided apparatus for testing a sample, comprising means (for instance a probe) for irradiating a sample, means (for instance a receiver) for receiving a response signal and means (for instance a processor) for analysing the response signal by identifying a component of the response signal representing a resonance response and by identifying a component of the response signal representing a non-white noise signal.

The apparatus may comprise means for selecting the component representing the noise signal in dependence upon a characteristic of the apparatus.

The apparatus may comprise means for determining the at least one characteristic of the apparatus in dependence upon at least one test measurement.

Preferably the apparatus comprises comparing means for comparing the response signal to a threshold and preferably comprises means for generating an alarm signal in dependence upon the comparison.

The comparing means may be adapted to generate an output in dependence upon the response signal and to compare the output to the threshold.

The comparing means may, additionally or alternatively, be adapted to compare the model to a threshold and preferably the means for generating an alarm signal is adapted to generate an alarm signal in dependence upon the comparison. Preferably the analysis means is adapted to generate an output in dependence upon the model and the comparing means is adapted to compare the output to the threshold.

Preferably the output comprises a test statistic.

The test statistic may represent the likelihood of the response signal including a resonance response signal.

Preferably the analysis means is adapted to generate the test statistic according to a generalised likelihood ratio test.

Preferably the comparing means is adapted to compare a component of the model representing a resonance response to the threshold.

Preferably the comparing means is adapted to compare the value of at least one of the resonance parameters to a threshold and preferably the apparatus comprises means for generating an alarm signal in dependence upon the comparison.

The variable environmental parameter may be one of temperature, pressure and magnetic field.

The response signal may be a time-dependent signal. In particular, the response signal may comprise a radio-frequency response signal and/or means for irradiating the sample may be adapted to apply radio-frequency excitation to the sample.

The response signal may comprises a resonance response signal, the resonance response signal being one of a nuclear quadrupole resonance (NQR) response signal, a nuclear magnetic resonance (NMR) response signal or an electron spin resonance (ESR) response signal.

Preferably the resonance parameters comprise at least one of frequency and relaxation time. In particular, the resonance parameters may comprise at least one of spin-lattice relaxation time and spin-spin relaxation time. The resonance parameters may comprise a plurality of resonance frequencies.

Each of the resonance frequencies may correspond to a respective resonance arising from the same substance.

Preferably the apparatus is adapted to detect the presence of a substance containing a given species of quadrupolar nucleus. The substance may be an explosive or a narcotic and preferably the substance is TNT or RDX.

Preferably the apparatus is adapted to detect the presence or absence of a buried or concealed sample, preferably a sample concealed in baggage.

The apparatus may be adapted to analyse the composition of the sample and preferably the sample comprises a pharmaceutical.

Preferably the apparatus is adapted to output the value of one or more of the parameters and preferably there is provided means for using the outputted values in a further measurement or in a further analysis procedure.

Background information concerning analysis of signals may be found in the following:- P. Stoica and R. Moses, Spectral Analysis of Signals, Prentice Hall, Upper Saddle River, NJ. 2005; S. Umesh and D. W. Tufts, "Estimation of Parameters of Exponentially Damped Sinusoids Using Fast Maximum Likelihood Estimation with Application to NMR Spectroscopy Data" IEEE Trans Signal Processing, vol. 44, no. 9, pp 2245 - 2259, September 1996; and S.M. Kay, Fundamentals of Statistical Signal Processing, Volume II: Detection Theory. Prentice Hall, Englewood Cliffs, NJ. 1998.

Aspects of the invention also provide for an apparatus for performing the above-mentioned independent and preferred methods of the invention.

The invention also provides a computer program and a computer program product for carrying out any of the methods described herein and/or for embodying any of the apparatus features described herein, and a computer readable medium having stored thereon a program for carrying out any of the methods described herein and/or for embodying any of the apparatus features described herein.

The invention also provides a signal embodying a computer program for carrying out any of the methods described herein and/or for embodying any of the apparatus features described herein, a method of transmitting such a signal, and a computer product having an operating system which supports a computer program for carrying out any of the methods described herein and/or for embodying any of the apparatus features described herein.

The invention extends to methods and/or apparatus substantially as herein described with reference to the accompanying drawings.

Any feature in one aspect of the invention may be applied to other aspects of the invention, in any appropriate combination. In particular, method aspects may be applied to apparatus aspects, and vice versa.

The present invention will now be described by way of non-limiting example with reference to the accompanying drawings, in which:

Figure 1 shows diagrammatically the ¹⁴N quadrupole energy levels; Figure 2 shows a pair of graphs illustrating the current state of the art in detection - the demodulation approach;

Figure 3 is a block diagram of apparatus for NQR testing;

Figure 4 shows the variation in frequency (arbitrary units) with temperature (Kelvin, horizontal axis) for a compound exhibiting nuclear quadrupole resonance; Figure 5 shows the performance gain as a function of the measured echo number for a real NQR signal using the NLS detector;

Figure 6 shows the performance gain as a function of the SNR for a simulated NQR signal using the NLS detector;

Figure 7 shows the receiver-operator curve (ROC) for real NQR data using the NLS detector;

Figure 8 shows two graphs of the two dimensional search conducted using the AML technique, here approximated as two, one-dimensional searches;

Figure 9 shows the results for frequency selective AML over the parameters of temperature and damping constant; Figure 10 shows the detection gain for AML,, and FSAML compared with prior art techniques;

Figure 11 and 12 show probability of false alarm versus probability of detection for AML and FSAML versus the prior art techniques at two different, very low signal to noise ratios; Figure 13 is an illustration of the real part of a typical echo train (the dashed line showing the overall damping of the signal reflects the η(τ) damping over consecutive echoes);

Figures 14a and 14b are plots of detection gain for a sample containing TNT and for a sample not containing TNT (from a set of 100 files each containing 20 summed echo trains);

Figures 15a and 15b are plots of detection gain for simulated NQR signals; ^' Figure 16a is a plot of receiver operator characteristic (ROC) curves for AML-based detectors, for simulated data without RFI components, at -32dB;

Figure 16b is ' a plot of receiver, operator characteristic (ROC) curves for DMA-based detectors, for simulated data without RFI components, at-32dB;

Figure 17a is a plot of receiver operator characteristic (ROC) curves for the ETAML, ETAML-s and ETAML-a detectors, without RFI components, at -33dB;

Figure 17b is a plot of receiver operator characteristic (ROC) curves for the FETAML detector with different amounts of zero padding, for simulated data without RFI, at -32dB;

Figure 18a is a plot of receiver operator characteristic (ROC) curves for the AML-based detectors, for simulated data with RFI components, at -3OdB; Figure 18b is a plot of receiver operator characteristic (ROC) curves for the DMA-based detectors, for simulated data with RFI components, at -3OdB;

, Figures 19a and 19b are plots of receiver operator characteristic (ROC) curves for the ETAML-s, FETAML-s, FSAML and AML detectors, for partially shielded measured data;

Figures 20a and 20b are plots of receiver operator characteristic (ROC) curves for the ETAML-s, FETAML-s, FSAML and AML detectors, for partially shielded measured data; and

Figures 21a and 21b are plots of receiver operator characteristic (ROC) curves for the DMA-based detectors, for partially shielded measured data.

Embodiments of the present invention may be used in relation to various resonance techniques, including nuclear quadrupole resonance (NQR), nuclear magnetic resonance (NMR), magnetic resonance imaging (MRI), and electron spin resonance (ESR). In the following, consideration is given in particular to use of embodiments of the invention in relation to nuclear quadrupole resonance.

The following description is structured as follows:- • firstly an introduction to nuclear quadrupole resonance is provided, and its relationship to both nuclear magnetic resonance and magnetic resonance imaging is discussed briefly

• an example of NQR test apparatus is described • a data model of free induction decay is derived and discussed

• next, various detectors based upon the free induction decay model are described (the first , using a non-linear least squares technique, the second using an Approximate Maximum

Likelihood (AML) technique, and the third using a frequency selective' Approximate Maximum Likelihood (FSAML) technique) and their performance discussed • then, data models for echoes and echo trains are derived and discussed

• various detectors based upon the echo train data model are described (the first using an Approximate Maximum Likelihood (AML) technique and being referred to as the ETAML detector and the second using a frequency selective Approximate Maximum Likelihood (FSAML) technique and being referred to as the FETAML detector) • finally, the performance of the ETAML and FETAML detectors is examined using both simulated and real NQR data.

Introduction to nuclear quadrupole resonance (NQR)

Nuclear Quadrupole Resonance (NQR) is a pulsed radio frequency (RF) technique which can be used to detect the presence of quadrupolar nuclei. The technique has recently attracted attention from the aviation industry as a method for screening baggage for narcotics and explosives, the pharmaceutical industry as a way to characterise various drugs, and the military as a technique for detecting landmines and unexploded ordnance.

NQR is related to both nuclear magnetic resonance (NMR) and magnetic resonance imaging (MRI), but unlike NMR and MRI, does not require a large static magnetic field to split the energy levels of the nucleus, making it attractive as a non-invasive technique for, for instance, landmine detection and a relatively inexpensive technique for the characterisation of various compounds. NQR can only be applied to compounds containing a quadrupolar nucleus, such as the ¹⁴N nucleus. A quadrupolar nucleus behaves as though it has a non-spherical charge distribution and therefore possesses an electric quadrupole moment. When this type of nucleus is placed in a nonzero ejectric field gradient (EFG), different quadrupole energy levels arise as different nuclear orientations are energetically more favourable than others. The EFG seen by the nucleus is a result of the neighbouring charges, both electrons and nuclei, and is directly related to the chemical structure of the compound. The compound specific nature of NQR is due to this relation between the EFG and the chemical structure of the compound, and its uniqueness means that there is little or no interference from other nitrogenous compounds that may be present. The symmetry of the EFG, the magnitude of the nuclear electric quadrupole moment and the spin-quantum number of the nucleus all determine how the energy levels will be split. The ¹⁴N nucleus has spin-quantum number /= 1 and the energy levels, obtained by solving the quadrupolar Hamiltonian for a spin-1 nucleus, are shown in Figure 1 with the three allowed transitions v₊, v. and V₀. The signals are acquired by applying pulsed RF radiation which drives transitions between the quadrupolar energy levels, and then measuring the response. Commonly, two types of signal are measured, the free induction decay (FID), which is the signal obtained immediately after an excitation pulse, and echoes, which are the signals obtained between a. string of pulses, the latter having the advantage that a larger number of useful signals can be collected in a given time.

Figure 2 illustrates the state of the art demodulation approach to analyzing NQR FID response_, signals. The drawback with that approach is that an accurate knowledge of temperature is assumed and this parameter is usually not known sufficiently well, as discussed in more detail below.

NQR test apparatus

An NQR test apparatus will consist of a transmitter (to irradiate the sample), a receiver (to listen for response signals), signal processing circuitry (to analyse the response signal) and a controller (for timing of, for instance, pulses used to irradiate the sample, and for selection of amplitude and phase of the pulses). A wide variety of apparatus is known for NQR testing, varying in accordance with whether it is intended for mine detection, baggage screening and so on. The present invention is applicable to all such apparatus - one example of which is shown and described in US 6208136 incorporated herein by reference.

Similarly, the present invention is also applicable to signals excited by any one of a range of pulse sequences as are known in the art. Examples of such sequences include single-frequency excitation, phase cycling, PAPS, NPAPS and so on.

Description of the apparatus in one embodiment is now described. Referring first to Figure 3, apparatus for NQR testing includes a radio-frequency source 11 connected via a phase/amplitude control 10 and a gate 12 to an r.f. power amplifier 13. The output of the latter is connected to an r.f. probe 14 which contains one or more r.f. coils disposed about or adjacent the sample to be tested (not shown), such that the sample can be irradiated with r.f. pulses at the appropriate frequency or frequencies to excite nuclear quadrupole resonance in the substance under test (for example, an explosive). The r.f. probe 14 is also connected to r.f. receiver and detection circuitry 15 for detecting nuclear quadrupole response signals. The detected signal is sent from circuitry 15 to a control computer 16 (or other control apparatus) for processing, and for signal addition or subtraction. The computer includes some means 17 for producing an alarm signal in dependence upon whether a given threshold of detection for the presence of the particular substance of interest has been exceeded. The alarm signal would normally be used to activate an audio or visual alarm to alert the operator to the presence of the substance under test.

The control computer 16 also controls all pulses, their radio frequency, time, length, amplitude and phase. In the context of the present invention all of these parameters may need to be adjusted precisely; for example, phase may need to be varied in order to be able to generate echo responses.

Re-tuning of the r.f. probe 14, alteration of its matching and alteration of its Q factor may all need to be carried out dependent upon the nature of the sample. These functions are carried out by the control computer 16 as follows. Firstly, the computer checks the tuning of the r.f. probe 14 by means of a pick-up coil 18 and r.f. monitor 19, making adjustments by means of the tuning control 20. Secondly, the matching to the r.f. power amplifier 13 is monitored by means of a directional coupler 21 (or directional wattmeter), which the computer responds to via a matching circuit 22, which in turn adjusts the r.f. probe 14 by means of a variable capacitance or inductance. The directional coupler 21 is switched out by the computer 16 when not required, via switch 23. Thirdly, the Q factor of the r.f. coil is monitored by a frequency-switch programme and adjusted by means of a Q-switch 24 which either changes the coil Q or alternatively alerts the computer to increase the number of measurements.

The control computer 16 may be programmed in various ways to reduce or eliminate the spurious interference described above by controlling the pulse amplitudes and phases by means of the control 10. These ways can involve the use of a comparator 25 for comparing the response signals from different pulses by making appropriate changes to the phase of the receiver and detection circuitry 15, and passing the resultant signals to the remainder of the control computer 16 for further processing.

The control computer is also used to analyse the response signals. The methods used for analysis of the response signals are described in more detail later. The control computer may also be used to sum and/or signal average response signals. Analysis may be carried out on individual response signals or on summed response signals or on averaged response signals.

Typically, appropriate software is run on the control computer in order to carry out one or more of the analysis methods described below.- In the preferred embodiment the user may select which analysis method is to be used, depending on the nature of the testing to be carried out.

As described in more detail below the analysis may comprise estimating the values of resonance parameters and/or one or more variable environmental parameters (for instance temperature), the estimation in the preferred embodiment being performed by fitting one or more response signals to a model. The control computer in fitting the response signal or signals to the model, selects values of the various parameters included in the model (including, in the preferred embodiment, selecting a value of temperature based upon the response signal or signals) and, in certain variants, also selects the form of the model and the components to be included in the model. Typically the control computer uses an iterative fitting procedure to fit the response signal or signals to the model. In particular uses of the preferred embodiment, the results of the analysis procedure, for instance one or more of the fitted values of the parameters of the model, are output. The outputted results may be used, for instance in a further analysis procedure (for instance in analysing, say, the purity of a sample) or in setting up a further measurement (for instance the fitted value of, say, temperature may be used in setting up another measurement on the sample, for instance a measurement using a different experimental technique).

The control computer may be set up so that the fitting of the response signal or signals is performed with a pre-determined level of accuracy.

In the case where the apparatus is used for the detection of the presence or absence of a given substance, it may not be necessary to find the best fit to the response signals; an approximate fit may be suitable (in that case, typically, the response signal may be fitted using a series of 1- dimensional fits of each parameter in turn rather than a full multi-dimensional fit of all parameters simultaneously).

In the case where the apparatus is used for analysis of a characteristic (for instance the purity) of the sample (or a particular substance included in the sample) under test then it may be important to find the best fit giving the most accurate values of the parameters. In that case, the response signal may be fitted using a full multi-dimensional fit of all parameters simultaneously.

Methods for searching for best or estimated fits are known in the art and may be used. Alternatively methods for searching for best or estimated fits described explicitly herein may be used.

In alternative embodiments response signals are passed to a further computer or processor for analysis. The analysis may be performed by dedicated circuitry or by suitable software. Shown diagrammatically in Figure 3 and designated as 27 is some means, such as a conveyor belt, for transporting a succession of samples to a region adjacent the r:f. probe 14. The computer 16 is arranged to time the application of the excitation pulses substantially simultaneously with the arrival of a particular sample adjacent the probe. In alternative embodiments, instead of the sample being carried on a conveyor belt, it may actually be a person, and the r.f. probe may be in the form of a walk-through gateway or a hand-held wand.

Although the apparatus described above would usually employ rectangular pulses, other pulse shapes may be employed. Furthermore although usually the radio-frequency probe would utilise a single coil for both transmission and reception of signals, any appropriate number of coils may be used, and different coils can be used for transmission and reception. Also, the apparatus would usually operate in the absence of any applied magnetic field.

Now we consider the data model of NQR FID responses as a background to the signal processing techniques that are exploited in the various embodiments of the present invention in relation to FID responses. As discussed in more detail later in some embodiments the data model of NQR FID responses is also used to model echo responses.

Data model - Free Induction Decay (FID)

Experimentally, we apply pulsed RF radiation at the proper frequency and geometry in order to drive transitions between the quadrupole energy levels. In a single pulse experiment where a single resonant line is excited, the magnetic component of the applied radiation couples with the nuclear magnetic moment and rotates the nucleus away from its equilibrium position in the electric field gradient (EFG), exciting the system. The system returns to equilibrium via various relaxation processes inducing a measurable signal, known as a free induction decay (FID), in the receiving antenna. The FID is observed to decay with time constant T^ , known as the FID or spin-phase memory decay time. Assuming that d resonant lines are excited, the measured signal can be well modelled as a sum of d damped sinusoids. The relative amplitudes of these d frequencies can be assumed to be known for a given experimental set-up. Relative amplitudes, between experimental setups, may change due to:

1. Variations in excitation bandwidth - caused by varying the excitation pulse width and the excitation frequency. 2. Changes in receiver bandwidth - caused by altering the Q-factor of the probe or the impulse response of the anti-aliasing filter.

3. Differing experimental repeat times, which may mean some of the ¹⁴N nuclei may not have fully relaxed before an excitation pulse is applied.

As mentioned above, the NQR signal can be well modeled as a sum of d damped sinusoids (see, for example, the prior patent referred to above)

where ot_k and, β_k denote the (complex) amplitude and the damping constant of the Mi sinusoid, respectively. We note that sampling can be made in an arbitrary, possibly non-uniform way, with t=to_> ... t_N-i. For simplicity, we will at this stage assume uniform sampling strategy starting at t_o=l, i.e. t=l, ..., N. Typically, all the spectral lines will have approximately the same damping constants, say βo, which may not vary significantly with temperature, but may vary between samples. In one embodiment, we do not exploit this fact in an effort to allow for cases when the damping constants are measurably different; the resulting general algorithm then simplifies in a natural way if we let β_k ~ β₀. Further, ω_k(τ) is the frequency shifting function of the

Mh sinusoidal component due to the (unknown) temperature of the explosive sample (where T represents the temperature) and ω(t) is an additive colored noise.

We note here that β_k may be related to the spin-phase memory decay time of the k^th frequency, here denoted T _2,k by the equation β_k=D_w/T _2,k, where D_w is the spectrometer dwell time.

One of the relaxation processes by which the nucleus returns to its equilibrium orientation in the EFG is by losing energy to the thermal motions of the solid. This process is governed by the spin- lattice relaxation time, here denoted 7";. It is worth noting that a wait time of 5T] is normally required to produce a fully relaxed system, i.e., after the acquisition of an FID, one must wait 5Ti before applying another excitation pulse. This is a limiting factor in the detection of TNT as 5Ti can be up to 30 seconds. One way to reduce this limitation is to use multiple pulse sequences to generate a train of echoes. The use of echo trains is considered further in subsequent sections.

Remaining with the FID model, an important point to note is that the number of damped sinusoids, as well as the frequency shifting function for each spectral line, ω_k(τ), may be assumed to be biown, whereas a^ β_k, as well as the temperature of the explosive sample, T, are unknown. For NQR signals of many explosive samples, particularly TNT, the frequency shifting function at likely temperatures can be well modeled by equation (2) - see for example PCT patent publication WO99/45409.

¹ where Ci]₀ and β_k, for k = 1, ..., d, are given constants. Often, the relative ratio between the modulus of the signal amplitudes, | (X_k | , are approximately known for a given explosive sample. To exploit this knowledge, we will let o^ = _pkδ_k, where δk denotes the a priori known scaling.

The variation of NQR v₊ frequencies with temperature in the monoclinic form of TNT is illustrated in Figure 4. Each of the lines on the graph relates to a single response frequency ("line" in the frequency domain) from the compound. It will be noted that the variation of frequency with temperature, though predictable, is somewhat different between the different responses. The consequence of a poor estimate of sample temperature can be visualised using this graph. If the NQR apparatus is searching for responses for a particular temperature but the sample is, in fact, at a different temperature, then the QR material will go undetected.

For completeness, Table 1 summarises the values of a_k, b_k and a^ for the four line region of monoclinic TNT, using an excitation frequency of 841.5 KHz, in the region of 830-860 KHz. In the table, a^ is referred to as K]₁. In practice, the particular values will depend on the actual experimental set-up.

Table 1

It should be noted that the present and later-described embodiments can also be used even in temperature regions where NQR resonance frequencies do not vary linearly with temperature. In that case the linear frequency shifting function expressed in equation (2) is replaced with a frequency shifting function appropriate for the temperature region in question. The temperature dependence of NQR resonance frequencies is discussed, for instance, in Advances in Nuclear Quadrupole Resonance, VoI 4, pl-69, 1980 which is herein incorporated by reference. Reference is made in particular to pages 2 to 13 of that document, including equation (9) on page 6.

Detector based upon the free induction decay model using a non-linear least squares technique

In a first embodiment of the present invention a nonlinear least squares technique is exploited. It is described here in mathematical terms but the skilled person will readily be capable of implementing this in data processing hardware and/or software, for instance on the control computer in the preferred embodiment.

For derivation purposes, we will herein model ω(t) as a complex white Gaussian noise. This might be thought to be a crude approximation, but we note that the NLS estimator will for sinusoidal estimation asymptotically achieve the same performance as the maximum likelihood estimator even in the colored noise case. We first rewrite equation (1) as

with (^.)T and

denoting the transpose and the Hadamard (element wise) product, respectively. Further, w_N is defined similarly to y_N. Using (3), the NLS estimate can be obtained as (see, e.g.,P. Stoica and R. Moses, Introduction to Spectral Analysis, Prentice Hall, Upper Saddle River, NJ. 1997)

where ||.||₂ denotes the 2-norm, yielding the least-squares estimate ofp as

where (•)* denotes the conjugate transpose, which inserted into (10) yields the maximization

where

nimizing (12) for the general case of d unknown damping constants results in a (d+1)dimensional search, over and τ and β, each requiring about O(Nd²) operations to compute. However, exploiting the fact that one often can approximate βk~ β₀, the maximization in (12) can be obtained by a 2-D search, over temperature and the common damping constant βo; initial estimates for both these parameters exist, and only a quite limited search region is required. The detection variable is thus selected as

where II_Dτ,β0 is defined as in (13), but with β = β₀. We note that the matrix inversion in (13) may be poorly conditioned for specific T due to the resulting closely spaced frequency components. To alleviate this problem, we employ a low rank approximation technique, noting that a least squares solution can be found using the singular value decomposition. Let

where Σ is a diagonal matrix containing the d singular values of Q on the diagonal, and where U and V are unitary matrices. Further, let o_/ denote the /th singular value of Q, and note that the solution minimizing ||Qx - D*_{T β}y_N||2 can be found in Golub and Van Loan, Matrix Computations, 3^rd edition, John Hopkins.

where U_/ and V_/ denote the /th column of U and V, respectively, and. where d is the rank of Q, or alternatively the selected low-rank approximation of Q. Using (16), (14) can be expressed as

If, the temperature region of interest includes very closely spaced, or even overlapping, spectral lines, one should select d- d - \, otherwise d = d. ln the numerical simulations below, we have used the former.

We now examine the performance of the technique using both simulated and real NQR data obtained by the NQR group at King's College London. A typical NQR measurement allows for the estimation of a number of consecutive decaying echoes following every transmitted RF pulse in a multi-pulse sequence (MPS), with each echo consisting of 256 data samples. The damping of the NQR signal depends on the explosive examined; typically RDX decays much more rapidly than TNT, but its signals are much stronger and it is therefore easier to detect. Here, the examined real NQR data is obtained as the NQR response from a TNT sample at T = 298 K; the simulated data has been generated to mimic such a signal, and is from (1), with a damping factor of βo = 0.015 and with (uniform) random initial phases, and an autoregressive noise model derived from real NQR data. Typically, current techniques only measure the response of a single a priori known resonance frequency, se e.g. Y. Tan, S. L. Tantum and L. M. Collins, "Cramer-Rao lower bound for Estimation Quadrupole resonance signals in non-Gaussian noise", IEEE Signal Processing letters, vol. 11, no. 5, pp. 490-493, May 2004; to ensure the most beneficial performance for this approach, we will herein allow it to have perfect knowledge of the sample temperature, so that the most dominating resonance frequency is exactly known. We denote this the demodulation approach with perfect temperature knowledge (DMA-p). Typically, it is difficult to estimate the sample temperature with more than 5 degrees (K) accuracy; as a comparison, we therefore also include the estimate for a sample with 5 degrees offset, terming this the demodulation approach with realistic temperature knowledge (DMA-r). Figure 5 illustrates the performance gain of the detectors as the gain factor between the detection thresholds for a sample containing TNT and for one without TNT, as a function of the measured echo number. The methods have been evaluated on a single echo of real NQR data; the figure illustrates how the gain decreases for the higher order (more damped) echoes. A gain larger than one indicates a probable detection of the TNT sample. Here, exact temperature and damping knowledge has been assumed for both the DMA-p and the NLS methods; as the temperature of the reflecting sample will typically not be exactly known, the gain of the DMA-r is given as a reference. We again stress that the NLS method will not suffer from such temperature uncertainties, and the gain for the NLS method evaluated with a temperature uncertainty will yield similar gain as if the temperature was exactly

known. Figure 6 shows the gain factor for simulated data, as a function of the signal to noise ratio

(SNR), defined as

where σω denotes the standard deviation of ω(/)- Both these figures indicates a strong gain of the NLS method as compared to the DMA method, even when the latter is allowed to have perfect knowledge of the sample temperature. Figure 7 illustrates the receiver-operator curve, showing the probability of a correct detection, as a function of the probability of false alarm, for the 8th echo of a real NQR signal.

We now consider the Approximate Maximum Likelihood (AML) detector.

It should be noted that in the preceding sections, equations have been numbered from (1) to (18). In the following sections, equations are numbered from (3) to (83). References to particular equations are clear in context.

Detector based upon the free induction decay model using an Approximate Maximum Likelihood Estimation technique

We now consider the Approximate Maximum Likelihood (AML) detector. It is described here in mathematical terms but the skilled person will readily be capable of implementing this in data processing hardware and/or software, for instance on the control computer in the preferred embodiment. Let

where w_N is defined similar to y_N, (.)^T denotes the transpose, and

where

and diag{κ} denotes a diagonal matrix with κ along the diagonal. Let

where β denotes a vecto formed from the d unknown damping constants, β_k, k = 1, ... , d. As is well known, the maximum likelihood estimator is found as (see, e.g., [19])

where ||P||² _w = P*W^-1P. Further, R_w denotes the noise covariance matrix

where E{.} and (.)* denote the expec ation and the conjugate transpose, respectively. As R_w is typically unknown, one is normally forced to use an estimate of R_w, say Ȓ_w in (11). Such an estimate can be formed in various ways; herein, we propose using a low-order approximative noise model estimated from, a collection of NQR signals as described in Appendix A. By applying a prewhitening filter, we form

where C_w denote the unique lower-triangular Cholasky factor of R_w ^-1. i.e.,

where Ȓ_w is formed using the ow-order autoregressive model derived in Appendix A. As filtering a (damped) sinusoidal signal through a linear filter will yield a scaled and phase shifted version of the signal, z_N can be expressed as

where e_N = C*_wW_N is a (zero-mean) white additive noise with variance σ_ε ². Further,

where H(λ_k ) denotes the ( autoregressive) prewhitening filter evaluated at the mode λ_k ; see Appendix A for further details on the autoregressive model. Using (15). the minimization in (11) can equivalently be written as

where || . ||_F denotes the Frobenius norm. Let Ψ = ρΦ and

Then, the unstructured least squares estimate of Ψ is found as

where (.)^† denotes the Moore-Penrose pseudoinverse yielding the estimated phase vector

where Zψ_k denotes th argument of the kth element of Ψ. By inserting ( 19 ) in ( 17), the least squares estimate of ρ can be found as

yielding, the maximization

where

We remark that the inversion of

in ( 18) may be poorly conditioned for specific r due to tho resulting closely spaced frequency components. However, as only the phase is estimated in ( 19 ), the possible adverse effects of the poor conditioning will only slightly affect the resulting ρ estimate. Furthermore, maximizing (21) for the general ease of d unknown damping constants results in a (d + 1)-dimensional search, over r and β. However, exploiting the fact that one often can approximate β_k ≈ β₀. the maximization in (21 ) can be obtained by a two-dimensional (2-D) search, over temperature and the common damping constant β₀; initial estimates for both these parameters exist. and only a quite, limited search region is required. Furthermore, similar to the discussion in [20], one can often simplify the 2-D search further by replacing, it with two 1-D searches, first searching over τ and then over β₀. Such an approach will dramatically lower the computational complexity of the estimation, as the maximization is typically evaluated over a fine grid. In the following, we will refer to the estimate obtained from this approximative maximization as AML-1, whereas the full 2-D search is termed AML-2. It is worth noting that one may easily include a refined search in AML-1 by performing a local search around the maximizing temperature value; further, such a procedure can be iterated to find the temperature over a very fine grid. In Section 5, we illustrate both the mentioned complexity gain, and the corresponding loss of performance for the simplified search approach. Using; either of these estimates. we form the test statistic, T(z_N), as the (approximative) generalized likelihood ratio test (GLRT) for an unknown noise model. i.e., [21]

Using ( 24), the signal component is deemed, present if and only if

and otherwise not, where γ is a predetermined threshold value reflecting the acceptable probability of false alarm (P_f) ; here,

where QF_{r p} denotes the complementary cumulative distribution function for a

F-distributϊon with r numerator degrees of freedom and p denominator degrees of freedom [21].

Figure 8 shows two graphs of the two dimensional search conducted using the AML technique, here approximated as two, one-dimensional searches

Detector based upon the free induction decay model using an frequency selective Approximate Maximum Likelihood Estimation (FSAML) technique

We will now consider the FSAML technique. It is described here in mathematical terms but the skilled person will readily be capable of implementing this in data processing hardware and/or software, for instance on the control computer in the preferred embodiment.11

As it is known that the signal of interest lies in an essentially known narrow band of frequencies, one can beneficially derive a frequency selective detector which only considers a band of frequencies. Using (15), the Fourier transformed (prewhitened) data, vector can be expressed as

Here, w_k = e^i2πk/ _N . Assuming that the d damped sinusoids belong to the set

where k₁, . . . , k_M are M g ven integers, (27) can be expressed as

where

and where Ě_M is defined similar to Ž_M . ParaJlelling (18), the unstructured least squares estimate of Ψ is found as

where

yielding an estimate

similar to (19). We note that from a computational point of view, one should exploit that the indices of

are geometric series: thus,

can be evaluated as

ehere

Similar to (20), the least squares estimate of ρ can be found as

where the estimated τ and β are obtained as

with

Similar to the AML estimator, it is possible to approximate β_k ≈ β₀ , reducing (36) to a 2-D search. Furthermore, it is here also possible to approximate the 2-D search by two 1-D searches; in the following, wu refer to the approximative approach as FSAML-1 , whereas the full 2-D search is termed FSAML-2. Using the estimated ρ, the test statistic is then formed using (24). The resulting FSAML detectors has the benefit of operating only on the frequencies of interest, while still fully exploiting the a priori knowledge.

Figure 9 shows the results for frequency selective AML over the parameters of temperature and damping constant.

We now compare the probability of an accurate detection for a given probability of false alarm. This is the standard way to evaluate detectors; a 45 degree line (diagonal from [0,0] to [100,100]) is equivalent to guessing. We compare the AML and FSAML techniques with a demodulation (DMA) detector, one that estimates a single dominant frequency and uses this as a detector. We use the term DMA-p for a detector based on perfect temperature knowledge, i.e., it is formed from the amplitude of the most dominant peak - without missing it even for very low SNR. We use the term DMA-r, i.e., the (realistic) DMA detector with a temperature estimate that is 5 degrees off the true one. DMA-r can therefore be expected to be the worst performing technique,

Figure 10 shows the detection gain of the four techniques. DMA-r follows the bottom of the graph and DMA-p is slightly better. At the top of the graph FSAML is just slightly poorer than AML. Figure 11 shows both a signal magnitude graph and the, probability graph for a signal to noise ratio (SNR) of -HdB. It can be seen that the DMA-r technique will, for this low SNR, be essentially as good as guessing. DMA-p gives somewhat of a detection gain as one can expect- it is after all using the true frequency of the most dominant peak. However, AML outperforms all of the other techniques with FSAML in second place.

If the signal to noise ratio improves slightly, as shown in Figure 12, to -9dB then the performance of AML and FSAML improves still further.

The preceding sections are based upon the FID data model. However, as mentioned above, a wait time of 5 T] is normally required to produce a fully relaxed system, i.e., after the acquisition of an FID, one must wait 57_/ before applying another excitation pulse. This is a limiting factor in the detection of TNT in particular as 5Ti can be up to 30 seconds. One way to reduce this limitation is to use multiple pulse sequences to generate a train of echoes.

One may also fonn an averaged echo by summing the echo signals from an entire, or part of, an echo train. The ETAML and FETAML algorithms will then simplify to detectors similar to the AML and FSAML algorithms, being formed on the echo instead of the FID.

The use of echoes and trains of echoes is now considered in detail, beginning with derivations of data models for echoes and trains of echoes, and then continuing with description of detectors based upon those data models.

In the following sections, similar mathematical expressions are used to those used in the preceding sections directed to the Free Induction Decay model and associated detectors. However, the symbols used to represent particular terms in the following sections are not necessarily. the same as those used to represent equivalent terms in the preceding sections. The meaning of particular symbols is clear in context.

Data model - echoes and echo trains

Echoes, produced by using multiple pulse sequences, are beneficial in the detection of explosives as they enable the NQR signal to be sustained for longer than an FID, before the wait time of 5Ti seconds has to be adhered to. There are several relaxation processes that contribute to the decay of the FID. One of these processes is caused by the inhomogeneous nature of the , , sample, which effectively means that for a given resonant line, there is not one but rather a range of resonant frequencies. As a result, after excitation the signals from the isochromats will become out of phase with each other and lose coherence, resulting in a loss of coherency of the transverse magnetisation. However, it is possible to use certain pulse sequences to refocus these signals and

5 bring them back into phase. This is the basis for the formation of echoes.

An echo can be produced by applying a second pulse, called a refocusing pulse, a time t_sp after the initial pulse, often with its phase shifted by 90 degrees (with respect to the first pulse). This has the effect of refocusing the de-phased signals; consequently, a time t_sp after the

10 refocusing pulse, the signals from the isochromats are all back in phase, corresponding to the peak of the echo. By applying successive refocusing pulses in a sequence known as a pulsed spin locking (PSL) sequence, an echo train is produced. The echo train cannot be sustained indefinitely as the nuclei dephase. The time constant for this process following two pulses is the spin-spin relaxation time, generally denoted T₂. However, in a PSL sequence, for example, an even longer

15 decaying echo train is often found with a time constant T_2e.

As an echo is a refocused FID, it can be assumed that firstly, an echo consists of a set of sinusoidal components with the same frequencies as those seen in the FID and secondly, that from the peak at the centre of an echo to the end of the echo, these sinusoids have the same damping

20 factors (related to T₂,_k ) as seen in the FID. Therefore, the second-half of each echo may be modelled as the FID. Furthermore, due to the way the echo is formed, it can be assumed that the way the signal expands from the beginning of the echo to the centre is related to the way it then decays in the second-half of the echo. If, for the moment, we ignore the loss in longitudinal magnetisation over an echo, then the expansion up to the peak of the echo can be viewed as due to 5 the recovery of the transverse magnetisation, and the following decay due to its subsequent loss. As a result, if it is assumed that there was no loss in longitudinal magnetisation, the echo may be modelled as being symmetric about the echo centre.

Certain embodiments take advantage of the fact that, as described in the preceding 30 paragraph, an FID model can be extended to model echo responses, if the loss in longitudinal magnetisation over an echo is ignored, and in those embodiments the detectors described above (detector based upon the free induction decay model using a non-linear least squares technique; detector based upon the free induction decay model using an Approximate Maximum Likelihood Estimation technique; and detector based upon the free induction decay model using an frequency selective Approximate Maximum Likelihood Estimation (FSAML) technique) based upon the free induction decay model are applied to echo response data.

Considering the modelling of echo responses in more detail, and going beyond the FID model, the region of expanding modes will be determined by the pulse sequence parameters, allowing for a precise modelling of the envelope of each echo. To also account for the loss of magnetisation over the whole echo train, each echo envelope should further be damped by the spin-echo decay time. To stress that the spin-echo decay time will be different for each resonant line, we denote it as T_2e,k(τ), also indicating its temperature dependent nature. We note that T_2e,k(τ) is also dependent on the echo spacing, 2t_sp. To simplify our notation, we use ηu(τ) to denote the observed echo train damping parameter for the kth resonant line, noting that it is inversely proportional to

i.e.,

Figure 13 illustrates that an echo train is not a continuous function of time.

Between consecutive echoes there is a delay, where no data is acquired, to allow for the RF refocusing pulse required to generate the next echo. Here, we shall denote the number of samples between the first sample of one echo and the first sample of the next echo by μ.

In summary, we propose that the m^th echo of an echo train may be well modelled as,

where t = t₀, .... t_N-1 is the echo sampling time, measured with respect to the centre of the refocusing pulse, not necessarily being consecutive instances, but typically 'starting at t₀ ≠ 0 to allow for the dead time between the pulse and the first measured sample (after the pulse). Accounting for the dead time is equivalent to accounting for a large contribution to the first order phase correction in the NQR spectrum. For simplicity we will hereafter assume a uniform sampling starting at to, but note that the detectors could easily be generalised to allow for nonuniform sampling. Note that the frequencies, initial phases and relative amplitudes of the d sinusoidal components can be assumed to be constant over each echo.

It should be noted that equation (40) uses similar notation to that used in equation (1) but that here we let α_k=pK_k , where p is the common scaling due to the SNR and K_k is the a priori known (complex) relative amplitude of the k^th resonance frequency.

Detectors based upon the echo train model are now described.

The ETAML Detector

We now consider the ETAML detector. It is described here in mathematical terms but the skilled person will readily be capable of implementing this in data processing hardware and/or software.

Using equation (40), the /«th echo of the echo train can be written as

where w_N is defined similar to y_N, t denotes the echo sampling time, t = t₀, ..., t_N_₁ with N denoting the echo length, and

with

_; and

and

denoting the transpose and the Schur-Hadamard (elementwise) product, respectively. Using equation (41), an echo train consisting of M echoes can be written as,

where the operation vec[X] stacks the columns of the matrix X on top of each other. Furthermore,

The maximum likelihood estimate of θ can be found as (see, for instance, P.Stoica and R. Moses, Spectral Analysis of Signals. Upper Saddle River N.J.:, Prentice Hall, 2005)

where

and

with β denoting a vector formed from the d unknown sinusoidal damping components, β_k. Furthermore, R_w denotes the noise covariance matrix

where E{.} and (.)* denote the expectation and the conjugate transpose, respectively. As

R_w is typically unknown, one is normally forced to use an estimate of R_w in equation (48). Herein, we will form such an estimate using an approximate low-order noise model derived from the measured data, allowing the additive noise to be approximately described using a sixth-order AR model. One may therefore form a pre-whitened data model such that

where the

is a zero-mean complex white additive noise with

It should be stressed that since the echoes are separate acquisitions, they are filtered

separately. Furthermore,

is formed from the

last rows of

and

where

with [ .]_k denoting the kth index, and

where t_n t₀ + n and and

. Furthermore,

denotes the AR prewhitening filter, defined as

We stress that due to the required n-tap prewhitening filter,

will only contain

samples. Using equation (51), the minimization in equation (48) can be written as

where ||.||_F denotes the Frobenius norm. Thus, the least squares estimate of ρ can be found as

where

denotes the Moore-Penrose pseudo-inverse, yielding the maximization

Using the (τ,β)-pair maximising, we proceed to form the test statistic,

, as the (approximative) generalized likelihood ratio test (GLRT) for an unknown noise model (see for instance S.M.Kay, Fundamentals of Statistical Signal Processing, Volume II: Detection Theory. Englewood Cliffs, NJ., Prentice-Hall, 1998).

Using equation (62), the signal component is deemed present if and only if

and otherwise not, where γ is a predetermined threshold value reflecting the acceptable probability of false alarm (pj); here,

where

denotes the complementary cumulative distribution function for a F- distribution with r numerator degrees of freedom and^> denominator degrees of freedom.

The FETAML Detector

We now consider the FETAML detector. It is described here in mathematical terms but the skilled person will readily be capable of implementing this in data processing hardware and/or software.

In this section, the above presented ETAML detector is extended to formulate a

Frequency selective ETAML (FETAML) detector. As the temperature of the sample can be assumed to lie in a known temperature range, we may, using equation (2) (with the factor of 2π omitted) determine the range of frequencies each sinusoidal component may be present in. Hence, a frequency selective detector that only considers these frequencies can be derived. ^'

Consider the frequency regions formed by

with k₁, ...k_L being L given, not necessarily consecutive, integers selected such that equation (65) only consists of the possible frequency grid points for each of the d signal components; each such region is given by the minimal and maximal frequency values for that component considering the measured temperature and the size of the expected temperature uncertainty region. Denoting the measured temperature

, and the temperature uncertainty region the minimal and maximal frequency values for each component can be determined using

equation (2) with

and

, respectively. It should be stressed that each echo should be Fourier transformed individually to avoid discontinuities in the data. The Fourier transformed (prewhitened) data vector for the mth echo and Mi frequency bin can be expressed as

(66)

where

represents the kth frequency bin of the prewhitened noise

sequence associated with the mth echo,

, and

with

. Thus, over the (possibly overlapping) frequency regions of interest, equation (66) can be expressed as

and where

is defined similar to

. We note that the pre-whitening can alternatively be performed in the frequency domain. Using equation (68) the data model for the whole echo train can be expressed as,

where E_LM is defined similar to Z_LM. Using equation (71), the minimization in equation (48) can be approximated as

Similar to equation (58), the least squares estimate of p can be found as

where the estimated rand β are obtained as

We note that from a computational point of view, one should exploit that the indices of

form geometric series; the gth index of

can be written as,

with [ x ] denoting the integer part of x, and

We remark that several simplifications can be made to the above presented detectors. Firstly, for the general case of d unknown sinusoidal damping constants, β_k, the maximisations in equation (59) and equation (74) result in (d+l)-dimensional searches, over r and β. However, exploiting the fact that one can often approximate β_k « βo, the maximisations can be reduced to two-dimensional (2-D) searches, over temperature and the common damping constant βo; initial estimates for these parameters exist and only a limited search region is required. Furthermore, it is possible to approximate the 2-D searches by two 1-D searches, first over rand then over βo. We refer to these approximative approaches as ETAML-I (or simply ETAML) and FETAML-I (or FETAML). The full 2-D searches are termed ETAML-2 and FETAML-2. It is worth noting that one may easily include refined searches in ETAML-I and FETAML-I by performing a local search around the maximising parameter values; such a procedure can be iterated to accurately find the parameters on a very fine grid.

Secondly, we note that T_2e,ic($ is dependent upon several experimental parameters, and that the echo damping shifting functions of temperature for all possible echo spacings, excitation bandwidths and frequencies may not be available. As an alternative, to avoid determining the echo damping shifting functions for a particular experimental setting, one may instead treat the echo damping parameters as unknown constants, each denoted' %, searching for the values that best fit the measured data. Assuming that the sinusoidal damping parameters have been approximated to be the same, this would lead to a (d+2) dimensional search, over the common sinusoidal damping constant (βo), the d unknown echo damping parameters and temperature. We denote the approach, where the d+2 dimensional search is approximated by d+2 1-D searches, the ETAML-s detector.

Thirdly, we may approximate all the echo damping constants to be the same. In this case, it may be that an overall echo damping shifting function is available, in which case the search is still a 2-D search over temperature and βo, and this detector is simply a special case of the ETAML detector. However, if the common echo damping shifting function is unknown, then we treat the echo damping parameter as a single unknown constant, where η_k (τ) ≈ ηo. Here, the search is 3-D over the common echo damping constant, ηo, the common sinusoidal damping constant, β₀ and the temperature, τ. We denote the approach, where the 3-D search is approximated by three 1-D searches, the ETAML-a detector. We note that similar approaches can be taken for the FETAML detector. Table 2 summarises the approaches.

Table 2: The different ETAML and FKCAML detectors.

Examples using simulated and real NQR data

In this section, we examine the performance of the proposed detectors using both simulated and real NQR data measured at King's College London. Two sets of real data were used in the evaluation of the algorithms.

The first data set consisted of 200 data files, 100 with TNT and 100 without, each containing 20 summed echo trains. The data was collected in a shielded environment. The SNR of this data was too high to fully evaluate the algorithms on, so another set was obtained, where each data file consisted of only 4 summed echo trains. The simulated data, designed to mimic the first data set, was generated using equation (40), Table 1 and a sixth-order AR noise model with AR coefficients, C₀= 1, C₁ = -1.24, C₂= 1.28, C₃ = -0.91, C₄= 0.54, C₅ = 0.30, C₆ = -0.16.

The detectors were compared using simulated data with and without RFI, representing the cases of unshielded and shielded data, respectively. The RFI is modelled using a simplistic model consisting of discrete sinusoids with random frequencies, uniformly distributed over the interval [-π,π], and normally distributed amplitudes. Typically, current techniques only measure the response of a single a priori known resonance frequency; to ensure the most beneficial performance for this approach, we will herein allow it to have perfect knowledge of the sample temperature, so that the most dominating resonance frequency is exactly known. We denote this the demodulation approach with perfect temperature knowledge (DMA-p). In general, it is difficult to estimate the sample temperature with more than 5 degrees (K) accuracy; as a comparison, we therefore also include the estimate for a sample with 5 degrees offset, terming this the demodulation approach with realistic temperature knowledge (DMA-r). Both of these approaches are based on evaluating the response of a single resonant frequency. Alternatively, one may allow for a search to find the most dominant peak within the relevant frequency region; such an approach will improve the robustness to errors in the temperature measurement. Hereafter, we denote this approach the DMA-s detector. The AML, FSAML and DMA approaches are generally applied to echoes, or echo trains that have been pre-processed to produce a single echo with stronger SNR, whilst ETAML and FETAML detectors are applied to unprocessed echo trains. Hereafter, we will apply the AML, FSAML and DMA approaches to the summed echo train, formed by adding all the M consecutive echoes, while the ETAML and FETAML are formed on the full echo train signals allowing for the fine structure between the different echoes. We stress that as the former estimators operate on a single echo, they cannot be formed on the full echo train.

In the described examples, the AML-based detectors and the DMA-s detector use a search region over temperature of 290K to 310K (in 100 steps). Furthermore, the AML- based detectors use a search region of β_o= 0 to 0.03 (in 100 steps). The FETAML-s and ETAML- s detectors search for each of the d echo train damping parameters, % = 0.0002 to 0.0004 (in 100 steps). Finally, the FETAML-a and ETAML-a use a search region of η₀ = 0.0002 to 0.0004 (in 100 steps). Figure 14 illustrates the detection gain of the detectors, expressed as a ratio between the detection thresholds for a sample containing TNT and for one without TNT, as a function of file number for the first data set.

Similarly, Figure 15. shows the detection gain for simulated data as a function of SNR, here defined as

where σj and σ/ denote the variance of w(t) and y(t) - w(t), respectively.

These figures are formed using 100 Monte-Carlo simulations. Both figures show that the detection gains for the ETAML and FETAML detectors are significantly greater than those for the other detectors, and that the gain for ETAML is slightly greater than for FETAML. However, it should be noted that Figures 14 and 15 only illustrate the detection gain of the methods, and that the performance of the detectors will also depend strongly on the variance of this gain.

As detection is the problem of interest, we proceed to examine the receiver operator characteristic (ROC) curves for the detectors. Figure 16 shows the ROC curves of the discussed detectors, for simulated data without RFI, using 1500 Monte-Carlo simulations at an SNR of - 32dB. The figure clearly illustrates the improved performance of the proposed detectors over the older AML approaches and the DMA approaches.

Figure 17(a) shows a zoomed portion of the ROC curve, comparing the ETAML, ETAML-s and ETAML-a detectors, for simulated data using 3000 Monte-Carlo simulations at an SNR of -33dB. The figure shows that there is hardly any difference in performance when searching over the d echo damping parameters as compared to using the known temperature shifting functions (similar results for FETAML, FETAML-s and FETAML-a detectors). Further, there is only a small gain in assuming d different echo dampings as compared to assuming a single echo damping parameter.

In a realistic measurement environment, it is not possible to fully shield the sample and the antenna from RFI. As discussed above, there are techniques available for mitigating RFI, however, it is likely that significant amounts of residual interference will still remain. Here, we evaluate the performance of the detectors in the presence of RFI. Figure 18 shows the ROC curves of the detectors, for simulated data with RFI present, using 1500 Monte-Carlo simulations at an SNR' of -3OdB. The figure clearly shows that the FETAML detector is superior to the other detectors.

We proceed to examine the effect of zeropadding the data on the performance of the FETAML detector. When no zero padding is used, we will call the detector the FETAML-zO detector, when N zeros are added to each echo, we shall call the detector the FETAML-zl detector, and when 3 N zeros are added to each echo, we shall call the detector the FETAML-z3 detector. Figure 17(b) shows the ROC curves of the FETAML-zO, FETAML-zl and FETAML-z3 detectors, for simulated data without RFI using 3000 Monte-Carlo simulations, at an SNR of -32 dB. The figure shows, that there are negligible differences in performance when using different degrees of zeropadding.

Next, we compare the computational complexity of the proposed algorithms to the older AML approaches. The discussion above indicates there is no significant gain in using zero padding with the FETAML detector. Table 3 shows the average execution time over 100 executions, using the earlier specified search spaces, normalised with respect to the AML algorithm.

Table 3

Finally, we evaluate the proposed detectors on real data. As mentioned earlier, the SNR of the first data set was too high to compare the algorithms on, so a second set was obtained, the collection of which is now described. The sample, which consisted of creamed monoclinic TNT and weighed 18Og, was placed inside a solenoidal coil. The Q-factor of the coil was set to 60 in order to ensure the bandwidth of the probe was sufficient to excite the four-line region of TNT using a single excitation. Only the coil was placed in a, shield, so the real NQR data contained some RFI components. The second data set consisted of 1000 data files, 500 with TNT and 500 without, each taking around one minute to acquire and consisting of four echo trains summed up and phase cycled to reduce baseline offset. The echo trains were generated using a PSL sequence and were made up of 32 echoes, each echo consisting of 256 samples. The first echo of the echo train was discarded before being input to the algorithms as it is significantly distorted by contributions from the FID produced by the preparation pulse in the PSL sequence. The excitation frequency was 841.5kHz and the temperature of the sample 30 IK. In order to reduce differences in gain between the set containing TNT and the set not containing TNT, each ' data file was normalised before being input to the algorithms. The echo. damping shifting functions were not available for the second data set, hence the (F)ETAML-s and (F)ETAML-a implementations of the algorithms were used. The AR coefficients, Co = 1, ci =0.03, C₂ = -0.06, c₃ = -0.06, C₄ = 0.24, C₅ =0.25, C₆ = -0.44 of the prewhitening filter, defined in equation (56), were derived from the noise data. The SNR of the second data set is clearly higher than the SNR of the simulated data, used in Figures 16 and 18, yet Figure 19 shows the beneficial performance of the proposed detectors and Figure 20 confirms that there is hardly any loss in performance when assuming all the echo damping factors to be constant. Most surprising is the poor performance of the FSAML detector, which appears to suffer from the averaging of the echoes. The differences in performance between the AML-based detectors and the DMA-based detectors are more pronounced in the ROC curves produced using simulated data (Figures 16 and 18) as compared to those produced with the real data (Figures 19, 20 and 21). This as the simulated data is based on the first data set which has the relative amplitude scalings given in Table 1, whilst the second data set has the following relative scalings: K₁ = 0.14, K₂ = 0.71, K₃ = 1 and K₄ = 0.2, showing a more concentrated energy distribution; such a concentration will benefit the DMA approaches which only use the most dominating component. The difference in scalings is due to the longer repeat time used in the first data set (30 and 15 seconds, respectively). In the second data set, a shorter repeat time was chosen to enable the collection of more data in a shorter time. The relative scalings show that in the second data set, the TNT signal is more concentrated around the excitation frequency and so the difference between DMA and AML models become smaller. An additional reason for the reduced difference is that the SNR of the second data set is much higher than the SNR used in the simulated data, and so all the algorithms are working well, making it more difficult to see the differences in performance using the ROC curve.

A detailed description of the NQR data model has been provided and detectors have been described that exploit the data model of an entire echo train and/or exploit the temperature dependencies of the NQR signals. The detectors ensure accurate detection even in the typical case where the temperature of the sample is unknown. Numerical evaluation using both real and simulated data show a significantly increased probability of detection, for a given probability of false alarm, for the presented ETAML and FETAML detectors over AML-based and demodulation approaches.

It will be understood that the present invention has been described above purely by way of example, and modification of detail can be made within the scope of the invention.

Each feature disclosed in the description, and (where appropriate) the claims and drawings may be provided independently or in any appropriate combination.

Appendix A

Using NQR noise data measured by the NQR group at King's College London, an approximate low-order noise model was estimated from 98 recorded NQR noise sequences, each, containing 8192 samples. The noise data was found to be well described as an autoregressive (AR) process⁵; this conclusion was also indicated independently in [ 10]. Estimating AR-models of varying model orders, and examining the percentage of unexplained output variance (PUV ) |22], as shown in Figure 4, indicate that the model order should bo chosen at least as b = 4. and possibly as b = 6. Further, examining the autocorrelation function.

Figure 4: The percentage of unexplained output variance for AR (6) models as a function of b. r _ε(τ). for the residuals (prediction errors). ε(t). as shown in Figure 5. for b = 4 aud b = 6 with 99% confidence intervals, indicates that b = 4 might be inadequate to describe the properties of w(t), whereas selecting b = 6 yields essentially white residuals. Estimating the linear prediction parameters, ĉ_k, for k = 1 ,...,6. for each data set. yields the parameters given in Table 2, with the corresponding standard deviations. The residual variance, estimatce d from the sum of the squared residuals. was

= 542.

Figure 5: The autocorrelation function r_ε(τ) for the residuals ε(t) for b = 4 (dashed) and b = 6 (solid) with 99% confidence intervals.

Table 2: The mean values and standard deviations for the 98 estimates of ĉ₁, . .. , ĉ₆ of tbe AR(6) parameters.

Claims

1. A method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by combining a plurality of resonance parameters preferably as a function of a variable environmental parameter.

2. A method according to Claim 1, wherein the step of analysing the response signal comprises selecting a value of the variable environmental parameter, 'preferably in dependence upon the response signal.

3. A method according to Claim 1 or 2, wherein the step of analysing the response signal comprises, for at least one of the resonance parameters, selecting a value of that resonance parameter, preferably in dependence upon the response signal.

4. A method according to any preceding claim, wherein the step of analysing the response signal comprises using a model of the response signal, the model combining the plurality of resonance parameters, and preferably at least one of the resonance parameters being a function of a variable environmental parameter.

5. A method according to Claim 4, wherein the step of using the model comprises fitting the model to the response signal, preferably by selecting a value of the variable environmental parameter and/or a selecting a value of at least one of the resonance parameters.

6. A method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by fitting a model of the response signal to the response signal, the model comprising a plurality of resonance parameters each being a function of a variable environmental parameter, and preferably the step of fitting the model of the response signal to the response signal comprises selecting a value of each of the resonance parameters and selecting a value of the variable environmental parameter.

7. A method according to any of Claims 2 to 6, wherein the step of selecting a value of the at least one resonance parameter and/or the value of the variable environmental parameter is carried out in dependence on a pre-determined relationship between the resonance , , parameter and the variable environmental parameter and/or a pre-determined relationship between the plurality of resonance parameters.

8. A method according to Claim 7, wherein the pre-determined relationship is a

5 relationship between resonance frequency and temperature, and is preferably a linear relationship.

' 9. A method according to Claim 7, wherein the pre-determined relationship is a relationship between decay constant and temperature.

10 10. A method according to any preceding claim, comprising selecting the value of at least one further parameter.

11. A method according to any of Claims 2 to 10 as dependent on Claim 4 or 6, comprising selecting the value of the variable environmental parameter and/or the value of the or

15 each resonance parameter and/or the value of the or each further parameter in dependence upon how well the response signal fits the model.

12. A method according to any of Claims 2 to 11, wherein the step of selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter

20 and/or the value of the or each further parameter is carried out using an iterative procedure.

13. A method according to Claim 11 as dependent on Claim 4 or 6, wherein the iterative procedure comprises varying the value of the variable environmental parameter and/or the value of the or each resonance parameters and/or the value of the or each further parameter

25 and comparing the model to the response signal until the model is determined to match the response signal and/or until a pre-determined number of iterations have been performed.

14. A method according to Claim 13, wherein the step of comparing the model to the response signal comprises generating a measure of how well the model matches the response

30 signal.

15. A method according to Claim 14, wherein the model is determined to match the response signal when the measure of how well the model matches the response signal is within a pre-determined limit.

35

16. A method according to any of Claims 2 to 15, comprising selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter using a least squares technique, preferably a non-linear least squares technique.

17. A method according to any of Claims 2 to 16, comprising selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter using a Maximum Likelihood Estimation technique.

18. A method according to Claim 17, wherein the Maximum Likelihood Estimation technique is a frequency selective Maximum Likelihood Estimation technique or a time selective Maximum Likelihood Estimation technique.

19. A method according to any of Claims 2 to 18, wherein the step of selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter comprises selecting the value of one of the parameters, fixing the value of that parameter at the selected value and then selecting the value of another of the parameters.

20. A method according to any of Claims 2 to 18, wherein the step of selecting the value of the variable environmental parameter and/or the value of the or each resonance parameter and/or the value of the or each further parameter comprises selecting the value of at least two of the parameters simultaneously.

21. A method according to any of Claims 2 to 20, wherein the method comprises selecting, for at least one of the resonance parameters and/or for the environmental parameter and/or for the at least one further parameter, a range or set of values and selecting the value of that parameter from that range or set of values.

22. A method according to any preceding claim, comprising selecting a portion of the response signal and performing the analysis step independently of any analysis of any other portion of the or a response signal, and preferably performing the analysis step only on the selected portion of the response signal.

23. A method according to Claim 22, wherein the step of selecting a portion of the response signal comprises selecting the portion of the response signal being within a selected range or set of frequencies.

24. A method according to Claim 22, wherein the step of selecting a portion of the response signal comprises selecting the portion of the response signal being within a selected range or set of times.

25. A method according to any of Claims 4 to 24 as dependent on Claim 4 or 6, wherein the step of using the model comprises selecting the components of the model.

26. A method according to any of Claims 4 to 25 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing a free induction decay.

27. A method according to any of Claims 4 to 26 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing an echo.

28. A method according to any of Claims 4 to 27 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing an echo decay.

29. A method according to any of Claims 4 to 28 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing a train of echoes.

30. A method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by identifying a component of the response signal representing a train of echoes.

31. A method according to Claim 29 or 30, wherein the component representing the train of echoes comprises a component representing the decay of the train of echoes.

32. A method according to Claim 31, wherein the component representing the decay of the train of echoes comprises a component representing the decay of peak echo amplitude.

33. A method according to any of Claims 30 to 32, wherein the model comprises at least one component representing a steady state signal associated with a train of echoes.

34. A method according to any of Claims 4 to 33 as dependent on Claim 4 or 6, wherein the model comprises a component which represents the response signal as at least one decaying sinusoid, and preferably the or each decaying sinusoid corresponds to a respective resonance.

35. A method according to Claim 34, wherein at least one of the resonance parameters comprises a decay constant of the or at least one of the decaying sinusoids.

36. A method according to any of Claims 4 to 35 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing an undesired signal.

37. A method according to any of Claims 4 to 36 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing radio-frequency interference.

38. A method according to any' of Claims 4 to 37 as dependent on Claim 4 or 6, wherein the model comprises at least one component representing a noise signal.

39. A method according to Claim 38, wherein the noise signal comprises a non- white noise signal.

40. A method of testing comprising irradiating a sample, receiving a response signal and analysing the response signal by identifying a component of the response signal representing a resonance response and by identifying a component of the response signal representing a non- white noise signal.

41. A method according to any of Claims 38 to 40, wherein the method comprises determining at least one characteristic of the apparatus used to perform the method and selecting the component representing the noise signal in dependence upon the characteristic.

42. A method according to Claim 41, wherein the method comprises performing at least one measurement using the apparatus and determining the at least one characteristic in dependence upon the at least one measurement.

43. A method according to any preceding claim, comparing the response signal to a threshold and preferably generating an alarm signal in dependence upon the comparison.

44. A method according to Claim 43, wherein the step of comparing the response signal to a threshold comprises generating an output in dependence upon the response signal and comparing the output to the threshold.

45. A method according to any of Claims 4 to 44 as dependent on Claim 4 or 6, comprising comparing the model to a threshold and preferably generating an alarm signal in dependence upon the comparison.

46. A method according to Claim 45, wherein the step of comparing the model to a threshold comprises generating an output in dependence upon the model and comparing the output to the threshold.

47. A method according to Claim 44 or 46, wherein the output comprises a test statistic.

48. A method according to Claim 47, wherein the test statistic represents the likelihood of the response signal including a resonance response signal.

49. A method according to Claim 48, comprising generating the test statistic according to a generalised likelihood ratio test.

50. A method according to Claim 45, wherein the step of comparing the model to a threshold comprises comparing a component of the model representing a resonance response to the threshold.

51. A method according to Claim 45, wherein the step of comparing the model to a threshold comprises comparing the value of at least one of the resonance parameters to a threshold and preferably generating an alarm signal in dependence upon the comparison.

52. A method according to any preceding claim, wherein the variable environmental parameter is one of temperature, pressure and magnetic field.

53. A method according to any preceding claim, wherein the response signal is a time-dependent signal.

54. A method according to any preceding claim, wherein the response signal comprises a radio-frequency response signal and/or the step of irradiating the sample comprises applying radio-frequency excitation to the sample.

55. A method according to Claim 54, wherein the excitation comprises pulsed excitation and preferably the excitation comprises a sequence of pulses.

56. A method according to any preceding claim, wherein the response signal comprises a resonance response signal, the resonance response signal being one of a nuclear quadrupole resonance (NQR) response signal, a nuclear magnetic resonance (NMR) response signal or an electron spin resonance (ESR) response signal.

57. A method according to any preceding claim, wherein the resonance parameters comprise at least one of frequency and relaxation time.

58. A method according to Claim 57, wherein the resonance parameters comprise at least one of spin-lattice relaxation time and spin-spin relaxation time.

59. A method according to any preceding claim, wherein the resonance parameters comprise a plurality of resonance frequencies.

60. A method according to Claim 59, wherein each of the resonance frequencies correspond to a respective resonance arising from the same substance.

61. A method according to any preceding claim, being a method for detecting the presence of a substance containing a given species of quadrupolar nucleus.

62. A method according to Claim 61, wherein the substance is an explosive or a narcotic and preferably the substance is TNT or RDX.

63. A method according to any preceding claim being for the detection of a buried or concealed sample, preferably being for the detection of a sample concealed in baggage.

64. A method according to any preceding claim, wherein the method is a method of analysis of the composition of the sample and preferably the sample comprises a pharmaceutical.

65. A method according to any preceding claim, further comprising outputting the value of one or more of the parameters and preferably using the outputted values in a further measurement or in a further analysis procedure.

66. Apparatus for testing a sample, comprising means for irradiating the sample, means for receiving a response signal and means for analysing the response signal by combining a plurality of resonance parameters preferably as a function of a variable environmental parameter.

67. Apparatus for testing a sample, comprising means for irradiating the sample, means for receiving a response signal and means for analysing the response signal by fitting a model of the response signal to the response signal, the model comprising a plurality of resonance parameters each being a function of a variable environmental parameter, wherein preferably the means for analysing the response signal is adapted to fit the model of the response signal to the response signal by selecting a value of each of the resonance parameters and selecting a value of the variable environmental parameter.

68. Apparatus for testing a sample, comprising means for irradiating the sample, means for receiving a response signal and means for analysing the response signal by identifying a component of the response signal representing a train of echoes.

69. Apparatus for testing a sample, comprising means for irradiating the sample, means for receiving a response signal and means for analysing the response signal by identifying a component of the response signal representing a resonance response and by identifying a component of the response signal representing a non- white noise signal.

70. A computer program or computer program product adapted to perform a method as claimed in any of Claims 1 to 65.

71. Apparatus substantially as described herein with reference to and/or as illustrated in any of the accompanying drawings.

72. A method substantially as described herein with reference to and/or as illustrated in any of the accompanying drawings.