Introduction to post-processing techniques

European Journal of Radiology 67 (2008) 202–217

Introduction to post-processing techniques Filip Jiru ∗ MR Unit, Department of Diagnostic and Interventional Radiology, Institute for Clinical and Experimental Medicine, Videnska 1958/9, 140 21 Prague, Czech Republic Received 29 February 2008; accepted 3 March 2008

Abstract The quality of data measured in in vivo MR spectroscopy is often insufficient due to a number of limitations such as low concentrations of observed metabolites and restricted measurement time resulting in a low signal-to-noise ratio. However, there are a variety of methods called post-processing techniques which allow the enhancement of the measured signal after measurement. In this review an introduction to the most important post-processing techniques for 1 H MR spectroscopy is given and practical examples are shown. In the first section the concept of FID and spectrum is introduced and the relationship between FID and spectrum is explained. Subsequently, the objectives and description of the following post-processing techniques are provided: eddy current correction, removal of an unwanted component (water), signal filtering for various purposes, zero filling, phase correction and baseline correction. © 2008 Elsevier Ireland Ltd. All rights reserved. Keywords: 1 H MR spectroscopy; Post-processing; Signal enhancement

1. Introduction Over the last few decades in vivo MR spectroscopy (MRS) has found an important place in experimental research as well as in clinical practice. To make use of MRS as a diagnostic tool a reliable quantitation of signals observed in MR spectra is necessary. Unfortunately, in vivo MRS examination has a number of limitations (e.g. low signal-to-noise ratio (S/N) due to usually low concentrations of observed metabolites and a restricted measurement time as well as poor homogeneity of magnetic fields) which make the measurement of high quality spectra challenging. There are, however, a variety of methods allowing the enhancement of the measured signal in different ways after measurement. Such methods, which will be further called postprocessing techniques, are the main topics of this paper. Most of the software used to quantify MR signal automatically include some post-processing techniques, therefore, the user often performs no additional corrections manually. However, for reliable analysis of in vivo MR spectra an understanding of the principles of post-processing techniques is necessary. This paper is organized as follows: firstly, parameters of free induction decay and spectrum and their correspondence ∗

Tel.: +420 23605 5235. E-mail address: [email protected].

0720-048X/$ – see front matter © 2008 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.ejrad.2008.03.005

are explained. Further, an overview of the most important post-processing techniques for 1 H MRS is given and practical examples are shown. The fundamentals of individual techniques are explained without an exact mathematical description to make the text more versatile for people with various mathematical backgrounds. Due to the limited scope of the paper only the main ideas of the techniques are presented and the readers are referred to the provided references for details. Reviews of post-processing techniques can be found in several papers [1–3]. 2. Free induction decay and spectrum The signal detected in the MR spectrometer in the receiving coil after exposing a sample to a radio frequency pulse is called Free Induction Decay (FID). In modern MR spectrometers the MR signal is detected using quadrature detection. As a result, the acquired MR signal is composed of two parts, often referred as real and imaginary parts of FID (the reason for this somewhat surprising notation is that complex numbers are used for the mathematical description of the signal), corresponding to the x and y components of the rotating magnetization vector M, respectively [4]. Since vector M is rotating around the z axis, the real and imaginary parts of FID are related to each other as sine and cosine (Fig. 1).

F. Jiru / European Journal of Radiology 67 (2008) 202–217

203

Fig. 1. Real (sx ) and imaginary (sy ) parts of FID (right) correspond to x and y components of the rotating magnetization vector M (left).

Let us now focus on FID composed of one magnetization component. There are four important parameters of such FID: amplitude A (which is proportional to the concentration of the studied component), frequency of oscillations ω, the function describing the decay of the FID’s amplitude (the envelope) characterized by the relaxation time T2* and the initial phase ␾0 (Fig. 2). Due to the way the signal is detected, the acquired FID oscillates with a frequency ω equal to the difference between the frequency ω (it should be noted that in this paper the word frequency stands for the angular frequency which is related to the frequency ν in Hertz as ω = 2πν) of the observed component and the adjustable reference frequency Ω of the phase detector (ω = ω − Ω). Therefore, FID with no oscillations can be observed whenever ω = Ω. The decay function is often assumed to be exponential in in vivo MRS. The relaxation time T2* then depends both on the natural relaxation time T2 of the observed component and on the homogeneity of the static magnetic field when poor homogeneity results in faster decay. The initial phase ␾0 of FID is related to the position of the M vector at the start of the FID acquisition and it is demonstrated as the shift of the FID (sine or cosine) along the time axis. It follows that, for example, for phase ␾0 = 90◦ the imaginary part equals the real part for the phase ␾0 = 0◦ .

Generally, the angle between the projection of the M vector to the transversal plane and the x axis is called the phase of FID at the given time. In Fig. 2 the real part of FID is shown; however, all parameters can be identified in the imaginary part analogously. For simplicity, only the real part of FID will be shown further in the text. The influence of FID parameters on its appearance is further demonstrated in Fig. 3 where FIDs with different values of the parameters are shown. If the sample contains more magnetization components (magnetization vectors differing in some of the parameters) the resulting FID equals the sum of particular FIDs corresponding to individual components, each with the given amplitude, frequency, T2* and initial phase. Many quantitation methods (so called time domain methods) enable analysis of the parameters of individual components directly from the resulting FID [2,5]. However, mainly for the visual inspection of a multi-component signal the visualization of FID is not useful. There is an alternative possibility to look at the signal: instead of observing time evolution of the signal in FID, a plot showing the individual components of the signal as a function of frequency can be drawn. Such plot is called a spectrum and the

Fig. 2. Parameters of one component FID: Amplitude A, frequency of oscillations ω, decay function with the decay constant T2* and the initial phase ␾0 . Rotation of the magnetization vector M in the transverse plane (left) and the real part of FID (right) are shown. Signal decay can be calculated using real (sx ) and imaginary (sy ) parts of FID as described in the figure.

204

F. Jiru / European Journal of Radiology 67 (2008) 202–217

