1 Introduction
Magnetic resonance imaging (MRI) is a widely used imaging modality in clinical practice due to its ability to produce a good contrast between soft tissues and to image a slice at any orientation. The MRI scanner acquires data in a Fourier domain called kspace. MRI captures an image by scanning through the kspace on a Cartesian or nonCartesian trajectory. This scanning process is time consuming and results in a long acquisition time and patient discomfort. Accelerating the data acquisition process is an active area of research in MRI, and Compressive Sensing (CS) is a promising solution that can improve the speed of data acquisition in MRI. CS [1, 2, 3, 4, 5] is a technique that permits the faithful reconstruction of the signal of interest from the data acquired below the Nyquist sampling rate. MRI is a ideal system [6] for CS applications as it acquires image already in encoded form rather than in pixel domain. The application of CS in MRI was first described in [7], where variable density random undersampling of phase encodes was suggested as a sampling strategy. Parallel imaging techniques [8, 9, 10] have also been combined with CS in [11, 12, 13, 14, 15, 16] to further accelerate MRI scans, and CSMRI has been applied for dynamic imaging, exploiting kt space sparsity [14, 17, 18, 19, 20].
The theory of CS provides a solution to an illposed inverse problem by exploiting prior knowledge of signal sparsity or compressibility. This theory guarantees perfect reconstruction of the signal from the undersampled data if certain conditions are satisfied [1, 2, 3, 4, 5, 21]: (i) sparsity or compressibility of the signal in the transform domain; (ii) Restricted Isometry Property (RIP) of the measurement matrix or incoherence between the measurement and sparsifying transform matrices; and (iii) a nonlinear reconstruction algorithm that promotes sparse representation of the image and enforces data consistency of reconstruction with the acquired data.
The sparsity or compressibility condition is satisfied by MR images as they are known to be sparse or compressible in the wavelet domain and the finite difference domain [7, 11, 12, 13]. However, the RIP is difficult to verify for a given deterministic measurement matrix since it is computationally NP hard [22]. An empirical solution to this problem in the CS literature is to use random measurement matrices. A randomly sampled frequency domain data can capture pertinent information from a sparse signal with fewer measurements and allows accurate reconstruction of the signal by the convex optimization program. This property was first proved mathematically for Gaussian matrices [5, 23] and has recently been extended to a wide class of random matrices [24]. Based on this property, [25]
proposed using spatially selective RF pulses to implement random encoding along the phase encode direction, with the entries of the random measurement matrix drawn from Gaussian distribution. This random encoding scheme attempts to approximate the sufficient conditions for perfect CS reconstruction, but as described in
[25], this measurement matrix is not unitary and results in noise amplification even after taking all the required measurements. Another problem with random encoding is computational complexity. Dense random matrices consume large amounts of memory and require computationally expensive matrix multiplications in CSreconstruction [3, 26]. This problem is partially alleviated in [25]by using fast Fourier transforms of the matrix multiplications, but still requires more memory and computations than those of structured/unitary measurement matrices.
MRI uses the Fourier basis to encode the excited region of interest. The Fourier measurement matrix is weakly incoherent with the wavelet sparsifying transform matrix, thus is suboptimal for CSMRI [27]. The incoherence is essentially a measure of the spread of sparse signal energy in the measurement domain [3]. Various attempts have been made in [25, 27, 28, 29, 30] to spread the energy of the MR signal in the measurement domain. In [27, 31], the spread spectrum technique was presented which convolves the kspace with the Fourier transform of a chirp function to spread the energy of the MR signal in the measurement domain. The chirp modulation is implemented through the use of second order shim coils. In [28, 29, 30, 32, 33, 34], other nonFourier encoding strategies were described for compressive sensing that aims to spread the energy of the MR signal in the measurement domain. While these encoding strategies can spread the signal energy to some extent, none of them has the theoretically proven maximal incoherence for the complete spread of the signal energy.
Noiselet bases [35, 36] are known to completely spread out the energy of the signal in the measurement domain, which is a desired property in CS. Noiselets are also known to be maximally incoherent with Haar wavelets that makes them the best suited bases for CS. Further, noiselet matrices are complex valued, symmetric and unitary, which simplifies the implementation of image reconstruction program in CSMRI. In the simulation study of [37], it is found that the noiselet measurement matrix outperforms the chirp modulation measurement matrix when the noise level is high. Also, as shown in Section III of this paper, the multichannel noiselet measurement matrix exhibits much better RIP than that of its Fourier counterpart. In order to take the advantage of maximal incoherence and better RIP provided by noiselet measurement matrix, we have investigated the use of noiselet encoding for CSMRI.
In order to take the advantage of multiple measurements provided by an MR scanner through the use of multiple channels, a Multichannel Compressive Sensing (MCS) framework is proposed in [11] for CS reconstructions. The MCS framework simultaneously uses data from the multiple channels to reconstruct the desired image instead of reconstructing separate images from each channel, resulting in higher acceleration factors and improved image quality. Therefore, in this paper we describe the theory and implementation details of using noiselet bases as the measurement matrix in MCSMRI. Considering the lack of analysis and sufficient understanding of MCSMRI in the literature, we also present an empirical RIP analysis of the multichannel noiselet measurement matrix in comparison with its Fourier counterpart. The results indicate that the multichannel noiselet measurement matrix outperforms its Fourier counterpart, and that noiselet encoding outperforms Fourier encoding in preserving image resolution for the same acceleration factors, and can achieve higher acceleration factors than the Fourier encoding scheme for the desired image quality and resolution.
The paper is organized as follows. In section II, we describe the background of CS, sufficient conditions for CS and develops a model for MCSMRI reconstruction. In section III, we describe the noiselet basis function, its properties and our motivation for using noiselets in MRI. A pulse sequence design to implement the proposed noiselet encoding scheme is also described in this section, followed by an empirical RIP analysis of the multichannel noiselet measurement matrix in comparison with its Fourier counterpart. In section IV, simulation studies comparing the performance of noiselet encoded and Fourier encoded MCSMRI for different acceleration factors are demonstrated on a brain image. The effect of the number of channels and level of noise on the reconstruction is also evaluated for both the encoding schemes. In section V, we demonstrate the feasibility of the proposed encoding scheme by acquiring noiselet encoded data from a phantom and a human brain. Retrospective undersampling is performed on the acquired noiselet encoded and the Fourier encoded data to simulate accelerated acquisition. The nonlinear conjugate method [6] with wavelet and total variation (TV) penalties is used to solve the minimization program for MCSMRI. In section VI, we discuss the findings, limitation and further extension of the technique.
2 Compressive Sensing
Compressive sensing is a mathematical theory describing how a sparse signals can be faithfully recovered after sampling projections well below the Nyquist sampling rate. Consider a signal in the dimensional complex space that can be sparsely represented in domain as , where is the sparsifying transform matrix. The signal is sparse, that is, only the coefficients in are nonzero. A measurement system measures signal in dimensional space by taking only projections of the signal as
(1) 
where , and is the measurement matrix. In MRI, is usually a (partially) randomly undersampled discrete Fourier transform matrix. Equation (1) can be further expressed as
(2) 
where represents the conjugate transpose operation and the signal is sparse in the domain. MR images can be sparsely represented in the wavelet domain using the wavelet transform matrix. Given the measurement and the matrices and , there exist many solutions satisfying (1) and recovering becomes an ill posed problem. The CS theory provides a unique solution to the illposed problem by solving the following optimization program:
(3) 
where is the norm of , with the th element of . Exact reconstruction of the signal is achievable if certain mathematical conditions hold.
2.1 Restricted Isometry and Incoherence in Compressive Sensing
An important sufficient condition for exact reconstruction of is the so called restricted isometry property (RIP) [2, 21, 23]. For a normalized measurement matrix with unit column norms, the RIP is given as
(4) 
where is called the RIP constant. The RIP (4) is equivalent to [2, 38, 39]
(5) 
where is the submatrix formed from distinct columns of , and and
are the minimum and maximum singular values
^{1}^{1}1In [2, 38, 39], the eigenvalue
is used for (5), where . of , respectively. The RIP constant is the smallest constant that satisfies the inequalities (5) for every sized submatrix of , and it is essentially a bound on the distance between unity and the singular values of all s. It is shown in [23] that if , then a sparse signal can be exactly reconstructed from the measurements of .While renders the exact reconstruction of , the value of determines the stability of reconstruction. In the presence of measurement noise , and the reconstructed signal satisfies (Section 5.2 [40])
(6) 
Thus, the smaller the , the smaller the reconstruction error, and vice versa. Since measurement noise always exists in practice, the value of is an important performance measure of a measurement matrix for both the reconstructability and reconstruction error. However, the computation of for a given is NP hard and hence intractable. Since the RIP constant is essentially a bound on the distance between 1 and the singular values of all s, the value of can be assessed by the distances from 1 to the s and s. The smaller the distance, the smaller the and hence the better performance of . Since (5) must hold for all the submatrices of , the statistics of and over randomly sampled s are used in [5, 38, 39] to assess the RIP performance of a given measurement matrix . This method is also adopted in this paper.
Another important sufficient condition for exact reconstruction of is the incoherence [3, 36]. For a pair of measurement matrix and sparsifying transform matrix , satisfying and , their incoherence is defined as
(7) 
where and are respectively the th and th columns of and , and . The value is termed as maximal incoherence. As shown in [3], if , where is a small constant, then a sparse signal can be exactly reconstructed. Thus, determines the minimum number of measurements needed for exact reconstruction of . The smaller the , the smaller the (the fewer the measurements) needed for exact reconstruction of .
It is important to note that both the RIP and the incoherence are sufficient conditions on the measurement matrix. So they are parallel and either or both of them can be used to design, analyze and assess the measurement matrix for exact reconstruction of .
2.2 Multichannel Compressive Sensing
MRI systems acquire multiple measurements of the desired signal through multiple channels. Given the multiple channels, the data acquisition process can be modeled as where is the complex valued sensitivity map matrix of the th receive channel, with being the sensitivity of the th channel at the
th pixel of the vectorized image,
is the number of receive channels, is the data acquired from the th receive coil and . In matrix form, the ’s can be written as(8) 
As seen from above, with receive channels, the measurement matrix for becomes , which has a column of measurement matrices ’s and dimension . It is important to note that s are complex valued and , in general. Hence, for and they can be independent of each other, depending on the specific values of and . As a result, the multichannel measurement matrix provides more independent measurements than that of the single channel , which may reduce the number of measurements, , needed at each channel for exact reconstruction of . This reasoning is confirmed by the empirical analysis of in Section III.
In light of the above discussion, the following MCS optimization is considered for reconstructing the desired image from the multichannel measurements of MRI.
(9) 
where is the wavelet transform operator and determines the allowed noise level in the reconstructed image.
MR images are also known to be sparse in the total variation (TV) domain. It is demonstrated in [41] that the TV penalty is critical to the performance of CSMRI, and that MR images can be recovered more efficiently with the use of TV penalty together with the wavelet penalty. Therefore, most of the CSMRI work[6, 7, 12, 13, 25, 28, 30, 33, 34] has used both TV and wavelet penalties for better reconstruction performance. To be consistent with this common practice in CSMRI, the TV penalty is included in the objective function (9) for MCSMRI reconstruction, together with the wavelet penalty. Equation (9) is a constrained optimization problem which is computationally intensive to solve. To relax the problem, (9) is converted to the unconstrained optimization problem with the inclusion of TV penalty.
(10) 
where is a 2D total variation operator and , are regularization parameters for wavelet and penalties, respectively.
Daubechies4 (db4) wavelet is usually used in CSMRI because of its superior performance in sparsifying the MR images. To be consistent with this fact and fair in comparison with the existing CSMRI results, the unconstrained objective function (10) with the of the db4 wavelet operator will be used throughout all the simulations and reconstructions in this work.
3 Noiselet Encoding in CSMRI
3.1 Noiselets
Noiselets are functions which are noiselike in the sense that they are totally incompressible by orthogonal wavelet packet methods [35, 36]. Noiselet basis functions are constructed similar to the wavelet basis functions, through a multiscale iteration of the mother bases function but with a twist. As wavelets are constructed by translates and dilates of the mother wavelet function, noiselets are constructed by twisting the translates and dilates [3]. The mother bases function can be defined as
The family of noiselet basis functions are generated in the interval as
(11) 
where and form the unitary basis for the vector space . An example of a 44 noiselet transform matrix is given below.
(12) 
Noiselets totally spread out the signal energy in the measurement domain and are known to be maximally incoherent with the Haar wavelet. The mutual incoherence parameter between the noiselet measurement matrix and the sparsifying Haar wavelet transform matrix has been shown to be equal to 1 [3], which is the minimum value possible for the incoherence. Therefore, theoretically, noiselets are the best suited measurement basis function for CSMRI when the wavelet is used as sparsifying transform matrix.
3.2 Motivation
The motivations behind using noiselets as a measurement matrix in MCSMRI are as follows:

Noiselet basis function is unitary and hence does not amplify noise as in the case of random encoding [25].

Noiselets completely spreads out the signal energy in the measurement domain and are maximally incoherent with wavelets.

Unlike random basis, noiselet basis has conjugate symmetry. Thus, this property of symmetry can be exploited by using the partial Fourier like technique.

Noiselets are derived in the same way as wavelets, therefore it can be modelled as a multiscale filterbank and can be applied in .
We proposed to use the noiselet encoding in the phase encode (PE) direction in 2D and 3D MR imaging. Therefore, the acquired data is noiselet encoded in the PE direction and Fourier encoded in the frequency encode (FE) direction/s.
3.3 Pulse Sequence Design for Noiselet Encoding
In conventional 2D MR imaging sequences, a spatially selective RF excitation pulse is used to select the slice and the linear spatial gradients are used to encode the excited slice onto the Fourier transform space. In [42, 43, 31, 32], it is demonstrated that the spatially selective RF excitation pulse can also be used to encode the imaging volume. In [44, 32, 45, 25], the wavelet, SVD and random encoding profiles have been implemented using the spatially selective RF excitation pulses. An analysis using the linear response model described in [42] provides a theoretical framework to design spatially selective RF excitation pulses for implementation of nonFourier encoding. Under the small flip angle ( 30) regime, the RF pulse envelope can be calculated directly by taking the Fourier transform of the desired excitation profile. However this method of designing an RF excitation pulse requires excellent RF and main field homogeneity.
To excite a noiselet profile during excitation, one can design an RF pulse envelope by directly taking the Fourier transform of the noiselet basis functions. For an image of size 256256, the noiselet measurement matrix has 256 rows and 256 columns (see (12) for the low dimensional example). The Fourier transformation of each row of the noiselet matrix will result in 256 RF excitation pulses.
A pulse sequence for the noiselet encoding of 2D MR imaging is shown in Fig. 1 (a). The pulse sequence is designed by tailoring the spin echo sequence. The RF excitation pulse in the conventional spin echo sequence is replaced by the noiselet RF pulse, and the slice select gradient is shifted to phase encoding axis. The 180 refocusing RF pulse is used in conjunction with the slice selection gradient to select the slice that refocuses the spins only in the desired slice. Spoilers are used after the readout gradient to remove any residual signal in the transverse plane. A new RF excitation pulse is used for every new TR to excite a new noiselet profile, and a total of 256 TR are required for excitation of the complete set of noiselet bases. The readout gradient strength determines the FOV in the readout direction, while the phase encoding gradient strength and duration of the RF excitation pulse determines the field of view (FOV) in phase encoding direction. The FOV in phase encoding direction is determined as
(13) 
where is the gradient strength in PE direction, is the gyromagnetic ratio, and is the dwell time of the RF pulse which is defined as = (Duration of RF pulse) / (Number of points in RF pulse). Equation (13) is used to calculate the gradient strength required in the phase encoding direction during execution of RF excitation pulse.
The method described above can also be used to design the pulse sequence for the noiselet encoding of 3D MR imaging as shown in Fig. 1 (b).
3.4 Undersampling in noiselet encoding
Noiselet transform is a type of HaarWalsh transform. The noiselet transform coefficients totally spread out the signal in scale and time (or spatial location) [35]
. As a result, each subset of the transform coefficients contains a certain information of the original signal at all the scales and times (spatial locations), and can be used alone with zero padding to reconstruct the original signal at a lower resolution. This important property is demonstrated by the example shown in Fig.
2.Fig. 2 shows a brain image of size , and the 3D magnitude map of the noiselet transform of the brain image along the phase encoding direction (all noiselet encodes). Fig. 2 (cf) shows the images reconstructed with the first, second, third and fourth 64 noiselet encodes by zero padding the rest. Each of these images are reconstructed using one quarter of the noiselet encodes and has low resolution than the original image. However, each of these images have complementary information about the original image and have approximately the same amount of energy and information because they are reconstructed using the same size of partial matrix from the original coefficient matrix.
Based on the above property of noiselet transform, we propose to undersample the noiselet encoded data along the phase encoding direction according to the uniform probability distribution function. One sampling mask using this scheme is shown in Fig.
3(a) where the white lines represent the sampled data points and the black lines represent the unsampled data points. Fig. 3(b) shows the sampling mask for Fourier encoding scheme drawn from a variable density probability distribution function shown in Fig. 3(c).3.5 Empirical RIP analysis of measurement matrix
According to the CS theory summarized in Section II, the measurement matrix is crucial to the performance of CS reconstruction, and the performance of a measurement matrix for a sparse signal can be assessed by the statistics of and over the s consisting of distinct columns of . To understand the behavior and advantage of the noiselet encoding proposed above, we have used this method to assess the noiselet measurement matrix in comparison to the conventional Fourier measurement matrix.
In the assessment, the size of the signal was , the number of measurements , the number of channels and , and the sparsity was varied from 5 to 100 with an increment of 5. For , the measurement matrices, s, were generated for the noiselet and Fourier encodings, respectively. For and , the measurement matrices s as given in (8) were generated for the noiselet and Fourier encodings, respectively. For each , 2,000 submatrices s were drawn uniformly randomly from the columns of , then the and of every were calculated. The same procedure is used to obtain the submatrices s from and to calculate the and of every . The statistics of s, s, s and s were accumulated from their respective 2,000 samples.
Fig. 4
shows the means and standard deviations of the minimum and maximum singular values of
s and s versus the sparsity for the Fourier and noiselet measurement matrices. As seen from the figure, in single channel case, the singular values of noiselet measurement matrix are closer to 1 than those of Fourier measurement matrix, but are not significantly different. As the number of channels increases, the singular values of noiselet and Fourier measurement matrices all move towards 1, but those of noiselet measurement matrix move much closer to 1 than those of Fourier measurement matrix. By the CS theory, when the maximum distance from 1 to the singular values is less than 1, it equals roughly the RIP constant . Therefore, the figure actually reveals two facts: 1) For both the noiselet and Fourier measurement matrices, the RIP constant decreases as the number of channels increases. 2) As the number of channels increases, the RIP constant of noiselet measurement matrix deceases much more than that of Fourier measurement matrix. According to the CS theory, these imply that the multichannel measurement matrix should generally outperform the single channel measurement matrix, and that the multichannel noiselet measurement matrix should generally outperform the multichannel Fourier measurement matrix.As a particular example consider the curves in Fig. 4 (d) for the noiselet measurement matrix. To facilitate discussion, the distances from 1 to the singular values of a measurement matrix will be called the distances here. In single channel case, the distances of noiselet measurement matrix are less than 1 for . By RIP, this implies that the single channel noiselet measurement matrix can guarantee the recovery of the signals with sparsity . When the number of channels is increased to 14, the distances are less than 1 for . So the 14 channel noiselet measurement matrix can guarantee the recovery of the signals with sparsity . The improvement in terms of sparsity is two folds. In contrast, for the 14 channel Fourier measurement matrix shown in Fig. 4 (c), its distances are less than 1 only for , so it can only guarantee the recovery of the signals with sparsity .
From the above assessment, it can be expected that the multichannel CS MRI will outperform the single channel CS MRI and that the noiselet encoding multichannel CS MRI will outperform the Fourier encoding multichannel CS MRI in practice. These are confirmed by the simulation and experiment results presented in the next sections.
It is important to note that the above empirical analysis results are obtained by using the complex valued sensitivity maps. If only the magnitudes of the sensitivity maps are used, the above results will not hold. Therefore, the complex valued sensitivity maps of multiple receive coils are the source for the significant performance improvement of the multichannel measurement matrix.
4 Simulation Study and Results
Simulations were performed on a (256
256) brain image to investigate the performance of noiselet encoded and Fourier encoded MRI. The simulation study was divided into two parts: (i) a simulation study with a single channel using uniform sensitivity and (ii) a simulation study with multiple channels where the sensitivity profiles were estimated from the acquired data.
4.1 Single Channel Simulation with a Uniform Sensitivity Profile
Fourier encoded CSMRI A Fourier transform of the image was taken in the PE direction to simulate Fourier encoding. Two types of sampling strategies were used to sample Fourier encoded data: (i) a variable density random sampling pattern as shown in Fig. 3 (a) where samples were taken in the PE direction according to a Gaussian distribution function, and (ii) a completely random sampling pattern as shown in Fig. 3
(b) where samples were taken in the PE direction according to the uniform density function. The nonlinear program of (
10) was solved to reconstruct the final image for acceleration factors of 2 and 3. In these cases the encoding matrix does not have any sensitivity information (i.e. ).Noiselet encoded CSMRI A noiselet transform of the image was taken in the PE direction to simulate noiselet encoding. A completely random sampling pattern was used to sample the noiselet encoded data in the PE direction and the nonlinear program of (10) was solved to reconstruct the final image for acceleration factors of 2 and 3. In these cases the encoding matrix does not have any sensitivity information (i.e. ).
Fig. 5 show the images reconstructed with Fourier encoding and noiselet encoding using variable density random undersampling and completely random undersampling pattern respectively. The noiselet encoded CSMRI performs similar to that of the Fourier encoded CSMRI. This is due to the fact that in the case of variable density random undersampling, the Fourier encoding judiciously exploits extra information about the data, namely the structure of kspace. The center of the kspace data has maximum energy and hence, by densely sampling the center of kspace, the Fourier encoding captures most of the signal energy and results in better performance.
In practice, the MR data is collected through the use of multiple channels, and data in each channel is slightly different from the other channels. The actual kspace data is convolved with the Fourier transform of the sensitivity profiles of the individual channel, making the data from each channel different from others. This sensitivity information can also be taken into consideration while performing the CS reconstruction, by applying the multichannel CS frame. Therefore, to further study the effect of sensitivity information on noiselet encoding and Fourier encoding, MCSMRI simulations were performed. To quantitatively compare the performance of both the encoding schemes, we used the relative error defined in (14) as a metric:
(14) 
First we investigated the effect of the number of channels on the reconstruction quality using the MCS framework. For a fixed number of measurements, the number of channels was varied and the mean of the relative error for 1000 such simulations was calculated. Fig. 6 shows the plot of the mean relative error versus the number of channels for the acceleration factors of 2 and 3. When the number of channels was two, the noiselet encoding scheme outperformed the Fourier encoding scheme for both the acceleration factors of 2 and 3. However, when number of channels was equal to one, the noiselet encoding outperformed the Fourier encoding for an acceleration factor of 2, but not for an acceleration factor of 3. It is interesting to note that noiselet encoding outperformed Fourier encoding for both acceleration factors when the number of channels was greater than one. These simulations suggest that noiselet encoding should take into account the sensitivity information while performing CS, and therefore noiselet encoding is potentially a better encoding scheme for MCSMRI. Based on the fact that noiselet encoding performs better than Fourier encoding for multichannel data, we investigated the performance of both the encoding schemes using mutichannel data.
4.2 Multiple Channel Simulation
A (256256) brain image was used to compare the performance of noiselet encoding and Fourier encoding in MCSMRI for different acceleration factors. Eight complex sensitivity maps (Fig. 7) obtained from the head coil of a Siemens Skyra 3T scanner were used to perform the simulations. For solving the minimization program in (10), we used the nonlinear conjugate gradient with the backtracking line search method [7]. The measurement matrix () was the discrete Fourier transform matrix while the daubechies4 wavelet transform matrix () and TV were used as sparsifying transforms.
Fourier encoded MCSMRI The reference brain image was multiplied by the sensitivity function to generate eight sensitivity encoded images. The Fourier transform of each these images was taken in the PE direction; only a few PEs were taken according to the Gaussian probability distribution function. MCSMRI reconstruction of (10) was solved using the nonlinear conjugate gradient on this data. An example sampling scheme for the Fourier encoded MCSMRI is shown in Fig. 3(a). Noiselet encoded MCSMRI A Noiselet transform of the sensitivity encoded images was taken in the PE direction, with only a few PE selected according to the uniform probability distribution function. MCSMRI reconstruction of (10) was solved using the nonlinear conjugate gradient on this data. An example of the sampling scheme for noiselet encoded MCSMRI is shown in Fig. 3(b).
For a noiseless simulation, the reconstructed images for different acceleration factors (4, 8 and 16) are shown in Fig. 8. The difference images in Fig. 8 (d)(f) and (j)(l) demonstrate that the error in noiselet encoding is always less than in Fourier encoding, and that the noiselet encoded MCSMRI reconstruction preserves spatial resolution better than the Fourier encoded MCSMRI. Fig. 8 (m) and (n) show the zoomed images reconstructed with Fourier encoding for acceleration factors of 8 and 16 respectively, while Fig. 8 (o) and (p) show the zoomed images reconstructed with noiselet encoding for an acceleration factors of 8 and 16 respectively. The zoomed images highlight that the spatial resolution of the noiselet encoded reconstructions outperforms the Fourier encoded reconstructions. Moreover, the spatial resolution provided by the noiselet encoding at an acceleration factor of 16 is comparable to that of the Fourier encoding at an acceleration factor of 8, suggesting that noiselet encoding performs approximately twice as good as Fourier encoding.
To measure the relative error, simulations were performed on the brain image for 1000 times by randomly generating a sampling mask each time. The mean of the relative errors was calculated after 1000 such reconstructions at every acceleration factor. The mean relative error versus the number of measurements is plotted in Fig. 9 and highlights that noiselet encoding outperforms Fourier encoding for all acceleration factors. The relative error for noiselet encoding at an acceleration factor of 16 was the same as the relative error for Fourier encoding at an acceleration factor of 8 indicating that higher acceleration factors are achievable with noiselet encoding compared to Fourier encoding.
In practice, MR data always has some noise and the level of noise depends upon many factors including the FOV, resolution, type of imaging sequence, magnetic field inhomogeneity and RF inhomogeneity. Therefore, simulations were carried out to evaluate the performance of both the noiselet encoding and Fourier encoding schemes in the presence of variable levels of noise. Different levels of random Gaussian noise were added to the measured kspace data, and MCSMRI reconstruction was performed for noiselet and Fourier encoding schemes. For every level of noise, 1000 simulations were performed and the mean of the relative error was calculated. Fig. 10 shows the mean relative error as a function of the Signal to Noise Ratio (SNR), demonstrating the comparative performance of noiselet encoding reconstructions and Fourier encoding reconstructions in the presence of noise. The plots demonstrate that noiselet encoding outperforms Fourier encoding for SNR above 20 dB for all acceleration factors, but does a poor job at SNR of 10 dB. However, as shown in the images in Fig. 10(d), only the images with SNR above 20 dB are adequate for diagnostic purposes. Thus, acquisitions with 10 dB SNR is not a viable scanner operation condition and the poor performance of noiselet encoding at 10 dB SNR is not a practical limitation. The poor performance of noiselet encoding at 10 dB SNR can be attributed to the fact that at extremely low SNR, most of the noiselet coefficients are severely corrupted by the noise since their magnitudes are approximately uniform. In contrast, the Fourier coefficients at the center of kspace have much larger magnitudes and hence are less affected by the noise at low SNR. These large magnitude coefficients are fully utilized in reconstruction because of the centralized variable density sampling scheme, hence Fourier encoding is less affected by the noise and performs better at low SNR.
5 Experiments
Experiments were carried out on Siemens Skyra 3T MRI scanner with a maximum gradient strength of 40 mT/m and a maximum slew rate of 200 mT/m/sec. Informed consent was taken from healthy volunteers in accordance with the Institution’s ethics policy. To validate the practical implementation of noiselet encoding, the pulse sequence shown in Fig. 1 was used to acquire noiselet encoded data. An RF excitation pulse with 256 points was used with the duration equal to 5.12 ms and flip angle of 10. We also acquired the Fourier encoded data using the spin echo (SE) sequence to compare the quality of the reconstructed image from the data acquired by the noiselet encoding sequence. An apodized slice selective sinc RF excitation pulse was used in the spin echo sequence with a duration of 2.56 ms and a flip angle of 10. The protocol parameters for the noiselet encoding sequence and the Fourier encoding SE sequence were (i) Phantom experiments FOV = 200 mm, TE/TR = 26/750 ms, averages = 2, image matrix = 256256; and (ii) In vivo experiments FOV = 240 mm, TE/TR = 26/750 ms, averages = 2, image matrix = 256256.
The performance of the MCSMRI reconstruction depends on the accuracy of the sensitivity matrix estimated. We used the regularized selfcalibrated estimation method [46] to estimate the sensitivity maps from the acquired data. This method estimates the sensitivity map of the th receive coil by using
(15) 
where and is a spatial roughness penalty function with weighting factor . The reference image can be obtained using the sum of squares of individual coil images ’s. For the experimental results presented below, the sensitivity maps were estimated from fully sampled images using (15). The data was acquired using a 32 channel head coil, but only the data from 14 channels with good SNR was used in the reconstruction.
NonFourier encoding in general is sensitive to field inhomogeneities, but careful design of the sequence and good shimming can result in high quality images. To reconstruct the noiselet encoded data the inverse Fourier transform was taken along the frequency encoding axis and the inverse noiselet transform was taken along the PE axis. To reconstruct the Fourier encoded data, an inverse Fourier transform was taken along both axes. Fig. 11 shows the images reconstructed from the noiselet encoded data and Fourier encoded data sets. These images demonstrate that the noiselet encoding reconstructions are practically feasible and produce artifact free images. Fig. 11(c) shows a zoomed portion of the noiselet encoded image, while Fig. 11(f) shows a zoomed portion of the Fourier encoded image. The zoomed images reveal that the resolution of the image from noiselet encoding with 256 noiselet excitation is the same as that of the image from Fourier encoding with 256 phase encodes. Fig. 11 (g) and (i) show the T2 weighted images for the brain with noiselet encoding and Fourier encoding, respectively. Fig. 11 (h) and (j) show the T1 weighted images for the brain with noiselet encoding and Fourier encoding, respectively. It is evident from the in vivo images that the proposed noiselet encoding is feasible in practice.
To validate the feasibility of the proposed reconstruction method, we performed retrospective undersampling on the acquired noiselet encoded data and Fourier encoded data to simulate accelerated data acquisition. After retrospective undersampling, the unconstrained optimization program (10) was solved using the nonlinear conjugate gradient method to reconstruct the desired image for different acceleration factors. Fig. 12 (a)(c) shows the reconstructed images for the acceleration factors of 4, 8 and 16 on the Fourier encoded data while Fig. 12 (d)(f) shows the corresponding difference images. Similarly, Fig. 12 (g)(i) shows the reconstructed images for the acceleration factors of 4, 8 and 16 on the noiselet encoded data, and Fig. 12 (j)(l) shows the corresponding difference images for noiselet encoded MCSMRI. These results on the acquired data are consistent with the simulation results and indicate that the noiselet encoding is superior to the Fourier encoding in preserving resolution.
Fig. 12 (AH) shows the zoomed portion of the reconstructed images with Fourier encoding and noiselet encoding. One can distinguish between the small dots in the zoomed images reconstructed with noiselet encoding while it is difficult to distinguish these dots in the images reconstructed with Fourier encoding. This demonstrates that noiselet encoding is able to preserve resolution better than the Fourier encoding. Fig. 13 show the images reconstructed with Fourier encoding and noiselet encoding for various acceleration factors on the data acquired for one axial slice of the brain. Since the SNR of the in vivo images is less than in the phantom images, reconstruction is shown only up to an acceleration factor of 8. The difference images demonstrate that noiselet encoding outperforms Fourier encoding for all acceleration factors. In particular, at the acceleration factor of 8 the image reconstructed with Fourier encoded data has significantly poorer resolution compared to the image reconstructed with noiselet encoded data.
Fig. 14 shows the images reconstructed with noiselet and Fourier encodings using 3D GRE sequence. In our implementation of noiselet encoding in 2D spin echo sequence, the flip angle of 10 was used, which results in loss of some available SNR. Therefore, we have implemented noiselet encoding in 3D Gradient Echo (GRE) sequence as shown in Fig. 1 (b), which uses noiselets encoding in one direction and Fourier encoding in other two directions. The noiselet encoding is performed using specially selective RF excitation pulse, while the Fourier encoding is performed using gradients. Phantom data was acquired using this sequence with the following parameters: FOV = 200200 mm, TE/TR = 6.5/13 ms, slices = 32, FOV of slice = 160 mm, flip angle = 5 and (noiselets) phase encodes = 256. A 3D Fourier encoded data was also acquired with exactly same parameters. It can be seen from the images that the noiselet encoding provides similar quality of images to that of Fourier encoding. These images are only shown here to demonstrate the feasibility that noiselet encoding can be implemented in 3D GRE sequence.
6 Discussion
We have presented a new method of encoding the MR data in PE direction with noiselet basis functions for accelerating MRI scans using MCSMRI reconstruction. The simulation results demonstrate that the proposed encoding gives rise to a multichannel measurement matrix with improved RIP, and the reconstruction method using the noiselet bases outperforms the conventional Fourier encoding scheme. The mean relative error for noiselet encoding at the acceleration factor of 16 is comparable to that of Fourier encoding at the acceleration factor of 8, demonstrating that higher acceleration factors can be achieved with noiselet encoding than the Fourier encoding in the MCS framework.
The reconstruction from the noiselet encoding scheme preserves image spatial resolution far better than the Fourier encoding scheme. The Fourier encoding scheme intelligently exploits the property of kspace since most of the energy is concentrated at the center of the kspace. Therefore, densely sampling the center and randomly undersampling the outer regions of the kspace retains most of the energy in the acquired data. However retaining most of the energy does not imply that most of the information is captured in the acquisition. The low energy (high frequency) component in the outer kspace contains the information about the fine features of the image that the variable density random undersampling pattern fails to capture. Therefore, the images reconstructed with the Fourier encoding scheme look visually good but have reduced resolution due to insufficient information about the high frequency components in the acquired data. On the other hand, noiselet encoding completely spreads out the energy of the signal in the measurement domain. Therefore each measurement in the noiselet domain has sufficient information to reconstruct the fine details of the image, thus preserving the resolution better than the Fourier encoding.
Noiselet basis functions have some interesting properties that can be exploited, for example noiselets are unitary basis functions and have complex conjugate symmetry. This conjugate symmetry property can be exploited in a way similar to that of the Fourier encoding for partial acquisition [47, 48]. Another interesting property of noiselet basis functions, as for the Fourier basis functions, is that regular undersampling in the noiselet domain results in aliasing in the image domain. Therefore in the case of regular undersampling in the noiselet domain, unaliasing with SENSE [8] alone can also be used for noiselet encoding.
In general the implementation of nonFourier encoding suffers from a few limitations, and hence the current implementation of noiselet encoding also suffers from these limitations as summarized below. (i) In 2D imaging implementation of noiselet encoding, the excitation of noiselet profile is not slice selective, thus a slice selective 180 pulse is required after excitation, which limits the noiselet encoding to spin echo type sequences. The spin echo sequence always have long TR, therefore the proposed noiselet encoding can only be used for applications requiring long TR, such as structural scans, but will be of little use in dynamic imaging. (ii) In the current implementation, to simplify the design of noiselet excitation pulse we have used direct Fourier transform method that limits the excitation to the low flip angle regime, resulting in the sacrifice of some available SNR. (iii) The duration of noiselet RF excitation pulse is long compared to the conventional sinc RF pulse and the noiselet encoding can only be implemented in one direction in the current implementation. (iv) Due to dielectric effects, etc, in practice, B1 field is always not perfectly homogeneous. Because B1 is used for spatial encoding in the proposed noiselet encoding scheme, B1 inhomogeneity may introduce some perturbations to the noiselet measurement matrix, which in turn may result in some image artifact if the perturbation is large.
The above limitations are not specific to noiselet encoding scheme but are common to all nonFourier encoding schemes. Here we discuss some probable solutions to the above mentioned problems. (i) The problem of slice selection can be alleviated using a 3D gradient echo (GRE) sequences where a 3D volume can be excited with a noiselet profile in one dimension and other two dimensions can be Fourier encoded. A demonstrative pulse sequence for this solution has been given in Section III and its scan result has been given in Section V, showing the feasibility of this solution. (ii) The low tip angle in current implementation is the limitation of the direct Fourier transform method used to compute the RF pulse. It is not an intrinsic limitation of noiselet encoding. Although difficult, computation of large tip angle noiselet RF pulse is possible by using nonlinear computation methods such as direct iterative solution of Bloch equation [49] and the SLR method[50, 51], which are currently being investigated. (iii) The duration of the RF pulse can be reduced by using paralleltransmit and multiple dimensional excitation of noiselet profiles to achieve encoding [52, 49]. This is our ongoing research. (iv) The perturbations to the measurement matrix induce an equivalent deterministic noise additive to the measured MR signals. When the inverse noiselet transform is applied directly to the fully sampled dataset to reconstruct the image, a structured artifact may show up if the perturbation is large. This problem can be alleviated when the CS reconstruction method as given in (10) is used for image reconstruction. This is because the CS reconstruction algorithm has inherent denoising capability, which can suppress small perturbations by enforcing the prescribed bound on the reconstruction error. See Section II A and the references therein for detailed discussions. For this reason and also because of the high quality of the new 3T scanner used in our experiments, we have not observed structured image artefacts in the experiments presented in Section V.
7 Conclusion
In this paper we have introduced a method of acquiring data in the noiselet domain and presented a method for the design and implementation of pulse sequences to acquire data in the noiselet domain. The performance of the noiselet encoding has been thoroughly evaluated by extensive numerical analysis, simulation and experiments. The results indicate that the multichannel noiselet measurement matrix has better RIP than that of its Fourier counterpart, and that the noiselet encoding scheme in MCSMRI outperforms the conventional Fourier encoding in preserving image resolution, and can achieve higher acceleration factors than the conventional Fourier encoding scheme. The implementation of noiselet encoding by tailoring spin echo and gradient echo sequences demonstrates that the proposed encoding scheme is pragmatic. The proposed technique has the potential to accelerate image acquisition in applications that require high resolution images.
References
 [1] D. L. Donoho, “Compressed sensing,” Information Theory, IEEE Transactions on, vol. 52, no. 4, pp. 1289–1306, 2006.
 [2] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” Information Theory, IEEE Transactions on, vol. 52, no. 2, pp. 489–509, 2006.
 [3] E. Candes and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse problems, vol. 23, no. 3, p. 969, 2007.
 [4] E. J. Candes and J. K. Romberg, “Signal recovery from random projections,” in Electronic Imaging 2005, pp. 76–86, International Society for Optics and Photonics, 2005.
 [5] E. J. Candes and T. Tao, “Nearoptimal signal recovery from random projections: Universal encoding strategies,” Information Theory, IEEE Transactions on, vol. 52, no. 12, pp. 5406–5425, 2006.
 [6] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing MRI,” Signal Processing Magazine, IEEE, vol. 25, no. 2, pp. 72–82, 2008.
 [7] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magnetic resonance in medicine, vol. 58, no. 6, pp. 1182–1195, 2007.
 [8] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, P. Boesiger, et al., “SENSE: sensitivity encoding for fast MRI,” Magnetic resonance in medicine, vol. 42, no. 5, pp. 952–962, 1999.
 [9] D. K. Sodickson and W. J. Manning, “Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays,” Magnetic Resonance in Medicine, vol. 38, no. 4, pp. 591–603, 1997.
 [10] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase, “Generalized autocalibrating partially parallel acquisitions (GRAPPA),” Magnetic Resonance in Medicine, vol. 47, no. 6, pp. 1202–1210, 2002.
 [11] R. Otazo and D. Sodickson, “Distributed compressed sensing for accelerated MRI,” in Proceedings of the 17th Annual Meeting of ISMRM, p. 378, 2009.
 [12] D. Liang, B. Liu, J. Wang, and L. Ying, “Accelerating SENSE using compressed sensing,” Magnetic Resonance in Medicine, vol. 62, no. 6, pp. 1574–1584, 2009.
 [13] D. Liang, B. Liu, and L. Ying, “Accelerating sensitivity encoding using compressed sensing,” in Engineering in Medicine and Biology Society, 2008. EMBS 2008. 30th Annual International Conference of the IEEE, pp. 1667–1670, IEEE, 2008.
 [14] R. Otazo, D. Kim, L. Axel, and D. K. Sodickson, “Combination of compressed sensing and parallel imaging for highly accelerated firstpass cardiac perfusion MRI,” Magnetic Resonance in Medicine, vol. 64, no. 3, pp. 767–776, 2010.
 [15] M. Lustig and J. M. Pauly, “SPIRiT: Iterative selfconsistent parallel imaging reconstruction from arbitrary kspace,” Magnetic Resonance in Medicine, vol. 64, no. 2, pp. 457–471, 2010.
 [16] J. X. Ji, C. Zhao, and T. Lang, “Compressed sensing parallel magnetic resonance imaging,” in Engineering in Medicine and Biology Society, 2008. EMBS 2008. 30th Annual International Conference of the IEEE, pp. 1671–1674, IEEE, 2008.
 [17] M. Lustig, J. M. Santos, D. L. Donoho, and J. M. Pauly, “kt SPARSE: High frame rate dynamic MRI exploiting spatiotemporal sparsity,” in Proceedings of the 13th Annual Meeting of ISMRM, Seattle, p. 2420, 2006.
 [18] H. Jung, K. Sung, K. S. Nayak, E. Y. Kim, and J. C. Ye, “kt FOCUSS: A general compressed sensing framework for high resolution dynamic MRI,” Magnetic Resonance in Medicine, vol. 61, no. 1, pp. 103–116, 2009.

