Pulse Reconstruction

The present invention relates to pulse reconstruction, and associated methods and apparatus. It has particular utility to the reconstruction of a pulse at terahertz frequencies, but is generally useful for all types of pulse reconstruction.

Terahertz electromagnetic radiation occupies a region of the electromagnetic spectrum between infra-red and microwave frequencies. Typically "terahertz" radiation is considered to include radiation having frequencies between 3x10 ¹¹ Hz and 3x10 ¹² Hz, corresponding to wavelength ranges from around 1 mm to 100 μm.

Terahertz radiation has previously attracted some attention for various imaging applications. In contrast to X-rays, terahertz radiation is non- ionising and so is safe to use on the human body. The radiation penetrates clothing, plastics, food and ceramics but is blocked by metal and water. Its absorption by water makes it attractive for use in medical imaging, as it can yield information on the relative water content of tissues. However, it has utility in a wide variety of imaging applications, including for example, surveillance, manufacturing applications and process control, satellite and high altitude telecommunications; as well as for other scientific applications including spectroscopy, study of condensed matter, and millimetre/sub-millimetre wave astronomy.

A limiting factor in the applicability of terahertz radiation imaging is the fact that it is strongly absorbed by water vapour in the atmosphere. This greatly reduces the range of a terahertz imaging system, and makes signal processing difficult due to spurious echoes and added noise.

Another reason that has hitherto prevented the widespread adoption of terahertz imaging is the difficulty in generating terahertz radiation, which requires the use of complex sources. These sources are either photoconductive or electro-optical in nature. For example, an electro- optical source such as a semiconductor can be triggered by a fast-pulsed laser, usually pulsed at a frequency in the order of femtoseconds. However, in this type of system there is little control over the time dependence of the output pulse. The initial pulsed electric field irradiated by the carrier acceleration is followed by a collection of after-runners generated by several processes, including subsequent dampening oscillations of the excited carriers and reflections from within the source structure which produced delayed and dispersed repeats of the initial output. The generation of these after-runners is unavoidable.

Furthermore, propagation through the atmosphere results in absorption and delayed re-emission of energy at the resonant frequencies of water molecules suspended in the atmosphere. These effects are generated both by the original excited pulse and by the after-runners.

The detector also adds similar echoes to the signal, and furthermore, if the terahertz pulse is transmitted through a dispersive region such as a waveguide, the various frequency components suffer varying phase changes and delays which further compounds the complexity of the signal.

These various factors combine to give a detected signal that is of a much more complicated form than would be expected from the carrier motion alone, which can result in the loss of finer detail and produce unwanted image artefacts.

For an imaging system using time delay measurements, the signal returned by the sample under observation may contain components with several different delays, causing confusion between echoes in the source output and overlapping pulses from the sample. If the returned signal includes both weak and strong echoes, the weak echoes may be obscured by the effects of the source after-runners.

Accordingly there is a need for an improved method for pulse reconstruction, in particular for terahertz imaging applications.

According to a first aspect of the present invention there is provided a method of pulse reconstruction for gaining information about a sample, comprising the steps of: obtaining a source pulse emitted by a radiation source; choosing a reconstructed input pulse; calculating a transformation generating function that produces the reconstructed input pulse upon convolution with the source pulse; introducing a sample and obtaining a signal emerging from the sample, said signal being subject to the same distorting effects as the source pulse; obtaining a reconstructed signal by convolving said transformation generating function with the signal emerging from the sample; and wherein: the step of calculating the transformation generating function comprises the steps of transforming the reconstructed input pulse and source pulse into frequency space, dividing them to gain the frequency spectrum of the transformation generating function, and then performing an inverse transform to obtain the transformation generating function.

A "pulse" is understood to mean a representation of how a signal's amplitude varies with time.

Preferably, the step of obtaining the source pulse comprises detecting the signal emitted by the source in the absence of any sample.

Preferably, the reconstructed input pulse is localised in time.

Note that when we say "localised" this is with the knowledge that perfect localisation will most likely never be achieved. There will almost inevitably be some fore-runners and/or after-runners in the reconstructed input pulse, however, the requirement for localisation in time corresponds only to a requirement that these fore-runners and/or after-runners are not of significant amplitude when compared with the amplitudes of the pulse in the chosen time localised region.