Fig. 3. The influence of FID parameters on its appearance. FIDs with decreasing relative frequency ω (first row), FIDs with increasing T2* (second row), and different initial phases of FID (third row).

tool to transform FID into a spectrum is called a Fourier transform (FT) [6]. FT decomposes FID into individual oscillating components and shows each component as a peak with a given amplitude and position on the frequency axis. Similarly to FID, the spectrum also has real and imaginary parts, called absorption and dispersion Lorentz lines, respectively (Fig. 4). The unique property of the Fourier transform is that there is a direct correspondence between the parameters of FID and the parameters of the peak in the spectrum. Similar to FID, parameters in the frequency domain are encoded both in real (absorption Lorentz line) and imaginary (dispersion Lorentz line) parts of the

spectrum. Because absorption lines only have positive intensity and narrower tails, spectral analysis is usually performed using absorption lines. In detail, the area under the absorption peak is proportional to the amplitude A of the corresponding component in FID, the position of the peak in the frequency axis corresponds to the frequency of oscillations ω and the full width at half maximum (FWHM) is inversely proportional to T2* . It can be shown that the exact shape of the peaks depends on the decay function (the envelope) of FID. An exponential decay function leads to the aforementioned Lorentzian peak shapes as shown in Fig. 4. However, in the case of an inhomogeneous static magnetic

Fig. 4. The relation between one component FID and a spectrum. Real and imaginary parts of FID are transformed to real and imaginary parts of spectrum using Fourier transform. Note that both real and imaginary parts are needed to calculate real (or imaginary) parts of the spectrum.

F. Jiru / European Journal of Radiology 67 (2008) 202–217

205

Fig. 5. Spectrum of multi-component FID. Thanks to the linearity of the Fourier transform, the spectrum of multi-component FID is given by the sum of individual peaks corresponding to the particular FID components.

field the decay function is generally no longer exponential which results in non-Lorentzian peak shapes. Only in the special case of Lorentzian distribution of inhomogeneities, the peaks remain Lorentzian, but with increased widths (this corresponds to the introduction of T2* < T2, instead of T2). It should be noted that in the absorption Lorentz lines the height of the peak depends both on the amplitude A and T2* , unlike the peak’s area which depends only on the amplitude A. Therefore, all quantitation methods determining FID parameters from the spectrum focus on the determination of peak area. For multi-component FID the resulting spectrum is given by the sum of individual peaks corresponding to the particular FID components as shown in Fig. 5. FT operates in both directions; this means that the original spectrum can be transformed back to FID without any loss of information. In this sense, the description of the signal in both domains is equivalent. It should be noted that there are also other interactions not mentioned here, such as J-coupling, dipole–dipole interactions and others which additionally alter the signal intensity and shapes of the peaks.

from FID and spectrum equivalently. However, as shown later, the exhibition of possible artifacts in the particular domain may favor one of the domains for signal quantitation. It follows that different post-processing techniques may be required depending on the domain the signal is quantified in. Also, some postprocessing can be applied in both domains with equivalent or similar outcome. Therefore, the separation of post-processing into time domain methods and frequency domain methods is not so straightforward. Further we will focus on the following postprocessing techniques: eddy current (EC) correction, removal of an unwanted component (water), signal filtering for various purposes, zero filling, phase correction and baseline correction. Although the post-processing may enhance the signal, two things should be kept in mind. Firstly, the acquisition parameters should be set optimally to avoid extensive post-processing if possible. Secondly, the improved visual appearance of the signal achieved by post-processing may be at the expense of the accuracy parameter estimates because some processing methods may alter the parameters of FID as explained later. 3.1. Eddy current correction

3. Introduction to post-processing methods As already mentioned in Section 1 by post-processing we mean any signal manipulation performed in order to improve the visual appearance of the MR signal or the accuracy of estimates of signal parameters. A variety of techniques have been proposed to achieve this goal. Below, we will focus on those techniques which are important for 1 H MRS and which are implemented in clinical scanners. Thanks to the correspondence of FID and the spectrum, the parameters of the MR signal can be calculated in principle both

The localization of the signal in measurement sequences is accomplished using RF pulses applied together with magnetic field gradients, as described in ref. [7]. In addition, crusher gradients eliminating unwanted signal contributions are used. Therefore, many events of switching the gradients on and off occur in the sequence. According to Faraday’s law of electromagnetic induction, each change of the magnetic field induces a current in the conductive part exposed to the field change. Such currents induced in the conductive parts of the MR scanner are called ‘eddy currents’. The exact distribution and time depen-

206

F. Jiru / European Journal of Radiology 67 (2008) 202–217

Fig. 6. An exhibition of eddy currents during signal measurement. The magnetic field induced by eddy currents adds to the gradient fields and interacts with the measured signal.

dence of the magnitude of the eddy currents depend both on the switching pattern of the gradients and the geometry of the conductor. However, after the field perturbation is finished the magnitude of the associated current starts to decay to zero. According to Ampere’s law each electric current is associated with a magnetic field. As a result, eddy currents produce additional time dependent magnetic fields which add to the gradient fields used for slice selection and interact with the measured signal (Fig. 6). The first effect leads to the deterioration of the slice profile and to changes in the FID’s amplitude. In further explanations this effect will be ignored. The presence of eddy currents during data acquisition exhibits as time dependent phase shifts in the FID. After Fourier transform typical artifacts can be observed in the spectrum, as shown in Fig. 7. Distorted lineshapes caused by heavy eddy currents therefore prevent the accurate quantitation of the signal. Eddy current correction consists of the elimination of the additional time dependent phase shifts from the FID. As described in the Section 2, when no EC are present the signal component with a resonant frequency equal to the reference frequency Ω (we will call such a component ‘on resonance’ further) shows no oscillations and hence no time evolution of the phase in FID (compare to Fig. 3). It follows that the phase evolution of the on resonance component directly shows an addi-

Fig. 8. Scheme of eddy current correction.