[19]
C. Qiu, W. Lu, and N. Vaswani, “Realtime dynamic MR image reconstruction using Kalman filtered compressed sensing,” in
Acoustics, Speech and Signal Processing, 2009. ICASSP 2009. IEEE International Conference on, pp. 393–396, IEEE, 2009.  [20] H. Jung and J. C. Ye, “Motion estimated and compensated compressed sensing dynamic magnetic resonance imaging: What we can learn from video compression techniques,” International Journal of Imaging Systems and Technology, vol. 20, no. 2, pp. 81–98, 2010.
 [21] E. J. Candès, “The restricted isometry property and its implications for compressed sensing,” Comptes Rendus Mathematique, vol. 346, no. 9, pp. 589–592, 2008.
 [22] A. M. Tillmann and M. E. Pfetsch, “The computational complexity of the restricted isometry property, the nullspace property, and related concepts in compressed sensing,” arXiv preprint arXiv:1205.2081, 2012.
 [23] E. J. Candes and T. Tao, “Decoding by linear programming,” Information Theory, IEEE Transactions on, vol. 51, no. 12, pp. 4203–4215, 2005.

[24]
M. Bayati, M. Lelarge, and A. Montanari, “Universality in polytope phase transitions and iterative algorithms,” in
Information Theory Proceedings (ISIT), 2012 IEEE International Symposium on, pp. 1643–1647, IEEE, 2012.  [25] J. P. Haldar, D. Hernando, and Z.P. Liang, “CompressedSensing MRI with random encoding,” Medical Imaging, IEEE Transactions on, vol. 30, no. 4, pp. 893–903, 2011.
 [26] J. D. Blanchard, “Towards deterministic compressed sensing,” in Proceedings of the National Academy of Sciences,, pp. 1147–1148, 2013.
 [27] G. Puy, J. P. Marques, R. Gruetter, J. Thiran, D. Van De Ville, P. Vandergheynst, and Y. Wiaux, “Spread spectrum magnetic resonance imaging,” Medical Imaging, IEEE Transactions on, vol. 31, no. 3, pp. 586–598, 2012.
 [28] D. Liang, G. Xu, H. Wang, K. F. King, D. Xu, and L. Ying, “Toeplitz random encoding mr imaging using compressed sensing,” in Biomedical Imaging: From Nano to Macro, 2009. ISBI’09. IEEE International Symposium on, pp. 270–273, IEEE, 2009.
 [29] F. Sebert, Y. Zou, and L. Ying, “Compressed Sensing MRI with random B1 field,” Proc ISMRM,(Toronto, Canada, 2008), p. 3151, 2008.
 [30] H. Wang, D. Liang, K. King, and L. Ying, “Toeplitz random encoding for reduced acquisition using compressed sensing,” Proc ISMRM,(Honolulu, Hawai, 2009), p. 2669, 2009.
 [31] X. Qu, Y. Chen, X. Zhuang, Z. Yan, D. Guo, and Z. Chen, “Spread spectrum compressed sensing MRI using chirp radio frequency pulses,” arXiv preprint arXiv:1301.5451, 2013.
 [32] Z. Liu, B. Nutter, and S. Mitra, “Compressive sampling in fast waveletencoded MRI,” in Image Analysis and Interpretation (SSIAI), 2012 IEEE Southwest Symposium on, pp. 137–140, IEEE, 2012.
 [33] H. Wang, D. Liang, K. F. King, and L. Ying, “Threedimensional hybridencoded MRI using compressed sensing,” in Biomedical Imaging (ISBI), 2012 9th IEEE International Symposium on, pp. 398–401, IEEE, 2012.
 [34] H. Wang, D. Liang, K. F. King, G. Nagarsekar, and L. Ying, “Threedimensional hybridencoded MRI using compressed sensing,” in International Society of Magnetic Resonance in Medicine Scientific Meeting, p. 73, 2012.
 [35] R. Coifman, F. Geshwind, and Y. Meyer, “Noiselets,” Applied and Computational Harmonic Analysis, vol. 10, no. 1, pp. 27–44, 2001.
 [36] T. Tuma, P. Hurley, et al., “On the incoherence of noiselet and Haar bases,” in SAMPTA’09, International Conference on Sampling Theory and Applications, 2009.
 [37] S. Datta, K. Ni, P. Mahanti, and S. Roudenko, “Stability of efficient deterministic compressed sensing for images with chirps and reedmuller sequences,” Applied Mathematics, vol. 4, p. 183, 2013.
 [38] K. Ni, S. Datta, P. Mahanti, S. Roudenko, and D. Cochran, “Efficient deterministic compressed sensing for images with chirps and reedmuller codes,” SIAM Journal on Imaging Sciences, vol. 4, no. 3, pp. 931–953, 2011.
 [39] S. S. R. C. L. Applebaum, S. D. Howard, “Chirp sensing codes: Deterministic compressed sensing measurements for fast recovery,” Applied and Computational Harmonic Analysis, vol. 26, pp. 283–290, 2009.
 [40] M. Elad, Sparse and redundant representations: from theory to applications in signal and image processing. Springer, 2010.
 [41] S. Ma, W. Yin, Y. Zhang, and A. Chakraborty, “An efficient algorithm for compressed MR imaging using total variation and wavelets,” in Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE Conference on, pp. 1–8, IEEE, 2008.
 [42] L. P. Panych, G. P. Zientara, and F. A. Jolesz, “MR image encoding by spatially selective rf excitation: An analysis using linear response models,” International journal of imaging systems and technology, vol. 10, no. 2, pp. 143–150, 1999.
 [43] D. Mitsouras, W. S. Hoge, F. J. Rybicki, W. E. Kyriakos, A. Edelman, and G. P. Zientara, “NonFourierencoded parallel MRI using multiple receiver coils,” Magnetic resonance in medicine, vol. 52, no. 2, pp. 321–328, 2004.
 [44] L. P. Panych, G. P. Zientara, P. Saiviroonporn, S.S. Yoo, and F. A. Jolesz, “Digital waveletencoded MRI: A new waveletencoding methodology,” Journal of Magnetic Resonance Imaging, vol. 8, no. 5, pp. 1135–1144, 1998.
 [45] L. P. Panych, C. Oesterle, G. P. Zientara, and J. Hennig, “Implementation of a fast gradientecho SVD encoding technique for dynamic imaging,” Magnetic resonance in medicine, vol. 35, no. 4, pp. 554–562, 1996.
 [46] Y. Kim, J. Fessler, and D. Noll, “Smoothing effect of sensitivity map on fMRI data using a novel regularized selfcalibrated estimation method,” in Proc. Intl. Soc. Mag. Res. Med, p. 1267, 2008.
 [47] P. M. Margosian, G. DeMeester, and H. Liu, “Partial Fourier acquisition in MRI,” eMagRes, 1996.
 [48] G. McGibney, M. Smith, S. Nichols, and A. Crawley, “Quantitative evaluation of several partial Fourier reconstruction algorithms used in MRI,” Magnetic resonance in medicine, vol. 30, no. 1, pp. 51–59, 1993.
 [49] D. Xu, K. F. King, Y. Zhu, G. C. McKinnon, and Z.P. Liang, “Designing multichannel, multidimensional, arbitrary flip angle RF pulses using an optimal control approach,” Magnetic Resonance in Medicine, vol. 59, no. 3, pp. 547–560, 2008.
 [50] J. Pauly, P. Le Roux, D. Nishimura, and A. Macovski, “Parameter relations for the ShinnarLe Roux selective excitation pulse design algorithm,” Medical Imaging, IEEE Transactions on, vol. 10, no. 1, pp. 53–65, 1991.
 [51] J. K. Barral, J. M. Pauly, and D. G. Nishimura, “SLR RF pulse design for arbitrarilyshaped excitation profiles,”
 [52] W. A. Grissom, D. Xu, A. B. Kerr, J. A. Fessler, and D. C. Noll, “Fast largetipangle multidimensional and parallel RF pulse design in mri,” Medical Imaging, IEEE Transactions on, vol. 28, no. 10, pp. 1548–1559, 2009.