Preferably, the length of the reconstructed input pulse is chosen between a maximum value corresponding to a shortest time interval between sample echoes, and a minimum value that yields a frequency spectrum that extends to a predetermined threshold beyond the frequency spectrum of the source pulse.

Preferably, the reconstructed input pulse is chosen to have a frequency spectrum with zero content outside a chosen cut-off frequency.

Preferably, the cut-off frequency is chosen to lie between a minimum value that corresponds to a chosen minimum acceptable resolution of the reconstructed signal, and a maximum value that corresponds to a chosen maximum acceptable level of noise in the reconstructed signal.

Preferably, the spectrum of the reconstructed input pulse approaches zero smoothly at the cut-off frequency.

Preferably, the reconstructed input pulse is calculated from binomially distributed amplitude coefficients.

Preferably, the reconstructed input pulse is generated from a seed function that comprises a collection of delta functions, which most preferably are chosen to have amplitudes related to the amplitude coefficients of the reconstructed input pulse.

Preferably, the frequency spectrum of the reconstructed input pulse comprises p real components and q imaginary components.

Preferably, p + q > 0.

Preferably, p is an integer between one and ten and q is set to zero or an integer between one and three.

Preferably, the frequency spectrum of the reconstructed input pulse has the form F(ω) = e ^!OK °[cos( — )] ^p [/sin( — )] ⁹ , or a mathematical equivalent

•J c -J c thereof.

Preferably, an exponential convergence factor is applied to one or more of the source pulse; the signal emerging from the sample; and the reconstructed signal: in order to remove spurious effects caused by sample echoes falling outside the time of the reconstructed input pulse.

Preferably, the source pulse and/or the signal emerging from the sample is multiplied by a decreasing exponential.

Preferably, the reconstructed signal is multiplied by an increasing exponential.

Preferably, an exponential factor is chosen so that the relevant pulses are reduced smoothly to zero at the end of the time period of the reconstructed input pulse.

Preferably, a filter is applied to remove anomalous values from the spectrum of the reconstructed signal.

Preferably, the filter compares a central value with the values of its nearest neighbours on opposite sides, and replaces the central value with the average of the neighbouring values if the product of the differences between the central value and a first neighbour and the central value and a second neighbour exceeds a predetermined threshold.

Preferably, the filter is applied more than once, so that in a first pass the major anomalous values are removed, and in a second pass and in any further passes, a finer filtering operation is performed.

Preferably, the predetermined threshold is varied between passes of the filter.

Preferably, the radiation source emits electromagnetic radiation at TeraHertz frequencies.

According to a second aspect of the invention, there is provided an imaging system comprising a radiation source and a radiation detector, means for receiving a sample to be imaged, and signal processing means arranged to carry out the method of pulse reconstruction according to the first aspect.

Preferably, the imaging system is a TeraHertz imaging system.

Preferably, the imaging system comprises a movable detector for imaging a fixed object from a plurality of different angles.

Preferably, the detector is in connection with a flexible optical fibre, with the signal processing means preferably arranged to compensate for signal artefacts introduced by the fibre.

Preferably, the imaging system is a tomography imaging system, which may or may not use TeraHertz radiation.

According to a third aspect of the invention, there is provided a computer apparatus arranged to carry out calculations needed for the performance of the pulse reconstruction of the first aspect. In a further aspect, there is provided a computer program product comprising code executable on a computer to produce the computer apparatus of the third aspect. The computer program product can be provided on a computer readable medium or made available for download, or provided as a signal.

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

Fig. 1 shows a pulse seed and reconstructed source pulse;

Fig. 2 shows a source pulse and pulses reflected from two PTFE slabs;

Fig. 3 shows recorded source and sample pulses;

Fig. 4 shows a reconstructed sample signal showing internal reflections in PTFE and simulated signal;

Fig. 5 shows measured and calculated delays of reconstructed reflections;

Fig. 6 shows recorded and simulated pulses for an example practical application;

Fig. 7 shows reconstructed signals for the application of Fig. 6;

Fig. 8 shows a plan schematic view of the operation of a back-projection technique; and