tional eddy current related phase. EC free spectrum can be then retrieved by the subtraction of the phase calculated from the on resonance signal from the phase of the FID corrupted by eddy currents [8]. The common practice in 1 H MRS is the acquisition of an additional reference FID without water suppression and with the reference frequency Ω set to the water frequency, since in the water unsuppressed FID the signal of other metabolites can be neglected (water has a much higher concentration). In this way, the acquired FID can be taken for a signal from one component on resonance. Because of the high water concentration no signal averaging is usually needed and the additional water scan takes several seconds. The phase of the water FID is determined in each time point and it is subtracted from the phase of the corrupted FID. Assuming the amplitude of the FID is not affected by EC the resulting FID is corrected and the corresponding spectrum shows no artifacts. The whole procedure is depicted in Fig. 8. The described correction procedure relies on several assumptions. Firstly, the position and size of the selected volume must be

Fig. 7. Manifestation of eddy currents in FID and in the spectrum. FID affected by eddy currents (left) and corresponding artifact in the spectrum (right). Using a broken line, the corresponding FID without eddy currents is depicted in the FID.

F. Jiru / European Journal of Radiology 67 (2008) 202–217

identical in the reference and EC corrupted FIDs. Secondly, the same gradient schemes must be employed in both sequences. Therefore, only water suppression and the number of acquisitions should be changed for the acquisition of the reference FID. Thirdly, only one signal component (the water signal in this case) must be present in the reference FID. If this is not the case (the presence of lipid signals or multiple water peaks) the described correction fails and an oscillatory ringing artifact can be observed in the corrected spectrum unless additional reference FID processing is performed [9,10]. Alternative methods for EC correction without the need for an additional reference FID have also been proposed making use of the strong signal in the eddy current corrupted spectrum resonating at a different frequency than the remaining components or using the acquisition of two sequences with opposite gradients [11]. It should be noted that with the advent of new generation MR scanners equipped with active gradient shielding and gradient pre-emphasis the EC do not usually represent a significant issue in 1 H MR spectroscopy. 3.2. Removal of the water signal The removal of selected signals from FID is often of advantage. Firstly, elimination of selected signals may reduce the complexity of the signal analysis and improve the accuracy of the parameter estimates. Secondly, artifacts are usually scaled with the signal strength and therefore, for example, even small distortions of the tails of a high peak may disturb signals of the remaining resonances with lower intensities. A typical example is the water peak in 1 H MRS. The concentration of water in the human brain is several orders higher than the concentration of other metabolites. Therefore a common practice in in vivo 1 H MRS is to suppress the water signal during measurement. However, there are several reasons why additional post acquisition water removal is necessary or of advantage. Due to magnetic field inhomogeneities and the multi-compartment nature of the water signal, water suppression is never perfect and a residual water signal remains in the spectrum. The residual water signal often has a complicated lineshape which is problematic when the spectrum is to be processed by a processing method

207

requiring the parameterization of the signal peaks. Recently, techniques enabling measurement of signals of metabolites from FID acquired without water suppression have been reported [12,13]. This is of great interest since the water signal can be used as an internal standard for quantitation and for EC correction without the need for any additional reference scans. For these applications the effective post acquisition water removal methods are of great importance. Several strategies to remove a signal from FID can be used taking advantage of the difference in some properties of the peak to be removed and the remaining signals. Concerning the water signal, three distinct approaches, also applicable with some modification for the removal of other signals, will be described below: (1) the approximation of the water signal and its subtraction from the FID, (2) elimination of the water signal using special filters, (3) removal of the broad water peak in the spectrum by baseline correction. The key step of the first method is the correct approximation of the water signal. If the reference unsuppressed water signal is available it can be scaled and subtracted from the original FID. The scaling can be performed by comparing the water peaks in the reference and original spectra. However, the procedure works well only if the residual water peak in the original water suppressed signal is not distorted, which is usually not the case when strong water suppression has been performed. Alternatively, the water signal can be determined by the analysis of FID. For this purpose fast non-parametric methods are used, such as HankelLanczos singular value decomposition (HLSVD) [14], matrix pencil method [15] or other methods. HLSVD is the fast modification of the singular value decomposition (SVD) method, which decomposes FID into a sum of exponentials. It follows that several exponential signals may correspond to the water peak if it has a complicated lineshape. Therefore, to remove the water component, exponential signals in the frequency range around the water frequency are summed and subtracted from the FID (Fig. 9). This method has proved to be very efficient and it has been implemented in various processing programs. SVD (HLSVD) is a ‘black box’ method requiring only a few parameters as input. The optimal parameters for water suppression using SVD can be found in [16]. Alternatively, only one or

Fig. 9. Water removal using the HLSVD method. Spectrum before (left) and after (right) water removal. Note that both spectra are subject to different scalings of vertical axis.

208

F. Jiru / European Journal of Radiology 67 (2008) 202–217

Fig. 10. Filtering of FID by the multiplication of FID by a filter function h(t). Multiplication of each time point of FID by the corresponding value of the filter function h(t) yields the filtered FID.

first few components with the highest signal contributions can be removed from FID using HLSVD regardless their resonating frequency [17]. This method is very fast but it fails when the water peak has a complicated lineshape or when another signal with the signal strength comparable to the strength of the water signal is present in FID. The second possibility is to utilize the fact that water resonates at a different frequency than other metabolites. Therefore, filters filtering out a defined frequency range can be applied [18]. Assuming water resonates at the zero frequency (or it can be shifted to zero frequency subsequently), the water signal is obtained by low-pass filtering and subtracted from the FID or, alternatively, high-pass filtering of the FID yields a water free spectrum. It should be noted that the selection of the filter is of great importance because improper filtering may influence the remaining signals. For this reason special filtering procedures have been proposed combining filtering and some other signal treatments [19]. As a consequence, such filters have the ability to completely remove the water spectral line including its frequency domain tails while minimally affecting signals resonating in the frequency range of the tails. This cannot be achieved with conventional frequency domain filtering. If the water signal is not subtracted using one of the described methods, the peaks of the metabolites are usually superimposed

