HARDWARE AND INSTRUMENTATION Full Paper Magnetic Resonance in Medicine 69:583–593 (2013) Gradient System Characterization by Impulse Response Measurements with a Dynamic Field Camera Signe J. Vannesjo,1 Maximilan Haeberlin,1 Lars Kasper,1,2 Matteo Pavan,1 Bertram J. Wilm,1 Christoph Barmet,1 and Klaas P. Pruessmann1* This work demonstrates a fast, sensitive method of characterizing the dynamic performance of MR gradient systems. The accuracy of gradient time-courses is often compromised by field imperfections of various causes, including eddy currents and mechanical oscillations. Characterizing these perturbations is instrumental for corrections by pre-emphasis or post hoc signal processing. Herein, a gradient chain is treated as a linear timeinvariant system, whose impulse response function is determined by measuring field responses to known gradient inputs. Triangular inputs are used to probe the system and response measurements are performed with a dynamic field camera consisting of NMR probes. In experiments on a whole-body MR system, it is shown that the proposed method yields impulse response functions of high temporal and spectral resolution. Besides basic properties such as bandwidth and delay, it also captures subtle features such as mechanically induced field oscillations. For validation, measured response functions were used to predict gradient field evolutions, which was achieved with an error below 0.2%. The field camera used records responses of various spatial orders simultaneously, rendering the method suitable also for studying cross-responses and dynamic shim systems. It thus holds promise for a range of applications, including pre-emphasis optimization, quality assurance, and image reconstruction. Magn Reson Med C 2012 Wiley Periodicals, Inc. 69:583–593, 2013. V Key words: gradient impulse response; GIRF; linear timeinvariant systems; mechanical oscillations; gradient calibration; dynamic shimming; magnetic field monitoring Time-varying magnetic field gradients are essential for signal preparation and encoding in magnetic resonance imaging and spectroscopy. Most MR methods rely on highly accurate gradient time-courses for correct signal encoding and suffer from artifacts when significant deviations from the prescribed time-courses occur. In practice, effective gradient waveforms usually do differ somewhat from the ideal shapes defined in the underlying pulse sequence. These deviations are largely due to a variety of 1 Institute for Biomedical Engineering, University and ETH Zurich, Zurich, Switzerland. 2 Laboratory for Social and Neural Systems Research, University of Zurich, Zurich, Switzerland. Grant sponsor: Philips Healthcare. *Correspondence to: Klaas P. Pruessmann, Ph.D., Institute for Biomedical Engineering, University and ETH Zurich, Gloriastrasse 35, CH-8092 Zurich, Switzerland. E-mail: pruessmann@biomed.ee.ethz.ch Received 23 November 2011; revised 3 February 2012; accepted 25 February 2012. DOI 10.1002/mrm.24263 Published online 12 April 2012 in Wiley Online Library (wileyonlinelibrary. com). C 2012 Wiley Periodicals, Inc. V hardware imperfections including bandwidth limitations of gradient amplifiers, eddy currents induced in gradient coils and in other conducting structures of the scanner (1,2), field fluctuations caused by mechanical vibrations after gradient switching (3,4), and thermal variation in hardware components (5). Slight perturbations can also stem from physiologically induced fields that originate in the subject under examination (6,7) or from magnetic sources and currents external to the MR system. Besides further hardware optimization, the most common ways of addressing dynamic field imperfections are precompensation of gradient waveforms (8–10) and postcorrection of acquired data (11). Both of these options are most feasible for mechanisms of perturbation that are reproducible and can be accurately modeled. Physiological and external field contributions as well as thermal effects are challenging to address in this way. However, eddy currents, which are usually the dominant causes of field imperfections, as well as mechanical effects have been successfully modeled and predicted. Eddy currents arise due to the self- and cross-inductance of gradient coils as well as the cross-inductance between gradient coils and other conductive structures of the scanner, such as shim coils or parts of the cryostat. To model eddy currents, the conductors involved are usually approximated as LR circuits, i.e., as combinations of an inductance and a resistance, in which induced currents decay exponentially (1,2,12). Using linear systems theory, the net step response of a gradient chain can then be modeled as a superposition of an ideal response, i.e., a Heaviside function, and a set of exponentially decaying field terms representing eddy currents of different sources. Response models of this kind are traditionally used to predict eddy current effects and to determine pre-emphasis settings for a given system (8,9,12,13). Field oscillations due to mechanical vibrations are a well-known disturbance in MR spectroscopy, where they can cause modulation sidebands particularly of unsuppressed water peaks. Similar to eddy current effects, the underlying mechanics and conversion from vibrations to magnetic fields have been viewed as linear mechanisms that contribute to a net linear impulse response (14). Using this approach, mechanically mediated field oscillations have been modeled as exponentially damped sinusoidal responses, which have been used to simulate sideband effects (4,15). Linear systems theory has also been used to characterize the acoustic response of a gradient system, relying on sound pressure measurements in the magnet bore (16–18). 583 584 So far, the models of field responses due to eddy currents and mechanical effects have relied on simplifying abstractions. For both mechanisms, it has been assumed that only a limited number of distinct instances occur, which obey an elementary pattern such that they can be represented by small sets of exponentials and exponential sinusoids, respectively. In general, this picture will hold only approximately. Due to the electrical and mechanical complexity of gradient systems and MR systems as a whole, both eddy-current and mechanical behavior may require expansion into very many components of different time constants. Eddy currents may be supported by structures that cannot be fully modeled as mere LR circuits. Moreover, coupling may occur both among structures that can carry eddy currents and among mechanical resonances, giving rise to yet more complex behavior. To better capture the variety of possible behavior, it is attractive to remove model assumptions on individual eddy-current processes or mechanical resonances and rather characterize a gradient system in terms of its complete impulse response function (19,20). According to linear systems theory, this approach should permit jointly representing all response mechanisms that are linear and time-invariant (LTI). A net gradient impulse response function (GIRF) should hence incorporate all influences on the gradient waveform between the console and the magnet bore. This would include amplifier and coil characteristics as well as eddy currents and vibration-induced fields, without the need to consider individual underlying mechanisms. Knowledge of the comprehensive GIRF could form the basis of advanced pre-emphasis and serve for quality assurance purposes. It could also yield more accurate estimates of effective kspace trajectories for image reconstruction and of other encoding parameters such as b-values in diffusion imaging or gradient moments in velocity mapping. The key challenge toward this goal is determining GIRFs accurately, with sufficient bandwidth and frequency resolution, and within reasonable measurement times. Probing the GIRF must generally involve observing a system’s response to given gradient input waveforms. This could be achieved with known NMR-based techniques for measuring gradients and k-space trajectories (21– 26). Partly, these methods rely on selective excitation of larger phantoms or imaging objects and signal detection with large receiver coils (22–24). Such detection methods offer only limited sensitivity, while the extent of the imaging object makes it prone to dephasing by potential inplane gradient responses. Methods that rely on a dedicated small sample are also limited in sensitivity unless the sample is mounted in a tight-fitting detector with a high filling factor (21,25,26). All of these methods have the disadvantage that each measurement needs to be repeated several times either to achieve temporal resolution or to distinguish different spatial field components. An efficient alternative is to record the field evolution with a dynamic field camera (27), i.e., with an array of miniature NMR probes that are operated simultaneously and positioned such as to distinguish different spatial components of interest (28,29). This approach was recently proposed for monitoring magnetic field evolutions concurrently with image acquisition (30). However, Vannesjo et al. it is equally suited for mere field observations as required for GIRF measurements. Relying on tightly wound detectors around small sample droplets, it combines high sensitivity with robustness against dephasing. The temporal resolution of the camera measurement is limited only by the acquisition bandwidth of the spectrometer used and the simultaneous recording at different positions immediately yields differentiation of a gradient’s self- and cross-responses. In this work, it is proposed to determine comprehensive gradient impulse response functions using field observations with a dynamic field camera. Starting from LTI systems theory, a strategy is derived for obtaining a full GIRF from suitable combinations of input functions. The method is demonstrated by GIRF measurements on a 3 T whole-body human MRI system and validated by comparing measured field evolutions with GIRF-based predictions. THEORY AND METHODS Probing an LTI System For a linear time-invariant system, the relation between the input to the system, i(t), and its output, o(t), is determined by the impulse response, h(t), as described by (31): oðtÞ ¼ Z1 FT iðtÞ  hðt
tÞdt $ OðvÞ ¼ IðvÞ  HðvÞ; ½1
1 where t and v denote time and angular frequency, and O(v), I(v), and H(v) represent the Fourier transforms (FT) of o(t), i(t), and h(t), respectively. Accordingly, the system response to any given input can be predicted based on the impulse response function. Knowing the impulse response function thus amounts to having full information about the behavior of the system, within its range of linearity, and the task of characterizing a linear system amounts to determining its impulse response. To do so by measurement, the output of the system must be observed for known input, requiring an appropriate measurement technique. The impulse response function can then be obtained by deconvolution of the measured output by the input, which can be efficiently performed by division in the frequency domain. Let ^ OðvÞ ¼ OðvÞ þ hðvÞ denote the observed output, including the measurement noise g(x) of standard deviation r(x). Then, the resulting estimate of HðvÞ is ^ ^ ðvÞ ¼ OðvÞ ¼ OðvÞ þ hðvÞ ¼ H ðvÞ þ hðvÞ H I ðvÞ I ðvÞ I ðvÞ ½2 with the signal-to-noise (SNR) ratio SNRH^ ðvÞ ¼ jHðvÞj  jIðvÞj : sðvÞ ½3 I(v) should therefore be large and ideally contain equal amplitudes of all frequencies in the bandwidth of interest. A perfect Dirac delta function would evenly cover an infinite bandwidth but is not a feasible input for a Gradient System Characterization With a Dynamic Field Camera 585 real-world gradient system. In modern MR systems, the gradient input is typically defined digitally with a finite dwell time and is subject to limitations in terms of the maximum amplitude and slope of the gradient waveform. Under these constraints, a feasible approximation of a Dirac delta function is a short triangle pulse. A triangular function can be described as the convolution of two identical boxcar functions: iðtÞ ¼ p  ðboxT  boxT ÞðtÞ ¼
p  ðT
jtjÞ; jtj  T 0; jtj > T ½4 where boxT has length T and amplitude 1, and p is the slope of the triangle (Fig. 1). In the frequency domain, a triangular input thus corresponds to a squared sinc function, IðvÞ / p  sin2 ðvT=2Þ ; v2 ½5 which has several relevant properties with respect to the choice of its parameters. First, as can be seen from Eq. 5, the envelope of the sinc-squared is determined by the slope, p, of the triangle. To optimize the sensitivity of the measurement, the slope should therefore generally be maximized. Second, the amplitude of the main lobe of I(v) scales with the total moment of the pulse: Zþ1 Ið0Þ / iðtÞdt ¼ p  T 2 ; ½6
1 and third, the sinc-squared input exhibits zeros with a spacing that scales inversely with the length of the triangle: vjI¼0 ¼ 2pn ; n ¼ 6 1; 2; 3:::: T ½7 These zeros amount to blind spots in the measurement. Therefore, in the choice of the triangle length, a trade-off needs to be made between the bandwidth covered by the main lobe, which is larger for a shorter pulse, and the sensitivity at low frequencies, which increases strongly with the pulse duration. This dilemma can be circumvented by using multiple triangular pulses of different length, giving complementary response information. When using multiple input functions, the data obtained in the different measure^ ments need to be combined in calculating HðvÞ. A leastsquares estimate at each frequency is obtained by ^ HðvÞ ¼  j Ij ðvÞ P P j ^ j ðvÞ O jIj ðvÞj 2 ¼ HðvÞ þ  j Ij ðvÞ P P j  hj ðvÞ jIj ðvÞj2 ; ½8 where I*(v) denotes the complex conjugate of I(v), and the index j counts the different input functions. Assuming equal noise levels and no noise correlation between the different measurements, the variance of the noise term in Eq. 8 sums up to s2 ðvÞ s2H^ ðvÞ ¼ P 2; j jIj ðvÞj ½9 FIG. 1. The set of triangular input functions used for GIRF measurements, shown in the time-domain (a), and corresponding sincsquared functions in the frequency-domain (b). The frequency-domain plot includes the envelope of the side-lobes falling off with p/v2, and the root-sum-of-squares of all inputs, indicating the combined sensitivity of the set of input functions. A zoom of the side-lobes is plotted in (c), showing how the different input functions complement each other in terms of frequency content. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] yielding the net SNR jHðvÞj  SNRH^ ðvÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 j jIj ðvÞj sðvÞ : ½10 As seen from Eq. 10, the SNR of the combined measurement is governed by the spectral root-sum-of-squares profile of the set of input functions and by the basic noise 586 Vannesjo et al. level of the individual measurements. Based on this relationship, the specific choice of input waveforms can be made according to the desired sensitivity profile. To boost net sensitivity, certain or all waveforms may be performed multiple times, yielding an additional SNR gain by implicit averaging through the summation in Eq. 8. As described so far, the proposed approach is geared to determining the self-response of any actuated field component, which may be a gradient field, a shim field, or else. However, it can readily be extended to incorporate crossresponses of undesired spatial field patterns as well. To this end, the magnetic field response in the magnet bore is expanded into a suitable set of spatial basis functions, which each constitute a separate output channel. Conversely, each control input of the system forms a separate input channel. For a gradient system, three input channels are considered, one each for the x, y, and z gradients. On this basis, individual impulse response functions are introduced for each pair of input channel l and output channel m, yielding the complete description om ðtÞ ¼ ¼ 1 X Z l X FT il ðtÞ  hl;m ðt
tÞdt $ Om ðvÞ
1 I ðvÞ l l  Hl;m ðvÞ: ½11 Given the capability to record the different outputs simultaneously, the multiple-channel impulse response can be determined simply by repeating the response measurements for each input channel. The resulting ^ ðlÞ ðvÞ yield observations O m;j ^ l;m ðvÞ ¼ H  j Il;j ðvÞ P P j ^ ðlÞ ðvÞ O m;j jIl;j ðvÞj2 ½12 with SNR as derived in Eq. 10. Dynamic Field Camera Measurements of magnetic field responses were performed with a third-order dynamic field camera (27), using the field monitoring methodology described in (29). In short, NMR probes (28,32) are used to sense the magnetic field dynamics at a set of different positions with high temporal resolution. The actual observable of each probe is the phase of its NMR signal, which is proportional to the integral of the magnetic field magnitude at the probe position. Based on knowledge of the probe positions, the spatiotemporal field evolution is then expressed in terms of the previously chosen spatial basis functions. The resulting time-varying coefficients, k^m ðtÞ, indicate the observed phase accumulated due to fields of the respective spatial structure. The timederivative of the phase coefficients gives corresponding field coefficients, which constitute the observed multiplechannel field output according to Eq. 11: ^m ðtÞ ¼ o ^m ðtÞ 1 dk ; g dt ½13 where g is the gyromagnetic ratio of the NMR-active nucleus in the probes. To fully determine the chosen spatial Table 1 Real-Valued Spherical Harmonics Used as Spatial Basis Functions Nr Spherical harmonic 0: 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 1 x y z xy zy 2z2
(x2 þ y2) xz x2
y2 3yx2
y3 xyz 5yz2
y(x2 þ y2 þ z2) 2z3
3z(x2 þ y2) 5xz2
x(x2 þ y2 þ z2) z(x2
y2) x3
3xy2 Order 0 1 2 3 model, the number of field probes must be equal to or larger than that of basis functions and the probes must be placed appropriately to permit a well-conditioned fit. The temporal resolution of the measurement is limited only by the acquisition bandwidth of the spectrometer used and the length of individual field observations is limited by the signal lifetime of the NMR probes. The dynamic field camera consisted of 16 unshielded transmit-receive NMR probes based on water samples of 0.8 mm in diameter, permitting field observations up to k values of  11,000 rad/m (28,29). The water was doped with 3.3 g/L of CuSO4 to reduce T1 for fast re-excitation, while keeping T2 long enough to maintain a signal lifetime of about 100 ms. The probes were evenly distributed on the surface of a sphere of 20 cm in diameter. Their positions were determined by measuring NMR frequency shifts under static gradients of 2.5 mT/m in the x, y, and z directions, respectively. GIRF Measurements In this work, real-valued spherical harmonics up to full third order were used to represent phase and field responses (Table 1). In total, this expansion comprises 16 basis functions, matching the number of available field probes. The single zeroth-order function captures the uniform component of phase and field, whereas the three first-order harmonics represent linear gradients in the x, y, and z directions. For better readability, the first^1 ; k^2 ; k^3 ) and field (^ ^2 ; o ^3 ) responses will o1 ; o order phase (k be equivalently referred to by the common notations kx ; ky ; kz and Gx ; Gy ; Gz . Measurements were performed on a 3 T Philips Achieva whole-body system with a maximum gradient strength of 40 mT/m and a slew rate limit of 200 mT/(mms). Gradient responses were acquired with and without built-in eddycurrent compensation and cross-term correction. The dwell time for discrete definition of the gradient waveforms was 6.4 ms. To avoid significant discretization errors, the shortest triangular gradient pulse used had a time-to-peak of 50 ms. Such a pulse has the first spectral zero at 20 kHz, Gradient System Characterization With a Dynamic Field Camera 587 FIG. 2. Magnitude (a, b, c) and phase (d, e, f) of measured GIRFs for all three gradient directions, with and without built-in eddy current compensation (ECC). The low-pass characteristic of the gradient system is apparent in (a). Details of the responses in (b) show how the eddy current compensation serves to broaden the response plateau and align the responses of the three different gradient channels. In a further zoom of the magnitude (c) and phase (f) plots, several mechanical resonances are apparent at low frequencies. Flat phase at low frequencies reflects near-zero net delay. where the system response is still non-negligible. Therefore, a combination of different pulse lengths was necessary. Overall, 12 triangular waveforms were used with a slope of p ¼ 180 mT/(m ms) and time-to-peak T ranging between 50 and 160 ms at increments of 10 ms (Fig. 1). The slope was chosen large for high sensitivity while keeping it well within the system’s specified range. The moment of the longest pulse was around 1200 rad/m and thus well within the observable range of the NMR field probes. In each of the response measurements, signal acquisition from the field probes was started 5 ms before the nominal start of the gradient pulse and had a duration of 68 ms, yielding a frequency resolution of about 14.7 Hz. The primary sampling bandwidth was 398 kHz and thus amply above the expected response bandwidth. The repetition time of the acquisitions was 2 s to allow for field settling between the individual measurements. Each of the 12 triangular inputs was performed 50 times per gradient channel. This results in a measurement time of 20 min for one gradient channel and a total of 1 h for all three gradient channels. From each individual set of probe signals, field outputs were calculated according to Eq. 13. All field outputs were then combined to yield the GIRF using Eq. 12. The whole procedure was performed separately for the x-, y-, and z-gradient channels. To verify the measured GIRFs, they were used to predict field responses to select gradient inputs according to Eq. 11. The predictions were then compared with direct measurements of the responses to the same gradient inputs. The gradient waveforms used for validation were a set of triangular pulses of 0.05 to 3.2 ms time-to-peak and 1 to 31 mT/m amplitude, and an echo-planar imaging (EPI) readout sequence (45 echoes, resolution ¼ 2.5 mm, acquisition time ¼ 62.5 ms). To actually probe the LTI assumption and accuracy of the GIRFs, the lengths of the triangular validation pulses were different from and partly much larger than those used for the GIRF measurement. The EPI was included to confirm that the linear model holds also for a full gradient sequence including both triangular and trapezoidal pulse shapes. RESULTS GIRF Measurements Magnitude and phase plots of the measured GIRFs in the frequency domain are shown in Fig. 2. By design, the gradient chains exhibit overall low-pass characteristics. The built-in eddy-current compensation broadens the response plateaus at low frequencies and aligns the responses of the three gradient subsystems. With compensation, the bandwidth at
3 dB is around 4.6 kHz, with the response tapering off toward zero at around 22– 25 kHz. As can be seen in the plot, the noise in the GIRF measurement increases considerably toward higher frequencies, where the power of the input waveforms is low. The phase response of the x and y gradients around DC (zero frequency) is essentially flat up to  400 Hz, with eddy current compensation. The corresponding zgradient phase response exhibits a minor slope at low frequencies indicating a residual group delay of
0.5 ms close to DC. The virtual absence of group delay in all three channels reflects appropriate delay calibration, which consists in constant shifts of the nominal time of the gradient chains relative to that of the acquisition subsystem. 588 FIG. 3. Time-domain representations of the GIRFs in all three gradient directions, with and without eddy current compensation (ECC). Vertical shifts (D) serve for better comparison of the GIRFs of the different gradient channels. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] At low frequencies, the GIRFs show channel-specific patterns of distinct peaks with amplitudes up to  1% of the full response. The x and y gradients exhibit several such peaks between 600 and 1800 Hz, while the GIRF of the z gradient shows mainly one rather narrow peak at 1300 Hz (Fig. 2). Based on their frequency range, amplitudes, and widths, these peaks are ascribed to mechanical resonances of the gradient coils. This interpretation is supported by a preceding acoustic observation of the single resonance at 1300 Hz (27). Acoustic responses at similar frequencies have also previously been noted in the literature (4,17). The prominent GIRF peaks at 25 kHz reflect the transistor switching frequency of the gradient amplifiers. The switching events cause slight modulation peaks at the switching frequency. However, these modulation peaks are not related to the actual input at 25 kHz and are consequently misinterpreted in the GIRF calculation. As the input at 25 kHz is extremely small, division according to Eq. 8 causes the small switching artifact to appear with large amplitude in the GIRF. Figure 3 shows plots of the three impulse responses, with and without eddy current compensation, in the time domain, filtered with a low-pass filter (raised cosine, 60 kHz full-width half-maximum (FWHM), 12 kHz transition band) to reduce high-frequency noise. The response kernels rise slightly more quickly than they fall and exhibit somewhat delayed second lobes, which are more pronounced for the y and z gradients. For the x and y gradients, the peak amplitude of the pre-emphasized responses is somewhat higher than without the eddy current compensation, whereas for the z gradient, the opposite is observed. From this representation, the width of the time-domain kernels was estimated at  40 ms with peaks at around
20 ms. The presence of a gradient response at t < 0 implies that a gradient field is observed before the nominal start of the gradient. This occurs as a result of the MR system’s aforementioned delay correction, which advances the nominal gradient time to play out gradient waveforms slightly sooner. Selected cross-responses are shown in Fig. 4, scaled to percent of the input. To compare field distributions of dif- Vannesjo et al. ferent spatial order, the basis functions were normalized to unit maximum field within a sphere of 20 cm diameter, centered at gradient isocenter (33). For instance, a linear field of 10 mT/m amounts to 1 mT in the normalized representation and is thus comparable with a zeroth-order field of 1 mT as they both produce the same maximum field shift within the defined sphere. The cross-responses were generally on the order of 0.1% or less for all three gradient chains, except at specific resonances and in the B0 term where it reached 1% for x-gradient input. In Fig. 4, several distinct cross-peaks can be seen at previously observed mechanical resonance frequencies of the respective gradient channel. In the x response to z-gradient input, there is a peak at the 1300 Hz resonance of the z gradient while the B0 and the higher-order 2z2
(x2 þ y2) (^ o6 ) terms both contain a resonance of the x gradient at o11 ) term shows 1780 Hz, and the 5yz2 – y(x2 þ y2 þ z2) (^ a resonance of the y gradient at 1700 Hz. The y response to x-gradient input shows peaks at 1060 Hz and 1750 Hz, both of which are less pronounced as resonances of the x gradient, the former however being a strong resonance in the self-response of the y gradient. Peaks at 1060 Hz appear in several higher-order field terms as well (data not shown), however only as responses to x-gradient inputs and not to the y gradient. A speculative explanation for this appearance of a resonance more strongly related to the y gradient than the driving x gradient may lie in the proximity of the two gradients and the shared structure embedding them both. Switching the x gradient may induce vibrations that propagate through the surrounding structure and decay more slowly in the y-gradient coil than in the x-gradient coil. The details, however, of how different vibrational modes are maintained, and how they affect the field inside the scanner remain to be investigated. The small distinct peaks seen at DC in several of the cross-terms largely stem from slow system drifts unrelated to gradient inputs as well as residual susceptibility broadening of the field probe samples (28), which mimics slight field variation at very low frequencies. GIRF-Based Predictions Figure 5 shows predictions of triangular gradient pulses with a time-to-peak T of 0.1, 0.4, and 1.05 ms, and peak amplitudes of 1 (Fig. 5a) and 31 mT/m (Fig. 5c), compared with nominal and directly measured (filtered to 60 kHz FWHM) pulses. Error plots of nominal vs. measured, and predicted vs. measured gradients are plotted for each pulse (Fig. 5b,d). For shorter pulses, the actual gradient deviates more from the nominal waveform than for longer ones, which is expected due to the weaker response at higher frequencies. The directly measured pulses follow the predictions closely with a maximal prediction error of about 0.2% of the pulse amplitude. Due to different dynamic ranges, the measurement noise appears more prominently in the upper plot. Increased noise toward the end of the 31 mT/m, 1.05 ms pulse, as seen in the error plot (Fig. 5d), was caused by nearly complete dephasing of the field probes in the validation measurement. The size of the NMR field probes thus limits the capability of validating predictions of large gradient pulses by single Gradient System Characterization With a Dynamic Field Camera 589 FIG. 4. Measured cross-responses of the x, y, z, B0, 2z2 - (x2 þ y2) and 5yz2 – y(x2 þ y2 þ z2) field terms due to operating the x- (blue), y- (green), and z-gradients (red), expressed in percent of the input. Several distinct peaks coincide with mechanical resonances previously identified in the self-responses (Fig. 2). Note the different scaling in the B0 plot. direct measurements. It does not, however, hinder GIRFbased predictions of large gradients per se. To study subtle long-term features of the gradient response, Fig. 6a shows a detail of the first 15 ms after a gradient pulse (30 mT/m amplitude, T ¼ 0.25 ms). In the prediction, a distinct oscillatory pattern is seen long after the pulse, reflecting the subtle resonance peaks in the underlying GIRF. The direct measurement, however, was too noisy to confirm the oscillations on this small scale. To see them more clearly, Fig. 6b shows the phase accumulated under the same gradient pulse. The underlying integration suppresses high-frequency noise and thus reveals the oscillations also in the direct measurement. It is now apparent that the predicted and the measured patterns again agree closely. The accumulated phase, expressed as the k-space position, is also the key parameter underlying MR image encoding and reconstruction. In Fig. 7, the measured GIRFs are used to predict a whole EPI readout sequence in terms of its effective gradient waveforms (left) and kspace trajectory (right). The directly measured trajectory again shows good accordance with the prediction, following it closely throughout, including the turns where it deviates significantly from the nominal trajectory. DISCUSSION The impulse response function of a linear time-invariant system is a well known and much used concept in engineering. Herein, we have investigated the feasibility of determining the impulse response functions of the gradient chains of an MR system. The investigation concerned each gradient chain as a whole between the two endpoints that are relevant for the user, i.e., between gradient waveform definition at the console and the actual magnetic fields seen by an object placed in the scanner bore. Using a third-order dynamic field camera based on NMR probes and a set of input functions that are easily implemented in standard scanner software, the impulse response measurement has proven a straightforward and fast procedure, yielding highly accurate results. When defined as above, GIRFs constitute a composite response reflecting all linear distortions that a gradient waveform undergoes between the software interface of the scanner and the resulting magnetic fields. This includes software signal processing, eddy-current compensation, gradient amplifier characteristics, potential cable effects, coil characteristics, coil coupling, eddy currents, and mechanical responses of the gradient system. 590 Vannesjo et al. stage of a gradient chain acts on the waveform, the input and output of that stage would need to be measured. Measuring the combined response, as done here, is a fast and comprehensive way of characterizing how a gradient system performs as one entity. Applications Knowledge of the gradient impulse response can be beneficial in a variety of ways. Most immediately, it forms a comprehensive basis of tuning the gradient response by pre-emphasis. Existing methods of eddy current compensation commonly rely on a limited number of exponentially decaying correction terms found by time-domain fitting. With the GIRFs, corresponding fits can be based on a complete picture of the native response in the frequency domain (34), including correction for mechanically induced field oscillations. Knowledge of the GIRF also permits preemphasis by a general filter equal to the frequency-domain ratio of the desired system response and the GIRF. In this fashion, LTI theory can be deployed to shape the net system response comprehensively within the technical limitations of the available gradient hardware. More generally, fast GIRF measurement may be useful for the calibration and quality assurance of gradient systems and as a diagnostic tool in cases of suboptimal gradient performance. As shown in this work, known GIRFs can also be used to calculate the actual gradient output that is obtained with any given pulse sequence. Given sufficient system stability, image reconstruction may thus be based on actual k-space trajectories and higher-order encoding (33) without the need for sequence-by-sequence field monitoring. Calculated field evolutions may also yield more accurate accounts of specialized signal preparation such as, e.g., in phase-contrast methods (35). In some situations, retrospective correction for field imperfections is not or hardly feasible. This is true, for instance, for spatially selective excitation based on k-space trajectories played out during radiofrequency transmission. In these cases, GIRF-based predictions could serve as a basis for designing trajectories such that they can be realized perfectly with a suitable precompensated demand. Model Limitations FIG. 5. GIRF-based predictions of gradient blips of 1 mT/m peak amplitude, T ¼ 0.1, 0.4, and 1.05 ms (a) and 31 mT/m peak amplitude, T ¼ 0.4 and 1.05 ms (c), compared with nominal and directly measured gradient waveforms. For each pulse, the differences between nominal and measured (black) as well as predicted and measured (blue) traces are plotted in (b, d). The prediction error amounts to maximally 0.2% of the gradient amplitude. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] Several features of the measured GIRFs have been linked to underlying phenomena. To fully understand how each The LTI approach is only valid to the extent that the system considered is actually LTI. Although typically minor, there are a number of influences on gradient responses that are nonlinear or vary over time. Power amplifiers generally exhibit nonlinear behavior to varying extent and procedures to linearize the output are a standard part of gradient amplifier engineering. Gradient responses also depend somewhat on the thermal state of the system. High-gradient-duty-cycle scans, such as for fMRI or diffusion imaging, heat up the gradient coils and supporting structures. Ensuing changes in electrical and mechanical properties will cause the gradient response to change slightly over time. Likewise, any other structures that support eddy currents may slightly heat up and thus change their eddy current characteristics. In as far as the heating is caused by driving the gradients themselves, it could be seen as an extension to the Gradient System Characterization With a Dynamic Field Camera 591 FIG. 6. Nominal, measured, and predicted waveforms of a triangular gradient pulse (30 mT/m peak amplitude, T ¼ 0.25 ms) in field (top) and accumulated-phase (bottom) representation. The zoomed insets on the right show details of the system response over about 15 ms after the pulse. As seen there, the GIRF predicts subtle oscillations after the pulse, which cannot be distinguished from noise in the direct field measurement. They are largely confirmed, however, by plots of the corresponding phase coefficient, ky , in which high-frequency noise is attenuated by integration over time. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] response of the system, which will however generally be nonlinear. There may also be external factors that influence the net field independently of how the system is driven and are thus inherently nonlinear. Any perturbation from the environment, such as fields produced by other electrical equipment, nearby traffic or power lines fall into this class. Nonetheless, as long as these model violations are suitably small, the LTI approach permits useful approximations of the system behavior. In this work, the validity of measured GIRFs and thus of the underlying LTI picture has been demonstrated by predictions of effective gradient waveforms. For a variety of gradient inputs, a high degree of prediction accuracy was achieved with a maximal error on the order of 0.2% of the pulse amplitude. Moreover, preliminary data indicate good stability of the system response over years. mechanical displacements in the system translate into field changes. Such studies are beyond the scope of this work. However, the observed agreements in vibration peaks between self-, cross-, and higher-order responses indicate that GIRF measurements will be instrumental in elucidating these mechanisms. Cross-responses via eddy currents are pronounced with unshielded gradient or shim coils, which couple strongly to the cryostat. In this situation, measurements of the cross-responses could form the basis of cross pre-emphasis and higher-order image reconstruction (33). In this work, a field camera with 16 NMR probes permitted GIRF measurements of full third order in terms of spherical harmonics. Further discrimination of field responses up to arbitrary spatial order could be readily achieved provided a suitable number and distribution of field probes. Mechanical Oscillations and Cross-terms Sensitivity and Resolution Among the features that the measured GIRFs capture and predict well are field oscillations due to mechanical vibrations. Judging by the achieved prediction accuracy, such oscillations seem to behave largely linearly. Understanding these effects in more detail will require models of how The sensitivity of the GIRF measurement has been shown to depend on the set of input functions used and on the noise level of individual field measurements. Triangular pulses offer the highest spectral sensitivity close to DC. Toward higher frequencies, their net sensitivity 592 Vannesjo et al. FIG. 7. Nominal, measured, and predicted gradient time-courses during an EPI readout sequence (left). Zoomed details show one of the phase-encoding blips and the settling of a frequencyencoding gradient (bottom left). The resulting full k-space trajectories are shown on the right, including a detail of one of the kspace corners including the end of the prephaser gradient. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] drops approximately as 1=v2 as seen from Eq. 5. The transition from phase data to field data further emphasizes high-frequency noise since taking the temporal derivative corresponds to multiplication by x in the frequency domain. Jointly, these effects caused the GIRF sensitivity to drop as 1=v3 overall. A straightforward way of increasing the overall SNR is to acquire averages of the field response measurements. To specifically increase sensitivity at higher frequencies, sinuisodally modulated pulses could be used in principle, shifting the spectral main lobe away from DC. To keep with slew rate limitations, however, the slope of a pulse to be modulated with a sinusoid would need to be kept lower than the system maximum, thereby reducing the total power of the pulse. For this reason, the sensitivity of a single measurement at a high frequency will necessarily be lower than what can be achieved around DC with the corresponding unmodulated pulse. The optimal choice of input pulses ultimately depends on the desired sensitivity profile for the system to be measured. Equations 8 and 12 describe how to obtain a least-squares estimate of the GIRFs based on the sensitivity profile of each input pulse, and hold for any combination of modulated and unmodulated pulses alike. Combining simple triangular pulses as done here offers the advantage of straightforward implementation on standard MR equipment. In this work, the gradient response measurements covered 68 ms, yielding a frequency resolution of 14.7 Hz. Particularly for unshielded gradient coils, it is known that significant eddy currents may persist for several hundreds of milliseconds (9), thus exceeding the range possible to capture in one single measurement with the field camera used in this work. To overcome this, multiple response measurements with interleaved timing can be acquired (34), which however comes at the expense of increasing the total time required for the measurements. With hardware able to perform continuous field monitoring, as described by (36), long-living responses could be captured straightforwardly. Measurement Time The time required for measuring GIRFs with the proposed approach depends on the total number of input pulses and the time allowed for gradient settling between pulses. The aim of the measurements presented here was to elucidate even subtle behavior of the gradient responses, to which end a large number of averages was acquired. For three gradient channels, the total measurement time amounted to  1 h. For a coarser measurement, even one observation of each pulse would suffice, reducing the measurement time to about 1 min. Initial experience suggests that even such coarser GIRFs yield useful predictions of gradients and trajectories. This is largely due to the fact that common input gradient waveforms have most power at relatively low frequencies, where the SNR of the measured GIRF is favorable. CONCLUSION GIRF measurements with a third-order dynamic field camera enable fast and easy determination of gradient system characteristics. Among others, GIRFs measured in Gradient System Characterization With a Dynamic Field Camera 593 this way readily reveal subtle field oscillations due to mechanical resonances. These have been found to be largely linear and time-independent in experiments without significant gradient heating. GIRFs also permit highly accurate predictions of net system behavior. Such predictions will be useful for image reconstruction and gradient waveform optimization. Notably, suitable sets of field probes straightforwardly yield GIRFs up to high spatial order along with the common gradient responses. This capability renders the proposed approach promising also for the characterization of dynamic shim systems. The GIRF concept assumes linearity and time-independence of the system to be characterized but is otherwise free of mechanistic assumptions. For more advanced models incorporating, e.g., thermal or nonlinear electromechanical effects, response measurements with a dynamic field camera will still be instrumental. 14. Barnett A. Comments on ‘‘Gradient-induced acoustic and magnetic field fluctuations in a 4 T whole-body MR imager.’’ Magn Reson Med 2001;46:207. 15. Nixon TW, McIntyre S, Rothman DL, de Graaf RA. Compensation of gradient-induced magnetic field perturbations. J Magn Reson 2008; 192:209–217. 16. Hedeen RA, Edelstein WA. Characterization and prediction of gradient acoustic noise in MR imagers. Magn Reson Med 1997;37:7–10. 17. Li W, Mechefske CK, Gazdzinski C, Rutt BK. Acoustic noise analysis and prediction in a 4-T MRI scanner. Concept Magn Res 2004;21B: 19–25. 18. Sierra CVR, Versluis MJ, Hoogduin JM, Duifhuis H. Acoustic fMRI noise: linear time-invariant system model. IEEE Trans Biomed Eng 2008;55:2115–2123. 19. Addy NO, Wu HH, Nishimura DG. Simple Method for MR Gradient System Characterization. In Proceedings of the 17th Annual Meeting of ISMRM, Honolulu, Hawaii, USA, 2009. p. 3068. 20. Vannesj€ o SJ, H€ aberlin M, Kasper L, Barmet C, Pruessmann KP. A Method for Characterizing the Magnetic Field Response of a Gradient System. In Proceedings of the 18th Annual Meeting of ISMRM, Stockholm, Sweden, 2010. p. 1536. 21. Wysong RE, Lowe IJ. A simple method of measuring gradient induced eddy currents to set compensation networks. Magn Reson Med 1993;29:119–121. 22. Terpstra M, Andersen PM, Gruetter R. Localized eddy current compensation using quantitative field mapping. J Magn Reson 1998;131: 139–143. 23. Duyn JH, Yang YH, Frank JA, van der Veen JW. Simple correction method for k-space trajectory deviations in MRI. J Magn Reson 1998; 132:150–153. 24. Alley MT, Glover GH, Pelc NJ. Gradient characterization using a Fourier-transform technique. Magn Reson Med 1998;39:581–587. 25. Goodyear DJ, Shea M, Beyea SD, Shah NJ, Balcom BJ. Single point measurements of magnetic field gradient waveform. J Magn Reson 2003;163:1–7. 26. Balcom BJ, Bogdan M, Armstrong RL. Single-point imaging of gradient rise, stabilization, and decay. J Magn Reson Ser A 1996;118: 122–125. 27. Barmet C, Wilm BJ, Pavan M, Pruessmann KP. A Third-Order Field Camera with Microsecond Resolution for MR System Diagnostics. In Proceedings of the 17th Annual Meeting of ISMRM, Honolulu, Hawaii, USA, 2009. p. 781. 28. De Zanche N, Barmet C, Nordmeyer-Massner JA, Pruessmann KP. NMR probes for measuring magnetic fields and field dynamics in MR systems. Magn Reson Med 2008;60:176–186. 29. Barmet C, De Zanche N, Pruessmann KP. Spatiotemporal magnetic field monitoring for MR. Magn Reson Med 2008;60:187–197. 30. Barmet C, Wilm BJ, Pavan M, Katsikatsos G, Keupp J, Mens G, Pruessmann KP. Concurrent Higher-Order Field Monitoring for Routine Head MRI: An Integrated Heteronuclear Setup. In Proceedings of the 18th Annual Meeting of ISMRM, Stockholm, Sweden, 2010. p. 216. 31. Lathi BP. Signal processing and linear systems, Indian Edition. New York: Oxford University Press; 2008. 32. Barmet C, De Zanche N, Wilm BJ, Pruessmann KP. A transmit/ receive system for magnetic field monitoring of in vivo MRI. Magn Reson Med 2009;62:269–276. 33. Wilm BJ, Barmet C, Pavan M, Pruessmann KP. Higher order reconstruction for MRI in the presence of spatiotemporal field perturbations. Magn Reson Med 2011;65:1690–1701. 34. Vannesj€ o J, Fillmer A, Barmet C, Boesiger P, Henning A, Pruessmann KP. Fast Characterization of Higher-Order Shim Dynamics by Impulse Response Measurements with a Dynamic Field Camera. In Proceedings of the 19th Annual Meeting of ISMRM, Montreal, Canada, 2011. p. 719. 35. Giese D, Haeberlin M, Barmet C, Pruessmann KP, Schaeffter T, Kozerke S. Analysis and correction of background velocity offsets in phase-contrast flow measurements using magnetic field monitoring. Magn Reson Med 2012;67:1294–1302. 36. Dietrich BE, Barmet C, Brunner D, Pruessmann KP. An Autonomous System for Continuous Field Monitoring with Interleaved Probe Sets. In Proceedings of the 19th Annual Meeting of ISMRM, Montreal, Canada, 2011. p. 1842. ACKNOWLEDGMENTS The authors thank Geran Peeren at Philips Healthcare, as well as Frederik Testud at the University of Freiburg, for reading and providing valuable comments on this work. Technical support from Philips Healthcare is gratefully acknowledged. REFERENCES 1. Boesch C, Gruetter R, Martin E. Temporal and spatial analysis of fields generated by eddy currents in superconducting magnets: optimization of corrections and quantitative characterization of magnet/ gradient systems. Magn Reson Med 1991;20:268–284. 2. Liu Q, Hughes DG, Allen PS. Quantitative characterization of the eddy current fields in a 40-cm bore superconducting magnet. Magn Reson Med 1994;31:73–76. 3. Wu YH, Chronik BA, Bowen C, Mechefske CK, Rutt BK. Gradientinduced acoustic and magnetic field fluctuations in a 4 T wholebody MR imager. Magn Reson Med 2000;44:532–536. 4. Clayton DB, Elliott MA, Leigh JS, Lenkinski RE. 1H Spectroscopy without solvent suppression: characterization of signal modulations at short echo times. J Magn Reson 2001;153:203–209. 5. Foerster BU, Tomasi D, Caparelli EC. Magnetic field shift due to mechanical vibration in functional magnetic resonance imaging. Magn Reson Med 2005;54:1261–1267. 6. Van de Moortele P, Pfeuffer J, Glover GH, Ugurbil K, Hu X. Respiration-induced B0 fluctuations and their spatial distribution in the human brain at 7 Tesla. Magn Reson Med 2002;47:888–895. 7. Versluis MJ, Peeters JM, van Rooden S, van der Grond J, van Buchem MA, Webb AG, van Osch MJP. Origin and reduction of motion and f0 artifacts in high resolution T2*-weighted magnetic resonance imaging: application in Alzheimer’s disease patients. NeuroImage 2010;51:1082–1088. 8. Jehenson P, Westphal M, Schuff N. Analytical method for the compensation of eddy-current effects induced by pulsed magnetic field gradients in NMR systems. J Magn Reson 1990;90:264–278. 9. van Vaals JJ, Bergman AH. Optimization of eddy-current compensation. J Magn Reson 1990;90:52–70. 10. Wysong RE, Madio DP, Lowe IJ. A novel eddy current compensation scheme for pulsed gradient systems. Magn Reson Med 1994;31:572–575. 11. Jezzard P, Barnett AS, Pierpaoli C. Characterization of and correction for eddy current artifacts in echo planar diffusion imaging. Magn Reson Med 1998;39:801–812. 12. Zur Y, Stokar S. An algorithm for eddy currents symmetrization and compensation. Magn Reson Med 1996;35:252–260. 13. Eccles CD, Crozier S, Westphal M, Doddrell DM. Temporal spherical-harmonic expansion and compensation of eddy-current fields produced by gradient pulses. J Magn Reson Ser A 1993;103: 135–141.