Fig. 9 shows a plan schematic view of an apparatus for performing the back projection technique illustrated in Fig. 8.

The pulse reconstruction algorithm described herein addresses the problems mentioned above by numerically processing the sample signal together with a reference signal which has been separately recorded, to remove the undesired effects of the source after-runners. The result is to show what the sample signal would be if the sample had received a single, simple pulse of a form that can be chosen appropriately. The reference signal is produced in such a way that it includes all the distorting effects, but without the influence of the sample. For example, in a reflection

measurement the reference signal could be obtained by replacing the sample with a mirror. These principles will be described in more detail.

Principle of Pulse Reconstruction

At the low power levels available from pulsed THz sources, the response of the sample is linearly proportional to the input electric field. The signal s(t) emerging from the sample can be expressed as a convolution of the electric field of the source pulse f(t) arriving at the sample with the impulse response function g(t) for the sample:

s(t) = g(t) ®f(t) (1 )

The impulse response function g(t) for a unit impulse at t=0 is in general a sum (or integral) of delayed impulses at times t _p (> 0) with amplitudes A _p

g(t) = ∑A _p δ(t -t _p ). (2)

P

The signal from the sample then has the form

s(t) = ∑A _p f(t -t _p ). (3)

P

The inventors have come up with a method of pulse reconstruction where a signal from a sample to be imaged is processed together with a separately recorded reference signal to remove the undesired effects of the source after-runners.

The first step is to generate a reference signal, which is produced in such a way that it includes all the distorting effects, but without the influence of

the sample to be imaged. In a reflection measurement the reference signal could for example be obtained by replacing the sample with a mirror, or in a transmission measurement the reference signal could for example be obtained by removing the sample so that the source pulse is transmitted through empty space (i.e. the atmosphere).

An ideal source f(t) would comprise a single, well defined pulse. However, as explained above, the actual source pulse has various after-runners, which themselves cause further undesired reflections within the source. These then result in a signal s(t) in which weaker reflections are indistinguishable. However, an idealised input pulse can be represented by an arbitrary input pulse fit) , and then the transformation generating function u(t) that acts on the measured source pulse f(t) can be calculated. The new reconstructed input pulse is then

Applying the same transformation u(t) to the sample signal s(t) produces the reconstructed signal

s(t) = <t) ®s(t) = ∑A _p f(t -t _p ). (5)

P

Comparing equations (3) and (5) shows that the reconstructed signal s is the response the sample would produce to the chosen input signal / .

Deconvolution of the Filtered Signal:

The relationship between the source pulse f(t) and the required transformation u(t) is given by equation 4 above. From this, the frequency

spectrum of u(t) can be derived from the frequency spectra of the measured source pulse f(t) and the new chosen input f(t) :

tf(ω) = -M (6)

The transformation generating function can therefore be calculated as the inverse transform of its frequency spectrum thus obtained.

The spectrum of the reconstructed sample pulse s(t) then follows from equation 5:

F(CO)

S(ω) = S(CO). F(co) U)

where the upper case symbols indicate Fourier Transforms.

It can be seen from the form of equations (6) and (7) that the calculation may become unstable if F(ω) takes small or zero values. This is a well- known difficulty with deconvolution and can arise in several ways:

For frequencies outside the effective range of the source, F(ω) is derived from random noise, but F(ω) can be devised to be zero at these frequencies, stabilizing the calculation.

In the presence of atmospheric water vapour, there is strong attenuation at the resonant frequencies of the molecules. F(ω) becomes small at these isolated frequencies, and consequently U(ω) may have anomalously large

or small values at isolated points. The effect of this on the reconstructed signal can be removed by the filter described below.

The reconstructed output pulse is given by the inverse Fourier Transform of ,S .

Choice of Input Pulse / (t):

The choice of the input pulse / (t) is arbitrary, but there are certain preferred embodiments where it is preferable that it meets certain criteria.

In order to make clear the observation of time separated return signals from the sample, the reconstructed input signal should be well localised in time, with no significant fore- or after-runners. The length of the input should be short enough to resolve the time intervals between sample echoes, but not so short that its frequency spectrum extends much beyond the range of the THz source. For example, for a source with a 3 THz maximum effective frequency, the reconstructed source pulse should have a length of ~ 0.3 ps or more, corresponding to values of f _c equal to 3 THz or less.