on the tail of the broad water peak. In this case the background water signal can be partially removed by baseline correction, as described below. 3.3. Signal filtering in time domain Some of the properties of FID can be changed by the multiplication of FID by special functions, called filters or windows (Fig. 10). The goal of the FID filtering (here filtering of FID is achieved by the multiplication by a filter, however, various mathematical operations can generally be used) is the enhancement or suppression of different parts of FID leading to improvement in the signal quality. In the following sections three different widely used filtering strategies will be described: sensitivity enhancement, resolution enhancement and ripple reduction (apodization). 3.3.1. Noise in the signal So far we have assumed that the recorded FID contains only the signal. However, as evident from any measurement, noise is always present in FID. By noise we mean a random function of time with fluctuations of certain standard deviation, which is added to the signal in FID. As a result, each point in the recorded FID is corrupted to some degree by the random noise

Fig. 11. Manifestation of noise in FID and in the spectrum. The noise in FID propagates into the spectrum through the Fourier transform. Although points at the early point of FID containing both noise and signal cannot be discarded, the acquisition of noise at the end of the FID should be avoided.

F. Jiru / European Journal of Radiology 67 (2008) 202–217

value. The origin of the noise in in vivo MRS is two-fold: (1) the thermal motion of the charged particles in the sample and (2) thermal motion of the electrons in the coil and reception paths. Under normal circumstances the latter is the major source of noise in the MR signal. The noise in the MR signal is assumed to be Gaussian noise with zero mean and constant standard deviation in time. To quantify the signal quality, the signal-to-noise ratio is widely used, and is usually defined as the amplitude of the signal divided by the standard deviation of the noise. The noise present in FID propagates after the FT into the spectrum. This is demonstrated in Fig. 11 where the same FID with increasing acquisition time tacq and the corresponding spectra are shown. From the figure it is clear that recording the FID after it has decayed to the noise level is just recording the noise and the corresponding spectrum shows a higher noise level the longer the noise is acquired. Therefore, acquisition time should always be adjusted so that no signal acquisition takes place when only noise remains in the FID. The presence of noise is an issue for signal quantification. The small resonances in the spectrum may get lost in the noise if the S/N is insufficient. Because S/N of the signal acquired in vivo is usually low, those methods which enable improvement of the S/N after measurement are of great importance. The post acquisition noise reduction by means of signal filtering is called sensitivity enhancement. 3.3.2. Sensitivity enhancement Although the signal has the greatest amplitude at the beginning of FID and decays to zero, the noise is constant throughout the entire FID (with constant standard deviation). It follows that the points at the end of FID carry more noise relative to the signal and less useful information than the points at the beginning. If the FID is multiplied by a filter equal to one at the

209

origin and slowly decaying to zero, the later parts of FID containing mainly noise will be suppressed but the early parts of the FID will be left almost unchanged. The value of the first point of the filter equal to 1 ensures that the amplitude A of the filtered FID remains unchanged, which is important for further quantitative analysis. Clear noise reduction in the spectrum is apparent in Fig. 12, where the described procedure is depicted. From Fig. 12 it is also apparent that noise reduction by signal filtering is at the expense of peak width. This is understandable because each multiplication of FID by a decaying function introduces faster decay of FID and hence broader peaks. Assuming the area under the peak remains constant, the broader peak leads to a lower peak. Therefore, the noise is reduced on the one hand, but the peak height is lowered on the other. As a result, an optimal filter exists which leads to the highest possible S/N in the spectrum. Such a filter is called a matched filter and it can be shown that the function describing the matched filter equals the function describing the signal decay (envelope) of the given FID [20]. For example, for exponential decay the matched filter is an exponential function with the decay function equal to T2* of the FID. This corresponds to the linewidth of the filtered peak equal to double the original value. The selection of the matched filter is not so straightforward for multi-compartment FID where compounds have different T2* values. In this case one value is selected and the condition of the matched filter is fulfilled only for one compound. Actually, exponential functions with slower decay than corresponding to the matched filter condition are often used in vivo because two-fold linewidth broadening is often not acceptable. Also, different filter functions, such as Gaussian function are often used for filtering. However, it should be noted that assuming exponential decay of the original FID, the only filter function ensuring Lorentzian peak shape after filtering is the exponential function. The selection of the other filter functions generally results in the modification of the line-

Fig. 12. Sensitivity enhancement of FID. Multiplication of noisy FID by exponential function leads to the reduction of noise in the resulting spectrum (right). However, the noise reduction is at the expense of line broadening as apparent from comparison of the original (left) and filtered (right) spectrum.

210

F. Jiru / European Journal of Radiology 67 (2008) 202–217

shapes, which must be taken into account during quantitative analysis. 3.3.3. Resolution enhancement Due to inoptimal homogeneity of the magnetic field in the examined volume, the linewidth of the peaks in vivo is much broader compared to the high resolution NMR. The broad linewidths are problematic especially for peaks resonating at close frequencies because the peaks overlap and the spectral resolution degrades. From the correspondence between the linewidth and decay function of the FID it follows that the linewidth (and the lineshape in general) can be modified by modification of the decay function. The artificial narrowing of linewidths in the spectrum by means of post acquisition filtering of the signal is called Resolution enhancement. There are several approaches to achieve narrower linewidths in the spectrum. The simple empirical filter used for resolution enhancement is called the convolution-difference filter. This method directly utilizes the fact that the low resolution of the spectrum is a consequence of the overlaps of broad resonances and it removes the components of the FID giving rise to the broad peaks in the spectrum (see Section 3.6 for details). The filter has the form   −t h(t) = 1 − a exp (1) b where a and b are empirical constants adjusting the sensitivity of the filter to the broad peak components. The next class of resolution enhancement methods aims to ‘undo’ the signal decay. To understand the idea imagine that one component FID is multiplied by the increasing exponential function with a parameter equal to T2* of the original FID. The filtered FID is then a rectangular FID with no decay and a very narrow linewidth in the spectrum. However, in practice two restrictions prevent doing this. Firstly, the discontinuity at the end of FID causes truncation artifacts leading to finite linewidth and oscillations in the spectrum as described further. Secondly, noise contributions in the later parts of the signal are excessively enhanced leading to unacceptably noisy spectra. To overcome the issues alternative filter functions increasing in the initial parts of FID and decaying to zero at the end are utilized. Generally, knowing the original lineshape, it is possible to achieve arbitrary lineshape in the spectrum through proper filter design. One of the commonly used techniques for changing lineshape is Lorentz–Gauss transformation. Providing the original FID decay function is exponential, the Lorentz–Gauss transformation transforms the original Lorentz peak to a Gaussian peak. The rationale behind this is the fact that Gaussian lines are much narrower at the base than Lorentzian lines with the same width at the half maximum (at 1% of the peak amplitude the Gaussian line is almost five times narrower). The Lorentz–Gauss filter has the form   t t2 h(t) = exp − (2) a b