However, shorter input signals result in increasing noise in the reconstructed sample signal. The choice of the input signal / (t) and corresponding transformation u(t) is therefore a compromise between time resolution and noise filtering. This compromise can be optimised by choosing an input / (t) that has a frequency spectrum with zero content beyond a chosen cut-off frequency, f _c . This eliminates higher frequency noise components, avoiding the artefacts that would arise from dividing by random noise values beyond the range of the source. At the same time

the subsidiary fore- and after-runners should be as small as possible, to avoid hiding a weak signal in the vicinity of a strong one.

Existing solutions for the problems of the apodisation of diffraction effects, antenna side-lobes and the design of RF channel filters can also be implemented.

In one embodiment, a set of input functions can be based on the binomial distribution. These have sufficient flexibility to allow for variations in the spectra of different THz sources, some of which have low frequency as well as high frequency cut-offs.

The amplitudes of the subsidiary fore- and after-runners of the input pulse / (t) are related to the behaviour of the spectrum of / (t) as it approaches the cut-off frequency, f _c . If the spectrum goes discontinuously to zero, it will result in relatively large subsidiaries to / (t), with amplitudes diminishing slowly away from the main part of the pulse. It is therefore preferable to design / (t) in such a way that its spectrum and the first few

(i.e. one, two, three or more) derivatives of the spectrum go smoothly to zero at the chosen cut-off frequency f _c . This corresponds to an input pulse that has only a few small subsidiaries close to the main part of the pulse. The frequency spectrum of / (t) is taken as

F(ω) = e™°[cos(^)F[/sin A]* ,|ω| < 2π/ _c (8)

and

F(ω) = 0, ω > 2π/ _c (9)

where t ₀ is the time of arrival of the main part of the source pulse, and the cut-off frequency f _c should not be taken beyond the effective range of the source.

The cut-off frequency f _c controls the time resolution limit ~ 1/f _c . This must be related to the maximum frequency of response of the THz source and detector. If f _c is chosen too high, F(ω) in equation 7 will have small random values within the range where F(ω) is non-zero. This will result in strong high frequency noise in the reconstruction. Choosing small values of f _c gives good noise reduction at the expense of time resolution.

Non-negative integers p and q are chosen to suit the requirements of a particular source spectrum, optionally with (p + q) > 0. For most sources p can be taken to be a small integer with q = 0. The parameter p with value 2 is suitable for most types of signal, as this gives amplitudes of the subsidiary oscillations of the source less than 3% of the main peak. If it is necessary to detect a weak signal in the vicinity of a strong pulse, larger values of p will reduce the fore- and after-runners of the strong pulse. The interpretation of the reconstructed signal is most transparent when q = 0. This gives an input pulse / (t) with a single maximum, and width proportional to f _c ^~1 , and relative signs of the pulses from the sample are immediately apparent. For example, if the sample has an air gap between layers, the reflected signals from the two faces of the gap have opposite signs. This is clearer to interpret for an input pulse with a single positive peak. However, if the source is deficient in low frequency power, taking q = 0 may result in the appearance of excessive noise at the low frequency end of the spectrum of the reconstructed signal. Larger values of q give incident source pulses with (q + 1) half-cycles of alternating signs. Such pulses have reduced low frequencies and result in less low frequency

noise, but can cause confusion in the interpretation of close or overlapping sample pulses.

Numerically, the spectrum F(ω) can be conveniently generated by using the fast Fourier Transform of a simple 'seed' function. The 'seed' can take any form, depending on the chosen form of f(t) , but as an example, it can be of a general form comprising a collection of delta functions with amplitudes chosen so that the first zero of the frequency spectrum is at f _c .

This results in an input function /(t)with width ~ Mf ₀ .

As a particular example, the seed can be generated by convoluting together p factors of type A with q factors of type B, where

A = [δ (t)+δ (t ~)]/2 (10)

J c

B = [δ(t)-δ (t-y)]/2 (11 )

J C

The seed function is in this example a collection of delta functions with amplitudes equal to binomial coefficients and oscillating signs when q > 0. Cutting off the spectrum of the seed function at ± f _c and inverting the

Fourier Transform gives the form of the input pulse / (t). An example of seed and pulse forms is shown in figure 1. For all graphs in figure 1 , the y-axis is an amplitude scale in arbitrary units, while for the graphs in the time domain the x-axis is plotted in units of s ^"12 and for the graphs in the frequency domain the frequency is plotted in units of terahertz. Fig. 1a shows seed pulse 10 which corresponds to p=2 and q=0 as defined above. The Fourier Transform 12 of seed pulse 10 is shown in Fig. 1 b. A

cut-off frequency f _c is defined, in this case 1THz. It can be seen from Fig. 1 b that the frequency spectrum of the seed pulse goes smoothly to zero at the positive and negative values of the cut-off frequency f _c . the cut-off frequency f _c is then used to truncate the spectrum, as shown in Fig. 1c. An inverse Fourier transform is then applied to the truncated spectrum 14 to give the reconstructed pulse fit) (16), shown in Fig. 1d. It can be seen here that, as has been mentioned above, the frequency cut-off used in the construction of fit) means that the subsequent signal reconstruction does not contain noise artefacts that would arise from dividing by random noise values beyond the range of the source. That is, the choice of fit) means that the deconvolution can be performed without diverging.

Elimination of Artefacts from the Reconstruction

Once the reconstructed signal s(t) has been obtained, further steps can be carried out to remove artefacts that remain in the reconstruction. These artefacts have various causes.

Signal Truncation

Spurious effects can appear in the reconstruction if the length of the time scan is insufficient to record the full range of source and sample signals. For example, if the sample produces a large delay, the effects of an echo in the source signal may be missed in recording the sample signal.

Therefore, the length of the time scan needs to be chosen carefully to make sure it is long enough. However there will be cases where it is unavoidable that a significant part of the source or sample signal continues beyond the range of time that is recorded, which as well as the

loss of information can cause unwanted effects in the reconstruction. Now, the fast Fourier Transform routines used in processing the signals have the implicit assumption of periodicity beyond the recorded range, meaning that the effects of signal truncation appear as spurious signals which can precede the incident pulse. These unwanted effects can be accounted for by use of an exponential convergence factor applied to one or more of the recorded reference signal, recorded source signal and reconstructed signal.

For example, a recorded reference (or source) signal and sample signal can be multiplied by a decreasing exponential, exp(-bt) and the reconstructed signal can be multiplied by an increasing exponential exp(bt), thus minimising the effects of the truncation. Parameter b is chosen so that the recorded signals are reduced smoothly towards zero at the end of the time scan, and usually has a value of a small multiple of the inverse time scan length, for example two or three times the inverse time scan length.

Water Vapour Absorption

Water vapour in the atmosphere has several sharp absorptions at frequencies in the THz range, which can reduce both the source and sample signals to noise at these frequencies. From equation 7 it can be seen that the reconstruction 5 ^* (ω) will contain narrow bandwidth noise at the water vapour resonances, but other frequencies will not be affected. The spectrum of the reconstructed pulse will have isolated noise peaks at the water vapour resonances, with lower noise over the remainder of the frequency range. This type of noise can be eliminated by application of a suitable filter.

After removing the phase factor exp(iωt ₀ ) from the Fourier transform of the reconstructed sample pulse F(ω) , the oscillations of the real and imaginary parts are reduced, except at the water vapour resonances where there may be anomalously large positive or negative values. The effect of these noise peaks can be removed from the reconstruction by replacing a(n) by /4[a(n+1 ) + a(n-1 )] when

[a(n) - a(n - \)][a(n) - a(n + \)] > C (12)

where a(n) is the nth component of the real or imaginary parts of the Fourier Transform, and a cut-off level C is chosen so that only the anomalous values are altered and the remaining values are unchanged. This form of filter removes sharp peaks representing noise, without reducing the information content of the reconstructed signal.

The value of C can be selected by averaging [a(n) - a(n-1)] over the length of the data and operating the cut-off when the local value exceeds the average by a chosen factor (for example, fifteen). The results of this filtering can be assessed by examination of the spectrum of the reconstructed signal. If the spectrum has a single large anomaly and several smaller anomalous values, a first application of the filter will only remove the largest anomaly. A second application will then remove the remaining smaller anomalies.