The parameter a should equal to the T2* of the FID to transform the original Lorentzian peak to a pure Gaussian peak. The parameter b then controls the degree of the peak narrowing. For multi-compartment FID, trade-off values of a, b must be selected. A pure Gaussian peak shape will then only be achieved for the signal component corresponding to the selected parameter a and lineshapes of the other peaks will be a mixture of Lorentz and Gaussian functions. Alternatively, repeated spectra processing with a, b adjusted for the given peak of interest can be performed. Comparison of the Lorentz–Gauss filter and increasing exponential function for resolution enhancement is shown in Fig. 13. As a result of an inhomogeneous magnetic field or the occurrence of possible artifacts, the measured peaks may have lineshapes different from the Lorentzian lineshapes. Various methods have been proposed to determine the experimental lineshape from a reference peak (usually the water peak) and to transform the lineshapes of the peaks in the spectrum to the desired shapes [21,22]. It should be noted that selection of filters for resolution enhancement and their constants may influence quantitative analysis of the signal. Therefore, they should be used with great care. Conversely, for visual inspection of badly shimmed spectra the resolution enhancement is beneficial. 3.3.4. Apodization As described in Section 3.3.1, for sensitivity reasons the acquisition time should be adjusted so that no signal acquisition takes place after the signal has decayed to zero. Sometimes the acquisition time is too short and the acquisition is finished before the signal decays to the noise level. We talk then about signal truncation which is common in 2D NMR and may occur in phantom measurements where T2* values are usually much longer then in vivo. Truncation is depicted in Fig. 14, where the spectrum calculated by the FT of the truncated FID is also shown. Typical oscillatory signal tails (‘ripple’) can be seen in the spectrum as a consequence of the signal truncation. The more severe the truncation, the higher the amplitude of the ripple is. In reality, no ripple is observed when Fourier transform (discrete Fourier transform) of the sampled FID signal is performed. This is because the calculated points in the spectrum correspond approximately (depending on the ratio tacq /T2* ) to the zero crossing of the oscillations as shown in Fig. 14. However, the ripple becomes visible when zero filling is performed (see next section). From the properties of the Fourier transform it follows that the oscillations stem from the discontinuity at the end of FID [6]. Therefore, the ripple can be suppressed by removing the discontinuity by FID filtering, which in this context is called apodization. To apodize the FID various filters equal to 1 at the first point and decaying to zero at the end of the FID can be used as in sensitivity enhancement. Also, similar to sensitivity enhancement, the ripple reduction by apodization is at the expense of increased linewidths. The optimal filter in terms of ripple reduction and line broadening cannot be described analytically [20] but it can be approximated well using the Hamming

F. Jiru / European Journal of Radiology 67 (2008) 202–217

211

Fig. 13. Comparison of signal filtering using increasing exponential function and Lorentz–Gauss filter for noise free (left hand trace) and noisy (right hand trace) spectra. Original high resolution spectra (first row), artificially broadened spectra mimicking bad shimming (second row), filtered spectra using increasing exponential function (third row), filtered spectra using Lorentz–Gauss filter (fourth row). Filters were applied to the broadened spectra shown in the second row. Different parameters of the filter were used for noisy and noise free spectra in the third row (increasing exponential function) to achieve optimal resolution enhancement. Note that the high noise level in the filtered spectrum limited achievable resolution despite the very high signal-to-noise ratio of the original broadened spectrum.

212

F. Jiru / European Journal of Radiology 67 (2008) 202–217

Fig. 14. The truncated FID (left) gives rise to the oscillatory signal tails in the spectrum (right). The solid line in the spectrum represents a continuous signal, which can be approached by unrestricted zero filling. The dots represent sampling points given by Fourier transform of the original FID without zero filling. It follows that the ripple is almost invisible in this case. However, the more zeroes that are added to FID by zero filling, the more the structure of oscillations become visible.

3.4. Zero filling

filter  h(t) = 0.54 + 0.46 cos

πt tacq

 (3)

where tacq is acquisition time of FID. As an alternative to the filtering, the maximum entropy [23] or linear prediction [24] methods can be used to extrapolate the signal after truncation. It follows that different filters are used to achieve different goals. Given that in practice only one filter is applied to the data, a compromise in the filter selection in terms of sensitivity enhancement, ripple reduction and resolution enhancement has to be made. It can be shown [20] that the signal filtering achieved by the multiplication of FID with the time dependent filter function corresponds to the linear filtering methods ensuring no interference effects between overlapping resonances during the filtering process takes place. In conclusion, it should be mentioned that signal filtering introduces the correlation of noise. As a result, some of the statistical methods used for signal fitting which assume that the noise in the signal is uncorrelated may give false results if strong filtering has been applied to the data prior to the fitting.