Applications

To illustrate the application of the above principles, some exemplary applications will now be described.

Transmission through PTFE Slabs with Air Gap

Figure 2 shows two PTFE slabs 18, 20 of thicknesses of 2.39 and 2.45 mm, with a 1 mm air gap 22 between them. A THz source 24 irradiates a source pulse f(t) 26, and reflected pulses 28-34 are generated by each slab surface, and a THz detector 36 is arranged for detection of the reflected pulses 28-34.

Figure 3 shows incident source signal 38 and transmitted sample signal 40. Both signals are strongly affected by atmospheric water vapour, so that only the time delay of the main transmitted pulse can be seen. The internal reflections inside the slabs are hidden by the water vapour emission.

Figure 4 then shows the reconstructed pulses 42 and simulated pulses 44, with an offset section illustrated at ten times the scale of the main graph for clarity of illustration. The reconstructed pulses 42 show the direct transmission at 15.5 ps together with five distinct delayed pulses above the noise level. Four of these delayed pulses, indicated by the bars 46-52, can be identified with the calculated transmission times for refractive index 1.39 with incident pulse at 8.8 ps. The origin of the pulse at 27.5 ps cannot be identified. The four surfaces of the two slabs produce six secondary internal reflections, but the similarity of the two thicknesses produces two pairs of pulses too close to be resolved.

The experiment was repeated for a range of separations up to 5 mm, and figure 5 shows the detect pulse delays, in clear agreement with the calculated times, assuming a refractive index of 1.39. The amplitudes of the internally reflected pulses are approximately 3% of the incident amplitude, which is consistent with the Fresnel formula for this index.

Simulation of Dispersion Compensation in a Waveguide or Endoscope

An example practical application will now be described. For internal examinations using a THz signal transmitted via an endoscope, it is necessary to compensate for the effects of dispersion. This can be simulated by treating the endoscope as a waveguide with a cut-off wavelength of approximately twice the endoscope diameter. The reference signal for the reconstruction could be obtained by recording the signal returned through the endoscope with a simple mirror at the remote end. Replacing the mirror with the sample then provides the signal for reconstruction. Figure 6 shows a recorded source signal 60 including the effects of water vapour, a simulated returned reference signal 62 from a mirror, the simulated signal 64 from a layer of transparent material approximately 0.1 mm thick. Figure 7 shows the reconstructed reference signal 66 and sample signal 68, showing reflections from the front and back surfaces. The waveguide had a length of 500 mm and cut-off wavelength of 4 mm.

Tomography

In tomography applications at THz frequency, information is carried by signals that have been scattered by the object. The information can be collected in the time domain for detector positions at a range of angles around the object, or alternatively the detector can remain fixed as the object is rotated (this is easier to achieve for smaller objects). In either case the information is in the form of angle dependent time delays, with separate parts of the object giving different delays. Pulse reconstruction can be applied to the signals at this stage, if it is necessary to remove interfering effects from source echoes, water vapour, etc. The image of

the object is then obtained by mathematical back projection. A second stage of reconstruction may then be applied to the resulting image after back projection is carried out. This type of image reconstruction can in principle be applied to the images produced by any other instrument, as described below.

An image forming instrument responds to some physical property of a surface or slice of the object described by the function obj(r), dependent on the two-dimensional vector r. The response of the instrument is given by the instrumental point spread function (psf) ip(r), and the resulting first image fim(r), before image processing, is given by the convolution

fim(r) = obj(r) ®ip(r) (13)

This first image may contain artefacts such as haloes and false lines unrelated to the objects. Signal reconstruction can suppress these artefacts by the use of a reference image produced by scanning a phantom consisting of a single point object against a uniform background. For a phantom object smaller than the instrumental resolution, this can be represented by δ(r) and the reference image is:

ref{r) = δ (r) ® ip(r) = ip{r) (14)

The instrumental psf is thus determined, but this cannot be used alone to reverse the convolution in equation 13 as it results in divergence.

The reconstructed image rim(r) is required to be a convolution of the object with a new psf cp(r) that is chosen to comply with the conditions given below.