During the acquisition of FID the signal is sampled point by point using an analogue–digital converter. Time spacing between successive points is determined by the desired bandwidth of the spectrum and the number of points by the decay rate of the FID, as pointed out in the Section 3.3.1. Because the Fourier transform preserves the number of points, the number of points (values) in the spectrum and in the FID are equal. It should be pointed out that plots of FID and of the spectrum are created by connecting calculated points. Therefore, if the number of points is not sufficient the representation of the signal fails. This is relevant especially for sharp peaks in the spectrum where a high density of points is necessary to correctly represent peak heights, shapes and positions, as demonstrated in Fig. 15. Moreover, the same misrepresentation happens during the signal fitting, which may lead to systematic quantitative errors. Although the number of points is dictated by the measurement parameters, a possibility exists to increase the number of points in FID artificially after measurement. The technique allowing this is called zero filling and it consists of the appending of zeros after the last measured point in FID. As a result the spectrum calculated by the Fourier transform of the FID with appended zeros is interpolated to the total number of points equal to the

Fig. 15. The importance of the number of points for accurate characterization of peaks in the spectrum. Low density of points may result in misrepresentation of peak height, shapes or position (left). When the sampling density is not sufficient the peak shape depends strongly on its position with respect to the sampled points (right)—observed peak shape (solid line) and the real peak shape (broken line) differ.

F. Jiru / European Journal of Radiology 67 (2008) 202–217

213

not able to make up for the missing signal, as can be seen in Fig. 17 where truncated FID was zero filled. The resulting ripple can be suppressed by apodization but only at the expense of line broadening and hence loss of resolution. Conversely, zero filling is always a method of choice when a longer acquisition time would result in the mere acquisition of noise. To utilize FFT, the fast algorithm for the calculation of the Fourier transform of discretely sampled values, the number of points in the FID must be a power of 2. Zero filling is therefore often used to increase the number of points to a value which is a power of 2 if an arbitrary number of time points is acquired. In conclusion, zero filling is a useful method allowing the improvement of digital resolution of the spectra which is of importance both when displaying spectra and for signal analysis. To avoid ripple of the baseline, apodization needs to be performed prior to zero filling if the FID has not decayed to the noise level at the end of the acquisition window. 3.5. Phase correction Fig. 16. Zero filling. Original FID and spectrum (top) and FID and spectrum after 1× zero filling (bottom). Appending of zeros after measured FID results in more spectral points and hence in better characterization of the peak shapes.

original number of points plus the number of appended zeros. An illustration of the technique is shown in Fig. 16. If the signal is properly sampled (a spectral range of the spectrum covers frequencies contained in FID) the sampling theorem states that the signal (FID in our case) is completely determined by its samples and can be fully reconstructed [6]. The procedure of zero filling represents a suitable interpolation process enabling the calculation of missing values from the acquired values. Appending zeros to the acquired FID brings no new information to the signal. In this respect zero filling only makes use of the information already present in the acquired signal. Therefore, zero filling does not improve the true resolution of the spectrum given by the acquisition parameters and the lineshape of the peaks. However, it does improve digital resolution and definitely helps to improve the characterization of the peak shapes in the spectrum, as apparent in Fig. 16. Resolution is generally given by the length of the acquisition time. Although zero filling allows artificially increasing acquisition time it is

In Section 2 it was shown that the Fourier transform of FID with zero initial phase yields a complex spectrum with absorption Lorentz peaks as a real part and dispersion Lorentz peaks as an imaginary part. However, the situation dramatically changes if the initial phase in the FID in non-zero. The appearance of the spectrum for different initial phases ␾0 is shown in Fig. 18. It is apparent that neither the real nor the imaginary part of the spectrum has absorption (or dispersion) lineshapes. Instead, both real and imaginary parts are a mixture of absorption and dispersion Lorentz lines. In measurements, the initial phases of different compounds in FID are generally non-zero. Therefore, in order to restore pure absorption and dispersion shapes a special correction called a phase correction must be applied. Mathematically the correction consists of the multiplication of the complex spectrum by the complex phase factor equal to the initial phase in FID. Since the value of the phase factor is a priori not known, the correction is performed interactively by observing the real (imaginary) part of the spectrum and multiplying the spectrum by the variable phase factor θ until the absorption (dispersion) shape is restored. It can be shown that the phase factor θ has two components, one which is not dependent on frequency and the second which

Fig. 17. Example of zero filling of truncated FID. Spectrum calculated from FID sampled to 1024 points (left), spectrum calculated from the same FID but sampled only to 32 points (middle), spectrum calculated from the FID sampled to 32 points and subsequently zero filled to 1024 points (right). Note that zero filling resulted in oscillations in the spectrum as a consequence of the signal truncation in the FID sampled to a low number of points.

214

F. Jiru / European Journal of Radiology 67 (2008) 202–217

Fig. 18. Appearance of real and imaginary parts of FID as well as real and imaginary parts of corresponding spectrum as a function of initial phase ␾0 of FID.

is frequency dependent and, therefore, different for each peak in the spectrum. The frequency independent part of the factor θ 0 is called the zero order phase or absolute phase. The frequency dependent factor θ 1 , which in a good approximation depends linearly on the frequency, is called the first order phase or linear phase. The overall phase factor θ required for the given frequency can be, therefore, written as θ = θ0 + θ1 (ω) = θ0 + τ(ω − ωcent )

(4)

where the constant τ determines how steeply the linear part of the factor changes with the frequency and the constant ωcent determines for which frequency ω the linear phase will be zero. It follows that to correctly phase the spectrum, the zero order and first order phases have to be adjusted. The adjustment of the overall phase is performed separately for zero and first order phases. In the first case we then talk about zero order (absolute) phase correction, in the second case about first order (linear) phase correction.

Fig. 19. Absolute and linear phase corrections of the spectrum. Absolute phase correction influences all peaks in the spectrum in the same way. Spectra corresponding to absolute phases 90◦ and 180◦ are shown (top). Linear phase correction influences the peaks depending on the central frequency ωcent and the parameter τ. Note that the peak with frequency ω = ωcent remains unchanged (bottom).

F. Jiru / European Journal of Radiology 67 (2008) 202–217

The impact of the phase parameters on the spectrum appearance is shown in Fig. 19. Although phasing of the spectrum can be achieved in many ways, one of the possibilities is the phasing of one selected peak using absolute phase correction only, then setting the parameter ωcent equal to the frequency of the phased peak and eventually phasing the whole spectrum by changing parameter τ. An alternative is the calculation of the linear phase from phases necessary to phase peaks at different frequencies. Automatic phase correction without operator intervention can be performed using several algorithms with varying degrees of success [25,26]. Two common problems (excluding the distortions caused by eddy currents) may be encountered when phasing the spectrum. Firstly, some of the peaks may show a different phase than the others, the second spectrum cannot be phased unless baseline distortion is introduced. The first problem is caused either by the existence of a peak with phase modulation due to spin–spin coupling (e.g. lactate) or the presence of a peak which has been affected during the magnetization preparation period. An example is a water signal composed of several differently suppressed signals corresponding to individual water compartments. The second problem is given by the fact that the first order phase correction does not work precisely when the linear phase is large. This is because the linear phase of the given peak depends on the centre frequency of the peak while the spectral frequency coordinate is used for phasing the spectra. As a result frequency dependent phasing does not work precisely in the wings of the peaks and in the regions where peaks overlap. In this case rolling baselines can be seen in spectra after phase correction. Methods to deal with the rolling baseline observed after phasing delayed acquisition spectra have been proposed [27,28]. The existence of constant and frequency dependent initial phases in FID has several reasons. The non-zero constant phase

215

can be caused by the maladjustment of the phase detector or it can be a result of phase cycling. First order phase is a consequence of the delay between the excitation and the detection of FID, the frequency response of the receiver filter and the off-resonance excitation. 3.6. Baseline correction Peaks in the spectrum often do not lie on a flat line with zero intensity. The unknown background in the spectra is called the baseline (Fig. 20). The existence of the baseline makes quantitative analysis difficult because the estimation of peak areas in the presence of an unknown background is problematic. The background signals are usually broad signals with very short T2* values. Therefore, most of the baseline signals can be identified as fast decaying components in the early parts of FID as illustrated in Fig. 21. In measurements two distinct reasons for the rough baseline can be generally found: hardware imperfections and signals originating from the sample. If there is an insufficient delay between the RF pulse execution and signal reception, the interaction of the transmitter with the receiver can occur resulting in the corruption of several initial points in FID. Although this artifact is due to the fast electronic of modern MR scanners not common, it is an illustrative example how the modification of a limited number of points in FID changes the appearance of the whole spectrum. In Fig. 22 spectra corresponding to the same FID but with changed first, first two and first three points are shown. As apparent, the alteration of one point in the FID causes a shift of all points in the spectrum above the zero line, two altered points introduce humps in the spectrum and three altered points result in baseline rolling. This behavior reflects the fact that the spectrum is calculated by Fourier transform where each single

Fig. 20. Baseline in the spectrum. Spectrum with an ideal (left) and rough (right) baseline.

Fig. 21. Identification of baseline signals in FID. The broad signal s2 in the spectrum (right) can be identified as a quickly decaying component of FID (left).

216

F. Jiru / European Journal of Radiology 67 (2008) 202–217

Fig. 22. Modification of few early points in FID leads to artifacts in the spectrum. One point in FID causes a shift of all points in the spectrum above the zero line (left), two altered points introduce humps in the spectrum (middle) and three altered points result in baseline rolling (right).

point in the spectrum is calculated using the information from the whole FID. A far more common source of rough baseline is the existence of broad signals in the measured sample. For example, macromolecules and lipids both have short T2* values giving rise to a typical baseline pattern in 1 H MRS of the brain. It should be pointed out that because components with short T2* contribute only very little at longer echo times, measurements with longer echo times lead to the simplification of the baseline. Also, inefficient water suppression causes the spectrum to be superimposed on the tail of the water peak unless post acquisition water removal is performed. The post acquisition removal of the background signal in the spectrum is called baseline correction. The correction consists of the subtraction of the function describing the course of the background signal, as depicted in Fig. 23. However, even though

the trend of the function may be obvious the exact determination of the function is difficult due to the existence of noise. In principle, several ways to remove the baseline are possible. A common approach is to mark several points in the spectrum belonging to the baseline and let the program fit a smooth line through the points. Modification of the method consists of marking areas where any peaks are present instead of baseline points. The processing program then calculates the baseline as an average value from the remaining noisy spectral points. Line fitting with splines or polynomials is usually employed. Although an almost arbitrary baseline can be fitted, the ambiguity of both the selection of the baseline points and of the peak marking makes this method highly subjective. Automated methods identifying the baseline in the spectrum have been proposed; however, selection of the threshold parameter is necessary.

Fig. 23. Baseline correction. The straight baseline in the spectrum (right) is accomplished by the subtraction of a function describing the baseline (middle) from the original spectrum (left).

Fig. 24. Incorporation of baseline signals into the model function used for signal fitting. Baseline (broken line) as a result of baseline correction when no background signals were included into the model function used for signal fitting (left) and with signals of lipids and macromolecules included in the model function (right). Note that the incorporation of baseline signals into the model function leads to the flatter baseline and that both approaches may lead to the different estimation of signals of metabolites of interest.

F. Jiru / European Journal of Radiology 67 (2008) 202–217