rim(r) = obj(r) ®cp(r) (15)

Transforming into two-dimensional k-space with kernel e ¹ - ¹ and using capitals to denote functions in k-space, equations 13 and 15 become

FIM ^' (K) = OBJ (k) ^■ IP{K) (16)

RIMiK) = OBJ (k) ^■ CP(K) (17)

Combining these to eliminate the unknown OBJ:

RIM(K) = FIM(K) ^■ CP(K) I IP(K) (18)

Equation 18 illustrates the potential, for divergent results. The range of IP is limited by the instrumental resolution. For k values beyond some limiting value kMAx IP(Js) reduces to zero or random noise. To avoid divergences, CP(J<) is chosen to be zero in regions of k-space where IP becomes small. The ratio CP/IP then remains finite over the whole of k- space.

Defining

T(K) = CP(K)IIP(K) (19)

the reconstructed image becomes

RIM(K) = FIM(K) ^■ T(K) (20)

Returning to spatial co-ordinates, the reconstructed image is given by a convolution of the first image with a reconstruction function t(r)

rim(r) = fim(r) ®t(r) (21 )

Choosing a suitable new psf cp(r) will remove artefacts if cp satisfies the conditions: (a) cp(r) is well localized in space with single maximum and no significant secondary peaks; (b) transform CP(k) vanishes outside the region limited by k _M Ax

The transformation defined by equation 21 can be described as a 'reconvolution', as the original instrumental psf has been removed and replaced by a more suitable psf chosen to comply with the instrumental characteristics.

For isotropic spatial resolution, it is convenient to define the new psf in terms of its transform in k-space. Suitable functions are given by