An alternative technique is the removal of the portion of FID containing broad signal components. A convolution-difference filter can be used for this purpose, as described in Section 3.3.3. Simple elimination of several initial points of FID removes baseline components but only at the expense of the occurrence of the rolling artifact (see explanation above) unless a special signal treatment, such as linear prediction, is employed. The main disadvantage of both approaches is that not only baseline signals but all signals in the FID will be influenced by this procedure. The most appealing and accurate way of coping with the baseline is the incorporation of baseline signals into the model function used for signal fitting (Fig. 24). This can be done in two ways. If individual background signals can be identified and parameterized they are added into the model function as additional metabolites [29]. Alternatively, the whole baseline pattern can be measured using some dedicated techniques suppressing signals of other metabolites and scaled in the process of the fitting. Baseline correction is often a source of systematic errors. Therefore, the proper baseline correction method is of importance when accurate quantification of signals is desired. 4. Conclusion Various post-processing methods allowing the enhancement of the measured signal after measurement have been described. The proper use of the techniques helps to improve the quality of acquired data, which is of great importance in clinical examinations. For an in-depth understanding of post-processing techniques both theoretical and practical approaches are necessary. Therefore, readers are encouraged to employ the described techniques in practice. Acknowledgements Supported by research projects MZ0IKEM2005 and LC554 of Ministry of Education, Youth and Sports, Czech Republic. References [1] Lindon JC, Ferrige AG. Digitisation and data processing in Fourier transform NMR Progress in NMR Spectroscopy 1980;14:27–66. [2] Vanhamme L, Sundin T, Hecke PV, Huffel SV. MR spectroscopy quantitation: a review of time-domain methods. NMR Biomed 2001;14(4):233–46. [3] Keeler J. Understanding NMR Spectroscopy http://www.spectroscopynow. com/coi/cda/landing.cda?chId=5&type=Education. [4] Haacke EM, Brown RW, Thompson MR, Venkatesan R. Magnetic resonance imaging, physical principles and sequence design. New York: John Wiley & Sons; 1999. [5] Helms G. The principles of quantification applied to in vivo proton MR spectroscopy. Eur J Radiol 2008;67(2):218–29. [6] Brigham E. Fast Fourier transform and its applications. Prentice Hall; 1988. [7] Klose U. Measurement sequences for single voxel protonMRspectroscopy. Eur J Radiol 2008;67(2):194–201. [8] Klose U. In vivo proton spectroscopy in presence of eddy currents. Magn Reson Med 1990;14(1):26–30.

217

[9] Wild JM. Artifacts introduced by zero order phase correction in proton NMR spectroscopy and a method of elimination by phase filtering. J Magn Reson 1999;137(2):430–6. [10] Simonetti AW, Melssen WJ, van der Graaf M, Heerschap A, Buydens LM. Automated correction of unwanted phase jumps in reference signals which corrupt MRSI spectra after eddy current correction. J Magn Reson 2002;159(2):151–7. [11] Lin C, Wendt 3rd RE, Evans HJ, Rowe RM, Hedrick TD, LeBlanc AD. Eddy current correction in volume-localized MR spectroscopy. J Magn Reson Imaging 1994;4(6):823–7. [12] Dong Z, Dreher W, Leibfritz D. Toward quantitative short-echo-time in vivo proton MR spectroscopy without water suppression. Magn Reson Med 2006;55(6):1441–6. [13] Dreher W, Leibfritz D. New method for the simultaneous detection of metabolites and water in localized in vivo 1H nuclear magnetic resonance spectroscopy. Magn Reson Med 2005;54(1):190–5. [14] Pijnappel WWF, van den Boogaart A, de Beer R, van Ormondt D. SVD-based quantification of magnetic resonance signals. J Magn Reson 1992;97(1):122–34. [15] Lin YY, Hodgkinson P, Ernst M, Pines A. A novel detection-estimation scheme for noisy NMR signals: applications to delayed acquisition data. J Magn Reson 1997;128:30–41. [16] Cabanes E, Confort-Gouny S, Le Fur Y, Simond G, Cozzone PJ. Optimization of residual water signal removal by HLSVD on simulated short echo time proton MR spectra of the human brain. J Magn Reson 2001;150(2):116–25. [17] Zhu G, Smith D, Hua Y. Post-acquisition solvent suppression by singularvalue decomposition. J Magn Reson 1997;124(1):286–9. [18] Coron A, Vanhamme L, Antoine JP, Van Hecke P, Van Huffel S. The filtering approach to solvent peak suppression in MRS: a critical review. J Magn Reson 2001;152(1):26–40. [19] Sundin T, Vanhamme L, Van Hecke P, Dologlou I, Van Huffel S. Accurate quantification of (1)H spectra: from finite impulse response filter design for solvent suppression to parameter estimation. J Magn Reson 1999;139(2):189–204. [20] Ernst RR, Bodenhausen G, Wokaun A. Principles of nuclear magnetic resonance in one and two dimensions. Oxford University Press; 1990. [21] Bartha R, Drost DJ, Menon RS, Williamson PC. Spectroscopic lineshape correction by QUECC: combined QUALITY deconvolution and eddy current correction. Magn Reson Med 2000;44(4):641–5. [22] Webb P, Spielman D, Macovski A. Inhomogeneity correction for in vivo spectroscopy by high-resolution water referencing. Magn Reson Med 1992;23(1):1–11. [23] Viti V, Massaro E, Guidoni L, Barons P. The use of the maximum entropy method in NMR spectroscopy. J Magn Reson 1986;70(3):379–93. [24] Barkhuijsen H, de Beer R, Bov´ee WMMJ, van Ormondt D. Retrieval of frequencies, amplitudes, damping factors, and phases from timedomain signals using a linear least-squares procedure. J Magn Reson 1985;61(3):465–81. [25] Chen L, Weng Z, Goh L, Garland M. An efficient algorithm for automatic phase correction of NMR spectra based on entropy minimization. J Magn Reson 2002;158:164–8. [26] Heuer A. A new algorithm for automatic phase correction by symmetrizing lines. J Magn Reson 1991;91:241–53. [27] Allman T, Holland GA, Lenkinski RE, Charles HC. A simple method for processing NMR spectra in which acquisition is delayed: applications to in vivo localized 31P NMR spectra acquired using the DRESS technique. Magn Reson Med 1988;7(1):88–94. [28] Saeed N, Menon DK. A knowledge-based approach to minimize baseline roll in chemical shift imaging. Magn Reson Med 1993;29(5):591–8. [29] Seeger U, Klose U, Mader I, Grodd W, N¨agele T. Parameterized evaluation of macromolecules and lipids in proton MR spectroscopy of brain diseases. Magn Reson Med 2003;49(1):19–28.