CP(K) = [cos(π \ k \ l2k _M )Y,\ k \< k

(22) CP(k) = 0,\ k \> k _M

with p a small positive integer. The spatial resolution is then ~ kMAx ^"1

The influence of random noise on the reconstructed image can be controlled in two ways. The choice of cut-off kMAx determines the effective bandwidth of the system. Reducing the chosen value of kMAx will reduce the spatial resolution in the final image, but also improve the signal-to- noise ratio.

A second effect can arise when the noise results in IP becoming small at isolated points in k-space inside the cut-off. The corresponding T(k) then has isolated large anomalous values. These can be numerically identified and replaced by a local average before transforming into the reconstruction function t.

Figures 8 and 9 illustrate the operation of mathematical back-projection for a tomography application having a fixed object and a detector that is rotated around the object to be imaged. Referring to Figure 8, a set up for the imaging of an object comprising two vertical rods 70,72 of metal or dielectric is shown.

The scattering from the objects 70,72 takes place in the horizontal plane and therefore information relating to the position of the scattering centres can be collected by moving the detector 74 over a horizontal arc 76 centred within the object area. The reconstructed signals then show pulse delays dependent on the coordinates of each scatterer and on the detector angle. The back-projection is then essentially a two dimensional problem.

For a scattering object, O with coordinates (x,y), irradiated by a THz beam with focus F at (0,R), and with the detector 74 at (Rsinθ, Rcosθ), the extra optical path length relative to an object at the origin is given by (OD - OF) where OD ² = (x - Rsinθ) ² + (y - Rcosθ) ² and OF ² = x ² + (y - R) ² . For each pixel in the field of view enclosing the sample and for each position of the detector the expected time delay can be calculated, and the corresponding value of the reconstructed signal for that time delay selected. Summing over the range of detector angles produces the back-projected image.

For small detector angles, where the detector is close to the incident beam focus, the time delays are relatively insensitive to the scatterer positions

and this part of the data does not lead to a well resolved image. As the angles are increased, the spatial resolution rapidly improves.

The spatial resolution is also dependent on the length of the THz pulse and bandwidth. If the detector is scanned over a range from +α to - α, and the bandwidth limited time resolution is t, the appropriate resolution limits can be set at a longitudinal resolution of δy ~ ct / (1 - cos α) and a transverse resolution of δx ~ ct / (sin α). These become equal for a ± 90° range of scanning.

As an example of this method, Figure 9 shows a plan view of an apparatus comprising a fibre-fed time domain spectroscopy (TDS) system. The output from a laser 80 is split in two, with 70% propagating through free space to generate the THz radiation, using a photoconductive switch, with the other 30% being used for detection. For THz detection, a photoconductive switch 88 is excited using the NIR probe pulse, which has propagated along a length of single mode optical fibre 82, which is three metres long in this example. Using the fibre 82 means that a stable detection unit 84 can be built that can be rotated about the target 86 of interest. The fibre 82 ensures that the probe laser beam is always focused to the correct point on the photoconductive receiver, a gap of area 250 μm ² . In both cases the photoconductive switches 88,90 are fabricated on low-temperature grown Gallium Arsenide. The pulse precompensator 92 partly compensates for the laser pulse dispersion in the optical fibre.

A major problem when sending short, intense light pulses along a length of optical fibre is the change in shape and lengthening of the pulse after propagation. This is because the intensity of the pulse causes non-linear effects within the fibre leading to both Group Velocity Dispersion (GVD)

and Self-Phase Modulation (SPM). GVD introduces a spectral dephasing, whilst SPM introduces a temporal de-phasing. These two effects within the fibre lead to a strong, linear positive chirp on the optical pulse.

SPM creates new frequencies "redder" than the central wavelength at the leading edge of the pulse and "bluer" frequencies than the central wavelength at the trailing edge. This increased spectrum breadth implies that the pulse can be reduced in duration. To compress the pulse the red wavelengths must be slowed down with respect to the blue ones; and this can be achieved using GVD. One way to do this is to use a diffraction grating which acts as a dispersive delay. Passing the pulse through the delay once spatially disperses the differing wavelength components, with a second pass then re-collimating the beam. The optical path length through the delay is then longer for longer wavelengths; and thus negative group velocity dispersion is induced on the pulse, which compensates for the positive chirp of the self-modulation.

As optical pulses become shorter, it becomes more difficult to compress them using the above technique. The limitation arises with the GVD: as more and more group velocity dispersion is required, so the size of the diffraction grating must grow, both in terms of the size of the grating and the optical path length required for each pass through the grating. However, the combination of GVD and SPM can be used to ensure that the femto-second pulses retain their shape and do not grow in length by too much on passing through a long optical fibre.

The present inventors have implemented this apparatus with pulses of average power 15 mW and duration of approximately 75 fs (assuming Gaussian shapes) launched into a 3m single mode optical fibre. The pulses were initially pre-chirped, using a Coherent Technologies mini-

GVD unit; and pulses of length 140 fs, and average power 6 mW, exit from the fibre. By securely mounting the optical fibre to the rotation stage, the pulses from it were focused down onto a bow-tie photoconductive detector. THz radiation from the emitter was focused onto the backside of the GaAs substrate. The current induced in the receiver was then measured using a lock-in referenced to the I kHz chopping frequency of the incident THz pulses. It was found that the stability of the detection unit could be much improved by inserting an amplifier circuit between the detector and the lock-in. By placing the amplifier as close as possible to the detector, the very low current output from the detector (of the order of a few picoAmps) could be readily amplified to a voltage of the order of micro-volts. This process amplifies the THz pulses, whilst having little effect on the background signal (due to the dark current in the receiver, laser fluctuations etc.) and thus improves the signal to noise ratio of the system as a whole. Pulses containing frequency components up to 1.5

THz, and with signal to noise ratio of a few hundred, were routinely output from this system. By rotating the detector, THz radiation scattered from various objects could be measured over the range of angles -90 to +140 degrees.

Various improvements and modifications can be made to the above without departing from the scope of the invention.

In particular, it is to be appreciated that the invention is not limited to electro-optical sources for terahertz radiation. It can be used with any type of radiation source, including any type of electro-optical or photoconductive source.

Furthermore, the invention has applicability in diverse fields such as specialist surveillance; 'filtered back projection'; the Radon transformation

as used in computerized tomography; scrutiny of mail; monitoring the progress (in location, time and chemical activity) of a reaction in a pharmaceutical process within a reactor vessel; monitoring the presence of impurities in a food mixture; or determining the nature and position of a foreign body within a fabricated plastic or composite structure. This list of examples is not exclusive or limiting.