WO1992003089A1 - Analyzing heart wall motion using spatial modulation of magnetization - Google Patents

Abstract

A technique for generating 'optimal' radio frequency pulses for use in creating SPAMM imaging stripes and for quantitatively determining the displacements of the imaging stripes during a pre-imaging interval. The imaging operator selects the desired stripe parameters (304) so that the resulting SPAMM imaging stripe will have the desired narrowness, sharpness, flatness and the like. These parameters are then used by an optimal pulse generating system (212) to generate the input radio frequency pulse sequence which will produce the desired stripes. Tagging of the heart wall or some other soft tissue structure with a grid of planes of altered magnetization is then used to provide finite element analysis to quantify regional heart wall motion and the like. The intersections of the tagging stripes are used as a set of fiducial marks within the heart wall, and the motion of these intersections through the heart cycle is used to track the motion of the underlying tissue.

Description

ANALYZING HEART WALL MOTION USING SPATIAL MODULATION

OF MAGNETIZATION

COPYRIGHT NOTICE AND AUTHORIZATION

A portion of the disclosure of this patent document

(APPENDICES A and B) contains material which is subject to copyright protection. The copyright owner has no objection to facsimile reproduction by anyone of the patent document or the patent disclosure, as it appears in the Patent and Trademark Office patent file or records, but otherwise reserves all copyright rights whatsoever.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation-in-part application of co-pending U.S. Patent Application Serial Number 07/570,207, filed August 17, 1990 and entitled "System and Method For Analyzing Heart Wall Motion Using Tagged Magnetic Resonance Images", which is a continuation-in-part application of co-pending U.S. Patent Application Serial Number 07/254,454, filed October 6, 1988 and entitled "Magnetic Resonance Imaging Using Spatial Modulation of Magnetization". This application is also related to co-pending patent application Serial Number 07/655,077, filed February 14, 1991 and entitled "Methods of Generating Pulses For Selectively Exciting Frequencies", which is a continuation-in-part application of U.S. Patent Application Serial Number 07/471,707, filed January 29, 1990, now abandoned, which is, in turn, a continuation application of U.S. Patent Application Serial Number 07/176,429, filed April 1, 1988, now abandoned.

REFERENCE TO MICROFICHE APPENDIX APPENDIX B is a computer program listing which is available in microfiche upon request.

BACKGROUND OF THE INVENTION

Field of the invention

The present invention relates to a method for generating optimal pulses for magnetic resonance imaging using spatial modulation of magnetization (SPAMM), and more particularly, to a method for specifying the desired magnetization in certain frequency ranges and the desired smoothness for each range and then using the specified values to generate a pulse sequence which will optimally achieve the desired magnetization and smoothness for imaging. A system and method is also disclosed which adapts the techniques of motion analysis from radiologic imaging of invasively implanted metallic myocardial markers to the techniques of forming tagged magnetic resonance images of heart wall motion using optimal pulse sequences.

Description of the Prior Art

Nuclear magnetic resonance (NMR) has many applications in chemistry and biology. Recent advances have made it possible to use NMR to image the human body. Because much of the body is water, which gives an excellent NMR signal, NMR can yield detailed images of the body and information about various disease processes.

The physical basis of NMR is that many nuclei, such as protons, have a magnetic spin. When an external field, B₀ (in the z axis direction), is applied to a system whose elements have a magnetic spin, the individual elements of the system precess around B₀ according to the Larmor frequency. The Larmor frequency, for typical values of the external field B₀ and for typical nuclei being investigated, is in the radio frequency (rf) range.

Generally, in NMR the system is perturbed by an alternating signal having a frequency ω at or near the Larmor frequency of the spins of interest, whose magnetic field (B₁) is normal to B₀. This field is called a pulse. Any spin with a Larmor frequency close to ω will interact with the applied field and rotate away from the z axis towards the xy plane. This interaction is called excitation of the spins. The precise axis of rotation and amount of rotation depend on the strength of the pulse, the duration of the pulse, and the difference between ω and ω₀. If B₁ is very strong in relation to the off-resonance effects (ω-ω₀) , then the rotation will be essentially around an axis in the xy plane and is called a hard pulse. Generally, the stronger the B_x field, the broader range of frequencies that it will perturb. The longer the B₁ field is on, and the stronger the B₁ field is, the more the spins will be perturbed.

In a typical application, the spins of interest will have different ω₀ values, either because they have a different environment, or because B₀ varies with location (a gradient). Furthermore, one needs to have different types of excitation of the spins depending on the application. The three most commonly needed types of excitation are described as a π/2 pulse, which rotates spins from the z axis to the xy plane, a π pulse, which rotates spins from pointing in the positive z direction towards the negative z direction, and a refocusing pulse, which rotates spins of interest 180° about an axis in the xy plane, the axis being the same for all spins of interest. However, other types of excitation profiles can be needed, and one may need to apply a different type of excitation to different frequencies.

NMR imaging schemes depend on applying magnetic field gradients so that different nuclei at different locations experience different magnetic fields and, therefore, have different frequencies. Accordingly, the location of the nucleus determines its frequency. One then applies a pulse, which rotates, or excites, the nuclei, and therefore, the magnetic dipole moment. One tries to excite only nuclei which have frequencies corresponding to the slice of tissue which it is desired to image, and to excite those nuclei to the same degree.

Thus, it is necessary to design pulses which perturb a particular range of frequencies to the same extent but that do not perturb all other frequencies. This "ideal" pulse is physically unrealizable. However, it has been possible to synthesize pulses which will do this more or less perfectly.

As a result, much work in the NMR imaging art has gone into synthesizing pulses which will give very "sharp" slices, for the sharper the slice, the better the resolution of the image.

Previously, so-called "hard" radio frequency pulses have been used to give sharp slices for imaging. As known by those skilled in the art, a "hard" pulse is a constant amplitude rectangular radio frequency pulse in which the strength of the applied alternating field is sufficiently large that the pulse can be assumed to affect all frequencies of interest equally. Selective excitation is thus achieved by applying several pulses in a row and varying pulse widths and inter-pulse delay times. Different frequencies precess for different amounts during the delays between the pulses and thus respond differently to the pulse sequence. Many pulse sequences using hard pulses have been devised for different applications in NMR imaging, as described by P.J. Hore in an article entitled "Solvent Suppression in Fourier Transform Nuclear Magnetic Resonance", Journal of Magnetic Resonance. Vol. 55, pp. 283-300 (1983), for example, and by Brandes in U.S. Patent No. 4,695,798.

However, hard pulse sequences have several major limitations. First the frequency response of hard pulses have "side lobes" or harmonics around the desired frequency range. The location of the side lobes depends on the delay between pulses. Second, if the frequency range of interest is large, as it is in NMR imaging, it may be difficult to create pulses that are strong enough to be considered "hard". Furthermore, the frequency response of hard pulse sequences currently in use is far from ideal.

Despite these limitations, hard pulses have been used effectively in the technique called spatial modulation of magnetization (SPAMM), which utilizes the motion sensitivity of magnetic resonance signals to generate images representative of heart wall motion. One of the principal sources of such motion sensitivity is that when the local magnetization of a material is altered, it carries the altered magnetization with it when it moves (within the limits of the relaxation times). It has been known to measure blood flow by locally altering the magnetization of blood and then detecting the passage of this "tagged" blood downstream. Analogous techniques have been proposed for measurement of blood flow with magnetic resonance imaging (MRI) by detecting the passage of blood with locally saturated or inverted magnetization through the region being imaged. Other techniques have been proposed by Zerhouni et al. in "Human Heart: Tagging With MRI Imaging - A Method For Non-Invasive Assessment of Myocardial Motion", Radiology, Vol. 169, 1988, pp. 59-63, for imaging of myocardial motion within the plane of the image by using selective saturation pulses to magnetically tag selected bands within the heart wall, thereby producing stripes of altered magnetization that move with the tagged tissue. However, such selective saturation for producing tagging stripes is limited by the number of tags that can be practically produced. Thus, although regional intramural heart wall motion has been followed by using the ability of MRI to produce and track regions of altered magnetization, it has heretofore been impractical for clinical use.

Axel et al. describe the use of SPAMM for following heart wall motion in a clinical setting in articles entitled "MR Imaging of Motion With Spatial Modulation of Magnetization", Radiology. Vol. 171, June 1989, pp. 841-845 and "Heart Wall Motion: Improved Method of Spatial Modulation of Magnetization for MR Imaging", Radiology. Vol. 172, August 1989, pp. 349-350. Heart wall motion may also be followed by producing signal phase shifts proportional to velocity as described by P. van Dijk in "Direct Cardiac NMR Imaging of Heart Wall and Blood Flow Velocity", Journal of Computer Assisted Tomography. Vol. 8, June 1984, pp. 429-436. Nayler et al. in "Blood Flow Imaging by Cine Magnetic Resonance", Journal of Computer Assisted Tomography, Vol. 10, September/October 1986, pp. 715-722 and Wedeen in U.S. Patent No.4,752,734 describe similar techniques for monitoring blood flow. Although, the above-referenced publications have described methods for production of images that can display effects of regional heart wall motion, they have not considered the quantitative analysis of the resulting images.

Prior art imaging approaches for analyzing regional heart wall motion have used contrast angiography or radioactive labeling of red cells to study the changing shape of the endocardial surface or tomographic imaging with ultrasound, x-rays or magnetic resonance to study the changes in shape of the inner and outer surfaces of the heart wall or its local change in thickness. However, the lack of any significant landmarks within the heart wall prevents using these techniques to study intramural distributions of motion or any components of motion other than radial components. Moreover, the motion of the curved heart wall through the fixed imaging plane often creates an artificial apparent motion within the imaging plane, thereby resulting in an artificial appearance of active thickening.

To overcome the limitations imposed by the restriction to imaging of the surfaces, rather than internal tissue, radiopaque markers have previously been implanted in the heart in order to permit the tracking of individual points within the myocardium during contraction of the heart. In particular, Meier et al. disclose in an article entitled "Kinematics of the Beating Heart", IEEE Transactions on Biomedical Engineering,. Vol. 27, No. 6, June 1980, pp. 319- 329 and in an article entitled "Contractile Function in Canine Right Ventricle", American Journal of Physiology. Vol. 239, 1980, pp. H794-H804, a technique whereby at least three small radiopaque markers are implanted in a living heart such that their movement can be ascertained. Meier et al. describe a theoretical framework of kinematics which accounts for the large time-varying displacements of the markers whereby, for arbitrary small segments of myocardium, displacements can be described completely by a rotation tensor and a stretch tensor. Similar results are disclosed for three-dimensional motion analysis of the myocardium by Fenton et al. in an article entitled "Transmural Myocardial Deformation in the Canine Left Ventricular Wall", American Journal of Physiology. Vol. 235, 1978, pp. H523-H530 and by Waldman et al. in an article entitled "Transmural Myocardial Deformation in the Canine Left Ventricle", Circulation Research. Vol. 57, 1985, pp. 152-163. However, the invasive nature of the implantation of the radiopaque markers in accordance with these techniques coupled with the difficulty of tracking a sufficient number of such markers to adequately define regional wall motion has largely restricted the techniques disclosed in the above- referenced articles to research applications. In other words, because of the invasiveness of these procedures, these procedures are presently unacceptable for clinical use.

SPAMM has been used to simultaneously create parallel sheets of altered magnetization that show up as stripes in subsequent images. Motion of the "tagged" tissue between the times of SPAMM application and imaging results in a corresponding displacement of the stripes, thereby permitting study of motions. It is desired to adapt the MR imaging techniques for monitoring regional motion with the heart wall as described above for use with the analysis techniques developed by Meier and others for imaging of invasively implanted radiopaque markers. In particular, it is desired to adapt the analysis techniques developed for imaging of invasively implanted radiopaque markers to MR studies of heart wall motion using the noninvasive MR imaging techniques of SPAMM whereby actual displacements within the myocardium of the heart may be measured. It is also desired to find pulse sequences with the desired frequency response and amplitude such that the best possible SPAMM imaging stripes can be used for such analysis. For example, in an article entitled "Heart Wall Motion: Improved Method of Spatial Modulation of Magnetization for MR Imaging", Radiology, Vol. 172, August 1989, pp. 349-350, one of the present inventors implements SPAMM using a binomial series of radio frequency pulses which are separated by equal gradient pulses. In this sequence, during the radio frequency pulses all spins are on resonance, so the pulses are "hard", i.e., their effect on the spins is independent of the spatial location. The frequency encoding then occurs during the gradient pulses. The result is a sequence of hard pulses with free precession between the pulses. Because of the nice stripes generated when such binomial pulses were used, binomial pulses have been preferred for use with SPAMM. However, binomial pulse sequences provide limited flexibility because the relationships of the amplitudes of the binomial pulses are predefined and cannot be adjusted by the SPAMM user.

Others have also expended considerable effort in devising pulses which have better frequency characteristics. One of the best of these is the hyperbolic secant pulse described by M.S. Silver et al. in an article entitled "Highly Selective π/2 and π Pulse Generation", Journal of Magnetic Resonance. Vol. 59, pp. 347-351 (1984). Unfortunately, even the response to this pulse does not have an ideal frequency response. Moreover, implementation of this pulse requires a spectrometer which can simultaneously modulate the amplitude and phase of the applied external field, which many spectrometers cannot do. Indeed, patents have issued on the instrumentation for the delivery of certain shaped pulses to the transmitter, such as the afore-mentioned patent to Brandes. Further continuing efforts to devise improved pulses are described by, for example, H. Yan and J. Gore in an article entitled "Improved Selective 180° Radio Frequency Pulses for Magnetization Inversion and Phase Reversal", Journal of Magnetic Resonance. Vol. 71, pp. 116-131 (1987). However, such other techniques of pulse generation for MR imaging do not provide selective pulses that have periodic excitations needed for creating the selective stripes of SPAMM, and, in any event, such pulse sequences are generally too long to be used with SPAMM. Better pulses for use with SPAMM are thus not suggested in the prior art.

In addition, the frequency response of a general system such as SPAMM has been studied by one of the present inventors in a sequence of papers: S. M. Eleff et al., "The Synthesis of Pulse Sequences Yielding Arbitrary Magnetization Vectors", Journal of Magnetic Resonance in Medicine. Vol. 12, pp 74-80 (1989); M. Shinnar and J.S. Leigh, "Frequency Response of Soft Pulses", Journal of Magnetic Resonance, Vol. 75, pp. 502-505 (1987); M. Shinnar and J.S. Leigh, "The Application of Spinors to Pulse Synthesis and Analysis", Journal of Magnetic Resonance in Medicine, Vol. 12, pp. 93- 98 (1989); M. Shinnar et al., "The Use Of Finite Impulse Response Filters in Pulse Design", Journal of Magnetic Resonance in Medicine, Vol. 12, pp. 81-87 (1989); M. Shinnar et al., "The Synthesis of Pulse Sequences Yielding Arbitrary Symmetric Magnetization Vectors", in Journal of Magnetic Resonance. Vol. 72, pp. 298-306; and M. Shinnar et al., "The Synthesis of Soft Pulses With a Specified Frequency Response", Journal of Magnetic Resonance in Medicine, Vol. 12, pp. 88-92 (1989). It was therein shown that the frequency response of a series of pulses can be written as a Fourier series, whose coefficients are nonlinear functions of the pulse amplitudes . A method was also described whereby, if one specified the desired M_z magnetization as a Fourier series, one could generate a pulse sequence which will actually yield that magnetization. Furthermore, the way of specifying the desired M_z magnetization as a Fourier series was realized to be identical to a solved problem in the design of finite impulse response filters. Using this technique, which has been described in detail in the afore-mentioned co-pending patent application Serial Number 07/655,077, one could specify the desired magnetization in certain frequency ranges, the desired smoothness for each range, and then generate a pulse sequence which would optimally achieve that desired magnetization.

The present inventors have joined forces to investigate whether the afore-mentioned pulse generation technique could be applied to the design of a pulse sequence or sequences which would be better than the afore-mentioned binomial sequence for SPAMM imaging. Since the stripes laid down for the SPAMM technique correspond to setting the z magnetization M_z in certain frequency ranges, it seemed ideally suited for this pulse generation technique. As will be described herein, the present inventors have in fact found that improved images can be produced by using pulses generated in accordance with the afore-mentioned pulse generation technique in place of binomial pulses and that the analysis techniques developed for imaging of invasively implanted radiopaque markers to MR studies of heart wall motion using the noninvasive MR imaging techniques of SPAMM may be used to quantitatively measure movement within the myocardium of the heart.

SUMMARY OF THE INVENTION

In accordance with a first aspect of the present invention, methods are provided of constructing pulse sequences to selectively excite frequency bands in NMR imaging using the techniques of SPAMM. The characteristics of the pulses are preferably predetermined by the SPAMM user and used by the pulse generating system of the invention to generate the pulses which will yield the desired imaging stripes. Such "optimal" pulses have been found to produce improved NMR images, for "optimal," as used herein, means that the pulse sequence will have the smallest possible error in approximating the desired pulse.

In accordance with this aspect of the invention, a sequence of N hard pulses, P₁-P_N, is selected which satisfies the following conditions: (1) Each pulse is assumed to be instantaneous, or sufficiently short that its duration can be ignored.

(2) The same amount of time, t, occurs between any two pulses. During this time the magnetization precesses around the z axis with frequency ω.

(3) The total duration of the pulse sequence, Nt, is sufficiently short that relaxation can be ignored.

Then, if it is assumed that the z magnetization, as a function of frequency, started out at equilibrium, with M_z(ω)=1, then the resultant z magnetization, assuming the system started in equilibrium, can be written as an (N-1)th order Fourier series in ωt, where ω is the off-resonance frequency. Furthermore, if all pulses have the same phase, then the z magnetization is symmetric in frequency (M_z(ω) = M_z(-ω)), and can be written as an Nth order Fourier cosine series in ωt. This Fourier series is not a Fourier transform of the hard pulse sequence, and the relationship between the pulses and the response is non-linear.

In preferred embodiments, the pulse generating method of the invention has the following characteristics:

(1) If one is given a Fourier series (in ωt) description of the desired z magnetization of a hard pulse sequence of N pulses, then one can actually do an inversion of this nonlinear problem and determine a hard pulse sequence which will actually yield that desired response.

(2) Furthermore, if one is given a Fourier cosine series, one can generate a series of hard pulses around a fixed axis which gives the desired response.

(3) In general, one starts with a specification of the desired z magnetization not as a Fourier series in ωt, but rather as having certain desired values over several frequency ranges. It was realized by one of the inventors that the techniques of digital filter design, specifically that of finite impulse response filters, not part of the art of NMR, could be applied to yield the desired Fourier series. Because the theory of digital filter design is quite advanced, one can generate Fourier series which optimally give the desired response. One thus can generate hard pulse sequences which have the optimal frequency response.

The method in accordance with the first aspect of the invention thus provides a technique of synthesizing a sequence of hard pulses for generating any desired realizable magnetization. The sequence of hard pulses is generated around a fixed axis which will generate any desired realizable magnetization which is symmetric in frequency. Preferably, the techniques of finite impulse response filters are used for specifying the desired magnetization. In accordance with a preferred embodiment of the invention, the generated hard pulse sequences are used in the pre-imaging SPAMM sequence to develop improved image slices of a patient's heart.

The novel advantages, features and objects of the first aspect of the present invention are specifically realized by a method of detecting motion of an image slice of a body portion of a patient using a magnetic resonance imaging device, comprising the steps of:

inputting parameters specifying the desired visual characteristics of the image slice;

synthesizing a radio frequency pulse sequence which, when applied to the body portion prior to imaging by the magnetic resonance imaging device, will cause the creation of an image slice having the desired visual characteristics;

applying to the body portion an external magnetic field so as to produce a resultant magnetization;

applying a pre-imaging sequence to the body portion so as to spatially modulate transverse and longitudinal components of the resultant magnetization to thereby produce a spatial modulation pattern, said pre-imaging sequence comprising the synthesized radio frequency pulse sequence and magnetic field gradient pulses;

applying an imaging sequence to the body portion to make visible said spatial modulation pattern of intensity; and

detecting motion of the body portion during a time interval between application of the pre-imaging and imaging sequences by examining displacements of the spatial modulation pattern of interest occurring during the time interval.

In accordance with a second aspect of the invention, the techniques developed for motion analysis from radiologic imaging of invasively implanted metallic myocardial markers are adapted to the analysis of tagged MR images of heart wall motion. The apparatus of the invention may be incorporated into existing MR imaging devices by way of a unit specially designed for this purpose, or a computer program may be implemented to provide motion analysis for 2-D or 3-D motion in accordance with this aspect of the invention so as to produce a visual display of the results. The applications include, but are not limited to, analysis of heart wall motion as described herein.

The novel advantages, features and objects of this aspect of the invention are realized by a noninvasive method of determining displacements of a body portion of a patient during a given time interval, comprising the steps of:

applying to the body portion an external magnetic field so as to produce a resultant magnetization;

applying a pre-imaging sequence to the body portion so as to divide the body portion into a grid of intersecting stripes having altered magnetization, the stripes intersecting at a plurality of intersecting points;

applying an imaging sequence to the body portion; quantitatively measuring physical displacements of predetermined ones of the plurality of intersecting points during the given time interval; and

providing an image representing the displacements. In a preferred embodiment of the method of this aspect of the invention, the pre-imaging sequence applying step comprises the step of specifying a triangular grid of the body portion wherein the predetermined ones of the plurality of intersecting points are respective vertices of respective triangles of the triangular grid. Also, the displacements measuring step may comprise the step of determining which of the intersecting points in the triangular grid at the beginning of the given time interval correspond to the predetermined ones of the plurality of intersecting points in a deformed triangular grid at the end of the given time interval. The displacements measuring step may also comprise the step of comparing the triangular grid at the beginning of the given time interval with the deformed triangular grid at the end of the given time interval. Accordingly, the image representing the displacements may comprise a comparison of an image representing the triangular grid at the beginning of the given time interval with an image representing the deformed triangular grid at the end of the given time interval. Preferably the body portion is the patient's heart wall.

BRIEF DESCRIPTION OF THE DRAWINGS

The above and other objects and advantages of the invention will become more apparent and more readily appreciated from the following detailed description of the presently preferred exemplary embodiment of the invention taken in conjunction with the accompanying drawings, of which:

FIGS. 1 (a)-(d) are timing diagrams of a binomial

SPAMM pulse sequence for a two-dimensional grid formation prior to imaging.

FIG. 2 is a schematic diagram of a magnetic resonance imaging system in accordance with an embodiment of the invention whereby optimal rf pulses are calculated for application during SPAMM imaging.

FIG. 3 is a flow chart illustrating a preferred embodiment of the SPAMM imaging technique of the invention.

FIG. 4 is a diagram illustrating the effect of choosing a narrower stripe on the smoothness of the interstripe region.

FIG. 5 is a diagram illustrating the stripe differences between a seven pulse binomial sequence and a seven pulse optimal pulse sequence generated in accordance with the techniques of the invention. FIG. 6 is a schematic representation of a 2-D transformation of a triplet of points (R1-R3) in the myocardium from an initial state (left) to a later state (right).

FIG. 7 is a schematic representation of a 3-D transformation of a quadruplet of points (R1-R4) in the myocardium from an initial state (left) to a later state (right).

FIG. 8 illustrates an overall flow chart for a program for analyzing motion from "tagged" MR images in accordance with the invention.

FIG. 9 illustrates two-dimensional tagged MR short axis images of the heart of a human immediately after tagging at the end of diastole (left) and at the same level in late systole (right).

FIG. 10 illustrates the points specified for stripe intersections for images in FIG. 9.

FIG. 11 illustrates the displacement of specified points between the times of the taking of two images.

FIG. 12 illustrates the corresponding triangulation for the points shown in FIG. 9.

FIG. 13 illustrates the magnitude of mean displacement of triangles between the imaging times of FIG. 9.

FIG. 14 illustrates eigenvectors of the transformation of triangles between the imaging times of FIG. 9, represented as a symbol with an initial unit circle and superimposed major and minor axes of a corresponding subsequent ellipse into which it would be transformed.

FIG. 15 illustrates the angle of the principal eigenvectors for the triangles of FIG. 14.

FIG. 16 illustrates a combined display of point displacements (vectors) and the magnitude of mean displacement of triangles derived from long axis views of the same subject and at the same times as in FIG. 9. DETAILED DESCRIPTION OF THE PRESENTLY PREFERRED EMBODIMENTS

A SPAMM imaging technique in accordance with presently preferred exemplary embodiments of the invention will be described below with reference to FIGS. 1-16. It will be appreciated by those of ordinary skill in the art that the description given herein with respect to those figures is for exemplary purposes only and is not intended in any way to limit the scope of the invention. All questions regarding the scope of the invention may be resolved by referring to the appended claims.

The SPAMM technique may be used to image the heart by simply using a three pulse sequence consisting of two nonselective RF excitation pulses separated by a constant amplitude magnetic field gradient pulse. The sequence is then followed after a predetermined delay by a conventional imaging sequence. The first RF pulse produces transverse magnetization which is initially phase synchronized. The first RF pulse is followed by a gradient pulse which produces a regular variation of the phase of the excited spins along the direction of the gradient, with a wavelength inversely proportional to the integral of the gradient pulse amplitude. The second RF pulse turns some of the phase-modulated transverse magnetization back into longitudinal magnetization with a corresponding modulation. The modulation takes the form of a sinusoidal dependence of the longitudinal magnetization on position along the direction of the modulating gradient. The amplitude of the modulation depends on the strength of the RF pulses. For example, two 45° pulses will produce modulation between the extremes of saturation and full magnetization. However, longitudinal relaxation between the SPAMM and imaging sequences reduces the amplitude of the modulation. Accordingly, the pulse angles of the applied radio frequency pulses may sum to more than 90° to account for the longitudinal relaxation.

In the subsequent imaging sequence, the modulated longitudinal magnetization results in a corresponding intensity modulation in the final image, appearing as regularly spaced stripes perpendicular to the direction of the modulating gradient. Then, since the magnetization of any material moves as the material moves, any motion of the material between the time of the initial SPAMM sequence and the subsequent imaging sequence will be reflected in the corresponding displacement of the stripes in the final image.

SPAMM is analogous to some techniques used to produce selective excitation of chemical shifts in nuclear

MR spectroscopy (e.g., for solvent suppression). In particular, the use of equal-strength nonselective RF pulses to produce selective saturation is roughly equivalent to the

1-1 case of binomial excitation pulses in spectroscopy as described by Hore in "Solvent Suppression in Fourier Transform

Nuclear Magnetic Resonance", Journal of Magnetic Resonance, Vol. 55, 1983, pp. 283-300, with the modulating gradient pulse taking the place of the chemical shift.

However, improved images have been obtained when more than two radio frequency pulses separated by gradients are used to generate the SPAMM stripes. For example, the selective excitation binomial pulse sequence has been shown by Axel et al. in the article entitled "Heart Wall Motion: Improved Method of Spatial Modulation of Magnetization for MR Imaging", Radiology. Vol. 172, 1989, pp. 349-350, to result in improved SPAMM heart wall images. In particular, Axel et al. therein disclose the use of a pulse sequence with a train of nonselective RF pulses, the relative amplitudes of which are distributed according to binomial coefficients such as 1-3-3-1 or 1-4-6-4-1, for example. The absolute amplitudes are adjusted for the total amplitude of modulation desired. Separating each of the RF pulses is a gradient pulse oriented perpendicularly to the desired stripe orientation and with a strength and duration dependent upon the desired stripe spacing. For the gradient amplitude G, the stripe spacing will be equal to the reciprocal of γ ∫ Gdt, where η is the gyromagnetic ratio.

FIGS. 1 (a)-(d) illustrate an example of a binomial SPAMM pulse sequence for two-dimensional grid formation prior to imaging, with RF pulse amplitudes distributed as 1-4-6-4-1. Preferably, there is a phase shift of 90° from RFI to RFQ to avoid artifacts, while gradients G1 and G2 are orthogonal so as to produce orthogonal sets of stripes. Gradient pulses used for modulation have an amplitude and duration such that the area under the pulse shape is inversely proportional to the desired stripe spacing. A brief delay between gradient pulses and RF pulses allows the gradients to settle. Also, to avoid perturbing the stripe pattern with the effects of local field variations or regional chemical shift differences (e.g., fat vs. water), the time interval between RF pulses should be kept short.

Thus, a two-dimensional array of stripes may be used to form a tagging grid by following the initial SPAMM sequence with a second one, with the gradient oriented in an appropriate new direction. To avoid producing artifacts from stimulated echo formation, the phase of the RF pulse is preferably shifted by at least 90" between the two sequences. For effective composite flip angles other than 90°, the degree of saturation in the stripe intersections may differ from that in the stripes.

The SPAMM sequencing for forming a two-dimensional array of stripes which forms a tagging grid as described above may be implemented on a conventional MR imaging system such as the Signa 1.5-T MR imaging system from GE Medical Systems. A simplified schematic diagram of such an MR imaging system is shown in FIG. 2. As shown, the MR imaging system has a main magnet 200 which provides a steady magnetic field B₀ which polarizes the nuclei of the protons of the specimen or subject for which an image is desired. Generally, within magnet 200 there is a space in which the specimen or individual to be examined is placed.

The MR imaging system of FIG. 2 also includes a set of three orthogonal direct current coils 202, 204, and 206 for producing spatial linear field gradients G_x, G_y, and G_z. These coils are driven by gradient generator 208, which in turn is controlled by a controller 210 which communicates with a host computer 212. Host computer 212 generates the pulse sequences and controls application of the field gradients. Host computer 212 thus implements the pulse generation and displacement measurement processes of the invention to be described in more detail below.

The MR imaging system of FIG. 2 also includes a radio frequency (RF) coil 214 which generates a radio frequency field in the specimen being analyzed and also senses a free induction decay or spin echo signal which is generated after termination of the radio frequency pulse. RF pulse generator 216 excites RF coil 214 in response to host computer 212 to apply the desired pulse sequences. The signal processor 218 receives the small microvoltage level spin echo signals, converts them into digital form and inputs them into host computer 212 for formation of an image (typically using Fourier transform techniques). The resulting images are then manipulated for display and stored in storage device 220 for later retrieval or immediately displayed on display device 222.

Hence, the apparatus of FIG. 2 may be used to implement the pre-imaging SPAMM sequence and to implement the imaging sequence for detecting heart wall movement. In a preferred embodiment, SPAMM pulses for heart in..ging may be integrated into a conventional cardiac-synchronized imaging sequence, where the SPAMM sequence can be started with a trigger pulse derived from the electrocardiogram. A two-dimensional grid of stripes can then be produced with two 1-4-6-4-1 RF pulse sequences, for example, along with gradient pulses. Although binomial sequences produce excellent SPAMM stripes, the present inventors have now found that "optimal" pulses may be obtained by allowing the SPAMM user to specify the desired stripe parameters and then generating the pulse sequence which will produce the desired stripes. Accordingly, this pulse generation technique will now be described in detail with respect to FIGS. 3-5.

FIG. 3 is a flow chart illustrating a preferred embodiment of the SPAMM imaging technique of the invention. As shown, the pulse generating system of the invention may be implemented as a process which runs on host computer 212 (FIG. 2). A preferred embodiment of such a process is given below as Appendix A and should be referred to in conjunction with the description of FIG. 3. However, before describing the pulse synthesis algorithm with respect to FIG. 3, it will first be explained how the units of the pulse synthesis algorithm correspond to the units of SPAMM imaging.

In particular, the algorithm to synthesize the pulse sequence specifies the desired magnetization in terms of frequency units and the duration of the portion of the pulse sequence during which the spins are freely processing. In MR imaging, on the other hand, one is interested in the magnetization at a desired location. The relationship between frequency and location is adjustable, depending on the gradient strength. However, one can easily convert the units so that one can see the relationship. In other words, since the magnetization response of such a hard pulse sequence is periodic, with a period 1/t, where t is the time the spins freely precess between any two pulses (not the total pulse sequence duration), this periodicity may be used to lay down a series of stripes.

The pulse synthesis algorithm of the invention allows one to specify a series of frequency bands, with the desired magnetization and desired smoothness being specified in each band. The specification of the frequency bands is most naturally done in terms of the periodic frequency, 1/t. For example, one may choose to have the z magnetizations of spins with frequencies between 0 and .1/t be 0 and those between .2/t and .8/t be 1. However, because the pulses are all around the same axis, one is guaranteed symmetry around zero frequency and therefore around .5/t. That is, Mz(ω) = Mz(-ω), where ω is the frequency.

Thus, to convert the above to imaging, it becomes apparent that this approach automatically specifies the stripe thickness in terms of the interstripe distance. The interstripe distance is set by determining how long the gradient will be on between pulses, which determines the interstripe frequency difference, and then the gradient strength, which determines the conversion of frequency to distance. Accordingly, the interstripe distance is ultimately determined by the product of the average gradient strength (H_z/cm) and the time the gradient is on between pulses. Increasing the gradient strength or the time will decrease the interstripe distance.

One skilled in the art will appreciate that the above discussion assumes the gradient is off during the application of the radio frequency pulses. This guarantees that, subject to the uniformity of the radio frequency field of the coil, the pulses are sufficiently hard to excite all spins equally. However, one may want to leave the gradient on during the pulses. This requires that enough radio frequency power be used that the pulses are sufficiently

"hard" to equally excite all frequencies of interest. For example, if the radio frequency pulses are very narrow, they must have higher amplitudes for a given gradient strength. Otherwise the procedure is unchanged from the approach where the gradient is off when the radio frequency pulses are applied.

As shown in FIG. 3, MR imaging starts at step 300 in accordance with the invention by setting the value of internal program constants (Appendix A lines 1-14) (step 302) and then requesting the user to input the desired pulse parameter data (Appendix A lines 15-57) (step 304). As noted above, this data generally includes the number of desired RF pulses to be used in the pre-imaging sequence, the desired stripe thickness and the size of the transition zone between the striped and non-striped regions. The Fourier cosine series which optimally fits the input specifications is then calculated (Appendix A lines 58-108) (step 306) using an algorithm similar to that previously disclosed by T.W. Parks et al. in an article entitled "A Program for the Design of

Linear Phase Finite Impulse Response Filters," IEEE

Transactions on Audio and Electroacoustics. Vol. AU-20, pp. 195-199 (1972). The mathematical theory behind application of this technique for pulse generation can be found in the afore-mentioned related U.S. Patent Application Serial No. 07/655,077 and is specifically incorporated herein by reference. The resulting Fourier cosine series can then be renormalized (Appendix A lines 109-130) (step 308) so that the absolute value of the Fourier cosine series is always less than one.

Once the Fourier cosine series is renormalized, it is then converted to a complex exponential Fourier series (Appendix A lines 135-142) (step 310). The resulting data is then formatted so that it can used by a pulse generating routine, or the data is fed directly into a pulse generating subroutine for synthesizing an optimal pulse sequence (Appendix A lines 131-143 and subroutine "PULS"). A Fourier series for the spinor components ((1+M_z)/2 = ¦p¦² and (1-

M_z)/2 = ¦Q¦² is then calculated from the Fourier series for M_z

(Appendix A subroutine "PULS" lines 16-38) (step 312). P and

Q are then calculated using cepstral deconvolution (Appendix A subroutine "PULS" lines 39-53) (step 314). For each set of P and Q, the pulse sequence which yields the spinor with P and Q as its spinor components is then calculated (Appendix A subroutine "PULS" lines 54-88) (step 316). This latter step is generally accomplished by attempting to reduce the Fourier series order one order at a time as described in detail in the afore-mentioned related U.S. Patent Application Serial No. 07/655,077 and then choosing which values of P and Q to use. The result will be the optimal RF pulse sequence which can then be output to a file or other media for storage (step 318).

Once the optimal pulse sequence has been generated, the pre-imaging and imaging sequences may be started. This is accomplished by applying the external magnetic field B₀ to the specimen (step 320) and feeding the optimal pulse sequence to the RF coil for application with the gradient pulses as the SPAMM pre-imaging sequence (step 322). The imaging sequence is then applied to detect the stripe displacements (step 324). Further details regarding the optimal pulse generating technique of the invention can be found upon closer inspection of the Fortran source code listing of Appendix A.

The disclosed pulse generating technique of the present invention thus obtains the B₁ field which will actually yield a desired excitation profile. If the desired excitation profile is physically realizable, the method of the invention will invert the nonlinear relationship between the input B₁ field and excitation to obtain the desired B₁ field. Moreover, if one starts out with a desired excitation profile which is not physically realizable, the method of the invention gives a procedure for a physically realizable profile which is closest to the desired profile, and then finds a B₁ which will yield the desired profile.

The following represents the mathematical formulation of a pulse sequence which will yield arbitrary symmetric magnetization vectors for use as RF pulses in SPAMM imaging in accordance with the invention.

As described by the one of the present inventors in an article entitled "The Synthesis of Pulse Sequences Yielding Arbitrary Symmetric Magnetization Vectors", in Journal of Magnetic Resonance, Vol. 72, pp. 298-306, the first step of the pulse generation method of the invention is to start with a physically realizable z magnetization which may be achieved by a hard pulse sequence and to determine a hard pulse sequence which will actually yield that z magnetization. In particular, an N hard pulse sequence to be applied around the x axis within a total duration of T has a magnetization response M_z which can be written (as a function of frequency) as a Fourier cosine series:

N-1

M_z (ω) = Σ a_jcos(jωT/(N-1))

j=0

Equation (1) with |M_z(ω)I always ≤1.

If a specification of a desired M_z(ω) is given, then we know, since magnetization is preserved by the pulses, that M² _x(ω) + M² _y(ω) + M² _z(ω) = 1. The first step of the invention is thus to find an appropriate solution of M_x and M_y which will satisfy this equation. For this purpose, we write M_xy(ω) = M_y(ω) - il%, (ω) . The above condition then becomes |M_xy(ω)|² + M² _z(ω) = 1.

This reduces to a phase retrieval or deconvolution problem, well studied in signal analysis, although generally not familiar to those skilled in the NMR art. There are a number of approaches to solving this problem. For example, one can convert M_z(ω) to a complex polynomial in s = exp(iωT/(N-1)) and then solve for the roots of the polynomial 1 - M² _z(ω) , or equivalently for the two polynomials (1 - M_z(ω), 1 + M_z(ω)). One then groups the roots of the polynomials by means of symmetry considerations and chooses half of them. There are also techniques such as cepstral deconvolution which allow for the solution of M_xy(ω) as an (N-1) degree polynomial in s = exp(iωT/(N-1)). This is a representation of M_xy(ω) as an N-1th degree complex Fourier series. One can guarantee that since the polynomial for M_z had real coefficients, the polynomial for M_xy will also have real coefficients. The Fourier coefficients may be generated by viewing M_z as a function of Ω, where Ω = ωT/(N-1), using a digital fast Fourier transform and looking at the cosine coefficients.

The following are the details of one possible approach to the resulting deconvolution problem:

Letting Ω=ωT/(N-1), we first express |M_xy(Ω)|² as a polynomial in cos(Ω). Namely, using Tschebyscheff polynomials, one can write cos(nΩ) = T_n(cos(Ω)) as a polynomial in cos(Ω). Thus, one can write:

N

M_z(Ω) = Σ b_j{cos(Ω))^j=Q(cos(Ω)).

j=0

Equation (2) We now assume that there is no loss of magnetization during the entire pulse sequence. That is, N_t<<T₂ for N pulses with a constant time spacing t. Then:

|M_z(Ω)|²+|M_xy(Ω)|²=|M_o|²,

Equation (3) where M_o is a constant which may be normalized to 1. We can therefore write M_xy(Ω) = sqrt(1-|M_z(Ω)|²)eⁱ'^ø, where ø is the phase relationship we are trying to construct. We shall use this relationship to write M_xy(Ω) as a complex Fourier series.

Since M_Z(Ω) is real, |M_Z(Ω)|² = (M_Z(Ω))². We can therefore write, for u=cos(Ω):

2N

P(u)=1-Q(u)²=1-|M_z(Ω)|²=Σ b_ju^j=|M_xy(Ω)|²

j=0

Equation (4)

Viewing P(u) as a polynomial in u, we now factor P(u) to solve for its roots:

2N

P(u) = b_2NII(u-u_j)

j=1

Equation (5)

We now define s=e^iΩ, where u=(s+s^-1)/2. We can therefore define rational function R(s) and give its factorization as follows: 4N

R(S) = P((s+s^-1)/2) = b_2N2-^2Ns^-2NII(s-s_j) = 1-|M_z(Ω)|² j=1

Equation )

Next, we group the roots into conjugate pairs noting that S_j is a root of R(s) if and only if (s_j+s_j ^-1)/2=u_k, for some k. Thus, if S_j is a root, s_j ^-1 is a root. Furthermore, because P(u) has real coefficients, if U_j is a root, u_j* is a root (with * denoting complex conjugate). So if s_j is a root, S_j ^* is a root as well. We also note that |s|=1 occurs only when u is real, and |u|≤1. Since P(u)=1-|M_z(Ω) |²≥0 for u=cos (ωt), thus any root U_j with u_j real and |u_d|<1 must occur with even multiplicity. Therefore, roots s_j can be grouped as follows:

1) r_j real, |r_j| = 1 can be grouped into r_j , r_j ^-1,

j=1,L.

2) f_j complex, |f_j| = 1. Then f_j, f_j ^-1, f_j ^-1* are

roots, j=1,M.

3) g_j real, Ig_j| = 1. Every g_j occurs with even

multiplicity, 2*K_j, j = 1, K.

4) h_j complex, |h_j| = 1. Again, h_j occurs with even multiplicity, 2*H_j. Furthermore, if h_j is a root, so is h_j ^* = h_j ^-1. Thus, we can group roots into h_j, h_j ^-1, with multiplicity 2*H_j, j = 1,H.

We next need to select a phase function for M_xy(Ω). Accordingly, for each set of roots, we choose a representative. Thus, for each set of real roots a_d we choose either r_j of r_j ^-1. For each set of complex roots, f_j , we choose a conjugate pair, either f_j, f_j ^*, of f_j ^-1, f_j ^-1*. This gives us ₂ ^(L+M) choices. Each choice gives us a different phase function for M_xy(Ω).

Using our choices in the previous step, we now define the following rational function to define M-^Ω) :

L M

G(s) = As^-N(II (s-r_j)II (s-f_j) (s-f_j ^*)) x

j=1 j=1

K H

(II (s-g_j)K_j II (s-h_j) (s-h_j ^-1).

j=1 j=1

Equation (7) L H K

with A² = b_2N2-^2N/ (II (r_j) II (f_jf_j*) II (h_j) K_j)

j=1 j=1 j=1

Equation (8) It can be shown that A is real, so G(s) has real coefficients. We now claim that R(s)=G(s)G(s^-1), which can be easily verified by multiplying through. We define:

M_xy(Ω)=G(s), for s=e^iΩ.

Equation (9) We claim that this M_xy(Ω) satisfies Equation (3), and since G(s) has real coefficients, G(s^*)=G(s)^*. Therefore, for |s|=1,

|M_xy (Ω)|²= G(S)G(S)^* = G(S)G(S^*) = G(S)G(S^-1) = R(S) =

1-|M_z(Ω)I².

Equation (10)

In defining G, we could choose either A or -A. Thus, there are 2^L+M+1 choices in defining G. M_Z(Ω) and M^fΩ) are thus expressed as finite Fourier series and the xy magnetization after the jth pulse is defined, M_xy(Ω), for j=0 and j=N. The next step is to find a pulse sequence which will yield that Fourier series.

One can express the effect of a pulse of f degrees around the x axis on the magnetization (M_x, M_y, M_z) as a matrix P(f). The effect of free precession for a time (T/(N-1)) at frequency ω can also be written as a matrix R(ω). Thus, the effect on the magnetization of a free precession and a pulse is (P(f) ·R), viewed as a matrix product.

We then apply R^-1P^-1(f) to the given (M_x, M_y, M_z). The resultant expression is still a Fourier series. We then solve for the f that will yield a Fourier series of order N- 2. We then repeat the procedure until we reduce the magnetization to a constant M_z and constant M_xy. We then apply a P^-1(f) which will yield M_z=1 and M_xy=0. The resultant series of pulses will yield, when applied in reverse to a system with M₂=1 and M_xy=0, a system with magnetization equal to the desired one.

That this part of the invention could be done was announced in an abstract published by one of the inventors entitled "An Exact Synthesis Procedure For Frequency Selective Pulses", in 5th Proceedings of the Magnetic Resonance in Medicine, pp. 1452-53, August, 1986. However, the details of how to obtain the pulse sequence were not given in that abstract. The technique for obtaining the pulse sequence first appeared in the aforementioned articled entitled "The Synthesis of Pulse Sequences Yielding Arbitrary Symmetric Magnetization Vectors" by one of the present inventors. The relevant portions of the calculations set forth in that paper for obtaining the hard pulse sequence will be reproduced here.

Having generated M_xy(Ω) , we have defined the xy magnetization after the jth pulse, M_xy(Ω)_j, for j=0 and j=N. We will develop a simple recursion relationship to allow the calculation of M_xy(Ω)_j-1 from M_xy(Ω)_j. That is, we will work backwards. In addition, a formula for calculating the jth pulse from M_xy(Ω)_j will be derived. Having accomplished these two substeps, the pulse sequence will be completely defined.

The following matrix definitions will be used for the derivations.

M = [m_x,m_y,m_z]^t , with ^τ denoting the transpose, and

R = [cosΩ -sinΩ 0] P = [1 0 0 ]

[ sinΩ cosΩ 0 ] [ 0 cosø -sintø ]

[ 0 0 1] [0 sinø cosø ] where M is the magnetization vector, P·M simulates a pulse and flips the spins around the x axis, and R·M rotates the spins in the xy plane. The sequence we will use applies a rotation and then a pulse. In matrix notation this is (P·R). The initial rotation with M=[0,0,1] has no effect, and the final pulse is not accompanied by a rotation since all spins are pre-rotated (that is, the operation is (P·R) not (R·P)). Spins are neither created nor destroyed during the application of the pulse sequence. Removing the effects of (P·R) from M_p will yield M_p-1.

This is accomplished by multiplying by (P·R) inverse,

(P·R)^-1:

Hence,

M_p-1 = {(P·R)_p ^-1} M_p = {(P·R)_p ^-1} {[m_x,m_y,m_z]_p ^t}

Equation (11)

The following derivations hinge on the fact that the only perturbations on the spin system are repeated applications of precession around the z-axis followed by a pulse around the x-axis. Rotation around the y-axis is not covered by the resulting formulae. Noting that M_xy = -M_y + iM_x, and substituting into Equation (11) yields:

[M_z]_p-1 = -sinø_p[M_y]_p + cosø_p[M_z]_p or

[M_z]_p-1 = (1/2)sinø_p[M_xy + M_xy*]_p + [M_z]_pcosø_p.

Equation (12)

Similarly,

[M_xy]_p-1 = [-M_y + iM_x]_p-1 = [M_xsinΩ - M_ycosøcosΩ

-M_zsinøcosΩ + iM_xcosΩ + iM_ycosøsinΩ + iM_zsinΩsinΩ]_p =

[M_x]_p(icosΩ + sinΩ) + [M_y]_pcosø(-cosΩ + isinΩ) + [M_z]_psinø(-cosΩ + isinΩ) = [M_x]_p (ie^-iΩ) + (M_x]_pcosø(-e^-iΩ) + [M_x]_psinø(-e^-iΩ) =

[M_z]_psinø(-e^-iΩ) = [M_z]_psinø(-e^-iΩ) + ( 1/2 ) e^-iΩ{ ( [M_xy]_p +

[M_xy*]_p)cosø+ ( [M_xy]_p- [M_xy*]_p) } . Thus ,

[M_xy]_p-1 = (1/2 ) (1+cosø_p) e^-iΩ[M_xy] _p - (1/2) (1-cosø_p) e^-iΩ [M_xy*]_p - sin ( ø_p )e^-iΩ[M_z]_p .

Equation ( 13 )

Numerically, we have [M_z]_N and [M_xy]_N as Fourier coefficients. We therefore need to convert Equations (12) and (13) into equivalent form. This can be accomplished by equating terms of equal order. Let a_r represent the rth Fourier coefficient for M_z which corresponds to e^irΩ and let f_r correspond to the rth Fourier coefficient for M_xy. We note that since:

M_z(Ω) = M_x(-Ω) and a_r = a_-r, that:

[a_r]_p-1 = (1/2)sinø_p([f_r]_p + [f_-r]_p) + cosø_p[a_r]_p

Equation (14) and

[f_r]_p-1 = ( 1/2 ) ( 1+cosøp) [ft_r+1] _p - (1/2) (1-cosø_p) [f-_(r+1)]_p-sinø_p[a_r+1]_p.

Equation (15)

For a sequence of N pulses, the Fourier coefficients range from -(N-1) to (N-1). That is, [f_r]_p = [a_r]_p = 0 for |r|≥p. We now set r = (p-1) in Equation (14) and get [a_p-1]_p-1 = 0 =sinø_p([f_p-1]_p + [f_-(p-1)]_p + cosø_p [a_p-1]_p. Rearranging yields: tanø_p = -[a_p-1]_p ([f_p-1]_p + [f-_(P-1)]_p).

Equation (16)

Setting r = -p in Equation (15) yields [f-_p]_p-1 = 0 = (1/2)(1+cosø_p)[f-_p+1]_p-(1/2(1-cosø_p)[f_-(p+1)]_p - sinø_p[a_-p+1]_p.

Since we know that a_r=a_-r, it can be shown that [f_p-1]_p and [f_1-p]_pare of opposite sign. Therefore, we can define: cosø_p = ( [f_p-1]_p + [f_-(p-1)]_p) / ( [f_p-1]_p " [f-_(p-1)]_p) .

Equation (17)

Substituting into Equation (16) we obtain the last of the recursion values:

sinø_p = -[a_p-1]_p/([f_p-1]_p - [f-_(p-1)]_p).

Equation (18)

The synthesis of the pulse sequence is completed by using the M_z Fourier coefficients, [a_r]_N, where M_z(Ω) is generated using a digital fast Fourier transform and from the

M_xy Fourier coefficients, and by using [f_r]_N from the generation of M_xy (Ω) to calculate the leading flip angles cosø_n and sinø_n from Equations (17) and (18) with p = N. We then decrease p, the current pulse, which started at N, and calculate [a]_p, [f]_p, cosø_p, and sinø_p for p = N-1, ..., 2. We can then calculate the first pulse by calculating the net angle required at M_z(Ω=0) and the sum of all angles calculated in the previous two steps, i.e., if M_z(Ω=0) = 0, and the sum of the ø_p' s for p=2,N was π/3, then the first pulse angle should be -π/3.

The above procedure can also be generalized to cases where M_z was specified as a complex Fourier series. In such a case, the above procedure could be repeated. However, the phase retrieval problem is more complicated, and the coefficients of M^ will in general be complex. One may then repeat the above inversion procedure. However, then a pulse is specified not only by the number of degrees it rotates, but also by the phase of the pulse in the xy plane. A solution of the inverse problem was published by one of the present inventors in an article entitled "The Synthesis of Pulse Sequences Yielding Arbitrary Magnetization Vectors", Journal of Magnetic Resonance in Medicine, Vol. 12, pp 74-80 (1989), the contents of which are hereby incorporated by reference. Hence, the present invention may be easily extended to cases where the magnetization does not correspond to a fixed axis.

As noted above with respect to FIG. 3, the preferred embodiment of the pulse generator of the invention uses spinors in the optimal pulse determination. Spinors allow the rotation operator to be monitored, as described, for example, by L.D. Landau and E.M. Lifschitz in Quantum Mechanics: Non- Relativistic Theory, pp. 188-196, Pergamon Press, London, 1965, and by I.V. Aleksandrov in The Theory of Nuclear Magnetic Resonance, pp. 169-181, Academic Press, NY 1966, and are thus presently preferred. Spinors are useful when one wishes to specify not only the final z-magnetization, but also the rotation operator. For example, in imaging, one may wish to specify a pulse that rotates 180° about an axis in the xy plane, like an inversion pulse. Furthermore, one may wish to specify that all frequencies get rotated about the same axis. This type of pulse is a refocusing pulse. Therefore, looking at the pulse as an operator on the magnetization, one would like to specify the components of this operator. The present invention may be extended to this area as well using spinors. Such an extension using spinors has been described by one of the present inventors in an abstract in August 1988 and later in more detail in a paper entitled "The Application of Spinors to Pulse Synthesis and Analysis", Journal of Magnetic Resonance in Medicine. Vol. 12, pp. 93-98 (1989), the contents of which are hereby incorporated by reference.

The present invention has been described as a system which generates "optimal" RF pulses for use in generating SPAMM images. As noted above, the RF pulses applied to the RF coil if they are the RF pulses requested by the operator of the MR imaging device. Accordingly, some comments regarding the features of such pulses will be given here to help the operator determine what pulse features may be desirable.

1) Interstripe distance. One wants to have stripes that are sufficiently close so that one can track many different points. However, the stripes should be sufficiently far apart so they can be easily separated. The major import of this interstripe distance is that if one can give narrower stripes, one might be able to space them closer together.

2) Narrowness. The stripes laid down by the pulses should be narrow. However, they should be at least a pixel wide, maybe even wider so that they can be easily seen on the image. Since the pixel size may change with the imaging technique and the gradient strength applied, the optimal stripe thickness may vary from specimen to specimen. This condition is easily adjustable by using the pulse generating algorithm of the invention to specify the line width.

3) Sharpness. The transition from stripe to normal image should be as sharp as possible, so as little of the picture is disturbed as is necessary. This is again a parameter which can be specified with the pulse generating algorithm of the invention. This parameter can be specified as a fraction of the interstripe distance. Also, the observed narrowness of the stripes is actually a function of the sharpness. That is, the stripe is observable not only on the band where it is specified, but also over part of the transition zone, where the z magnetization is between the specified value and 1. As long as the z magnetization is sufficiently small, it will be observed as a dark line.

Thus, for the sharpest lines, one can specify an extremely narrow line and use the transition zone to make it sufficiently wide so as to be visible.

4) Flatness. Away from the stripes, the image intensity should be minimally disturbed, so as not to interfere with the interpretation of the image. This can be specified in the pulse generation algorithm of the invention. Also, as shown in FIG. 4, there may be a tradeoff between flatness and sharpness, as increasingly sharp transitions lead to ripples. In other words, as shown in FIG. 4, the flatter pulse 400 is less sharp than the rippled pulse 402. The effect of a sharp transition can be modeled, so that one can choose an appropriate tradeoff.

5) Persistence. The stripes laid down should persist through the cardiac cycle, at least to the end of systole. This may be accomplished by using a flip angle of more than 90° to generate the stripes. This will lead to the stripes initially having a white center, as the middle is inverted more than 90°. The white stripe gradually disappears in the cardiac cycle. For example, a flip angle of 120° after a short delay will provide the same stripe as a 90° flip angle with a lesser delay.

6) Determination of the center. One would like to be able to measure the center of the stripe w-th some precision. The presence of a white line in the middle may help in this determination.

7) Duration of the pulse sequence. The pulse sequence should not be too long. The precise determination of what too long means can change, but one would like to devote not more than roughly 30 msecs to generating the stripes. That way one can lay them down at end diastole, with minimal cardiac motion during the generation of the stripes. In order to apply the pulses, one has to give the pulse, then cycle the gradients on and off, which takes roughly 2 msecs. This means that given limitations of the equipment, and that one is generating two orthogonal sets of stripes, one cannot use more than a 7 pulse sequence. Sometimes, even that can be too long, and one has to use a shorter sequence. However, one skilled in the art will appreciate that an instrument which can be cycled faster could use even longer pulse sequences. As shown in FIG. 4, if one could use more pulses, one could make even sharper stripes.

Sequences of 5 and 7 pulse trains were synthesized using the technique described herein for different frequency constraints. The pulse trains were then implemented on a GE Signa Research magnet as the encoding pulse train of a 1 dimensional SPAMM experiment, and applied to a copper sulfate phantom. The gradient pulses were applied between each pulse as in FIG. 1, and a dephasing gradient was applied at the end of the pulse train. An imaging sequence was then done. The images were photographed, their intensities measured, and the stripe width measured. The stripes resulting from one such pulse sequence are shown by way of example as curve 500 in FIG. 5. Curve 500 was generated as a seven pulse sequence having the following respective flip angles (summing to 120°) : 12, 15.5, 20.7, 22.8, 20.8, 15.6 and 12.6. For comparison purposes, the stripes which resulted when a seven pulse binomial pulse sequence having the following respective flip angles (also summing to 120°): 1.9, 11.3, 28.1, 37.5, 28.1, 11.3 and 1.9, was applied to the same arrangement are shown in FIG. 5 as curve 502. From FIG. 5, it can be readily appreciated that the synthesized pulses were substantially better than the binomial pulses (i.e., sharp and flat).

Thus, the optimized pulses of the invention can significantly improve the SPAMM technique for cardiac imaging. Also, because the pulse values can be adjusted, a large family of pulse sequences with fairly similar performance can be generated by simply adjusting the relative weightings of the narrowness, sharpness, and flatness criteria. The operator is also given the opportunity to specify how much ripple or how little sharpness is tolerable. Furthermore, for any given z magnetization, one can synthesize many different pulse sequences (up to 2^N-1, where N is the number of pulses) which yield the same z magnetization, depending on root selection in the phase retrieval part of the pulse generation algorithm. This can also be done independent of the gradient strength.

The ability of the MR imaging techniques using optimal pulses in conjunction with periodic spatial modulation of magnetization (SPAMM) to define regional motion within the heart wall opens up the possibility of adapting the aforementioned analysis techniques for imaging of invasively implanted radiopaque markers to the analysis of MRI studies of heart wall motion. In particular, the techniques of SPAMM may be used in accordance with the invention prior to imaging to "tag" the heart wall for motion analysis. For example, as will be described below with reference to FIGS. 6-16, acquisition of two sets of tagged images in orthogonal planes will permit definition of the full strain tensor throughout the imaged myocardium as a function of cardiac cycle phase, thereby allowing regional myocardial function to be quantitatively assessed.

As will be apparent from the following description, the two-dimensional tagging of the heart wall with a grid of planes of altered magnetization offers the possibility of using finite element analysis to quantify regional heart wall motion. This is accomplished in accordance with the present invention by using the intersections of the tagging stripes as a set of fiducial marks within the heart wall. The motion of these intersections through the heart cycle tracks the motion of the underlying tissue. This grid of reference points provides the basis for a triangular tiling of the heart wall between them. Then, by approximating the regional heart wall motion between each such triangular finite element as being locally homogeneous, the motion of the triangular boundary will completely define the motion within it. In particular, the motion of the region may be decomposed into its equivalent rigid body (translation and rotation) and deformation (strain) components as described by Meier et al. in the aforementioned articles entitled "Kinematics of the Beating Heart", IEEE Transactions on Biomedical Engineering_f Vol. 27, 1980, pp. 319-329, and "Contractile Function in Canine Right Ventricle", American Journal of Physiology. Vol. 239, 1980, pp. H794-H804, the contents of which are incorporated herein by reference as if set forth in their entirety herein. The strain tensor then can be decomposed, in turn, into its principal strains or eigenvectors. Such a decomposition of the motion has the useful property that it is independent of the particular orientation of the tagging grid from which it is derived.

Meier et al. describe that the motion of tissue between groups of markers (or, as described herein, the intersections of the tagging stripes) can be analyzed using the concepts of continuum mechanics. For example, for non-collinear triplets of markers or intersections in accordance with the invention, the assumption of locally homogeneous motion in the triangle between them permits complete characterization of that motion from the time dependence of the coordinates of the markers. In other words, a hypothetical unit circle embedded in the initial state of the tissue will be transformed into an ellipse by the motion, and eigenvalues of the strain give the major and minor axes of that ellipse as shown in FIG. 6. Similarly, these concepts can be extended to three dimensions with groups of four markers or intersections defining a tetrahedron, if the motion is again considered locally homogeneous within the tetrahedron. In other words, as shown in FIG. 7, a hypothetical unit sphere embedded in the initial state would generally be transformed into an ellipsoid by the motion, and the eigenvalues of the transformation would give the axes of the ellipsoid. One method for analyzing and mapping the regional variation in wall motion (hence contractile strength), is to divide the motion of the heart wall into rigid body displacement and regional deformation. Rigid body motion can easily be visualized by a vector representing the displacement of the centroid of a region of the left ventricular wall. Deformation cannot be visualized as easily since it has many components. Thus, one method of displaying information about the regional deformation in the heart wall is through tne calculation of the eigenvalues and eigenvectors of deformation. This method is based on the principles of continuum mechanics and was first used to describe heart wall motion by Meier et al., the assumptions and methods of which will be summarized here.

Even though the deformation of the entire heart varies temporally throughout the heart cycle, the deformation of a small region of the heart can be represented as a linear transformation, provided that the vectors that it operates on are "small enough". In two dimensions, this linear mapping is obtained by calculating the change in lengths of both components of two vectors in the same region, where a vector is defined as:

X = / X1 \, which is mapped into a vector X' = / X1'\ \ X2 / \ X2'/ by a tensor:

F = / F11 F12 \

\ F21 F22 /, where F is actually calculated from a tensor of column vectors, A and B, in the undeformed and deformed state, respectively, where:

A = / X1 X2 \

\ Y1 Y2 /

B = / X1' X2" \

\ Y1' Y2' /. Hence, F = Ainv * B.

By the Polar Decomposition Theorem, F can be separated into two operators, R and U, where R is an orthogonal rotation matrix and U is a positive symmetric stretch tensor which has real positive eigenvalues which are the principal strains of the transformation. The angle a of R is the rigid body rotation. Thus, R can be written as:

R = / cos (α) sin (α) \

\-sin (α) cos (α) /, where a is calculated such that the off-diagonal elements of U are equal. The resulting formula for α is: α = ( F21 - F12) / (F11 + F22).

Equation (19)

The resulting positive symmetric tensor U can be decomposed into its eigenvalues and eigenvectors by setting the determinant of U-I(λ) = 0. This equation for two dimensions results in a quadratic equation for λ, which can be solved leaving two eigenvalues for U. The eigenvectors are obtained by solving the following equation for the components of the eigenvector V1 and V2:

(U) (V) = (A) (V).

Equation (20)

Other types of displays may also give insight into the regional variation by providing scalar representations of the strain vectors, namely the strain invariants. The first invariant, calculated directly from U, equals the sum of the diagonal components of U, i.e., I1 = U11 + U22 in two dimensions. The second invariant equals to the determinant of the elements of U, i.e.: 12 = {(U11)(U22) - (U12)(U21)} (in 2 dimensions). In practice, the vectors connect two intersections of SPAMM stripes so that they range from 5 to 10mm in length.

More details of this analysis technique may be found by referring to the aforementioned articles by Meier et al., as incorporated by reference above.

In order to apply the above-mentioned analysis techniques to "tagged" MR images of the heart in accordance with the invention, an initial approximation for a two- dimensional motion analysis is made that the intersections of stripes in the image of the heart define unique tagged points in the myocardium. In a preferred embodiment, portions of a suitable interactive computer program such as that available on microfiche as Appendix B (which is a listing of an embodiment of an analysis program in accordance with the invention which is written in the programming language "C" for use on a Sun workstation using Sunview display tools) is used for determining the locations of the tagging stripe intersections within the heart wall. These intersections are identified and recorded for a series of images at a given location at different phases of the cardiac cycle. The sequential positions at later phases corresponding to each such identified point in the initial image of the series are then linked. Linked lists of positions can then be used to provide the basis for a triangular tiling of the heart wall by choosing appropriate groups of three neighboring points. Generally, only triangles completely contained within the heart wall can be unambiguously used. The motion of the regions delineated by the resulting sets of serial triangular configurations can then be analyzed as described above for the motion analysis with respect to the serial locations of triplets of implanted metallic markers. This analysis yields a set of variables characterizing the rigid body motion and deformation of each triangular region defined for any pair of cardiac images at different phases. Preferably, for any pair of cardiac images at different phases the initial reference image is defined in such a pair so as to be of a well-defined state such as end diastole. Appendix B sets forth a preferred embodiment of a computer program which performs the above functions and is available on microfiche upon request.

Having measured the location of many such points and calculated the variables characterizing the motion of the corresponding triangular regions they define in accordance with the techniques taught by Meier et al. and set forth above, the resulting large data set must be presented in a readily analyzable form so that the diagnostician can properly analyze the imaged data. This is most readily done visually using computer graphics to either overlay suitable points, lines or shaded or colored areas on the original heart images from which they were derived, or simply to display the derived graphics by themselves as "functional" images of the regional heart wall motion. Preferably, the contours of the inner and outer surfaces of the heart wall are also displayed for reference, particularly when displaying the functional images alone. If desired, the actual computed values of the motion variables for any region in the functional image can still be interactively displayed. Many different variables can be chosen for such displays, but examples of point or line overlays would include the locations of the selected points, the sequential positions of these points (shown as displacement vectors), and the configurations of the triangles defined by the points. The triangles can be shaded or colored according to selected scalar motion variables, such as the magnitudes of the displacement of the centroids, the eigenvalues, or the angle of the eigenvectors. Higher order variables such as the strain tensor can be represented with appropriate symbols.

A generalized flow chart of the computer program of Appendix B running on host computer 212 or a peripheral computer so as to perform the above-identified functions of the preferred embodiment is shown in FIG. 8. As shown, host computer 212 or the peripheral computer reads MR images from the clinical scanner's archive magnetic tapes (storage 220) at step 800. The images are then sorted by location and cardiac cycle phase and filed for subsequent display and analysis at step 802. The images are then flexibly displayed (e.g., with zooming) on display monitor 222 at step 804 so as to aid in picking the locations of the points of intersection of the tagging stripes in the heart wall at step 806. Once the desired stripe intersection points are selected at step 806, triangles are selected at step 808 in accordance with the techniques described above. The necessary steps for motion analysis in accordance with the techniques described by Meier et al. and set forth above, wherein eigenvalues and eigenvectors are calculated, are then performed at step 810. The resulting motion is then displayed at step 812 using a graphics package or some other display device whereby the necessary information is relayed to the diagnostician. Of course, any display technique which relays the desired information may be used by one skilled in the art. For example, the display may compare the triangular files just before imaging to the deformed triangular tiles immediately after imaging.

The analysis in accordance with the present invention can be extended to three-dimensional motion in a variety of different ways. For example, two sets of tagged images can be acquired in different orientations, and the resulting sets of two-dimensional components of the motion may be combined to find the three-dimensional motion of the underlying overlapped volume. In other words, information on the through-plane motion into the MR images is incorporated using an orthogonal set of tagged images. Alternatively, one set of two-dimensional tagged images can be acquired to study in-plane motion, with the addition of phase shift sensitization to through-plane motion. Thus, one set of images can contain three-dimensional motion information, with in-plane motion stored as stripe displacements, thereby integrating the motion between the time of tagging and imaging, and through-plane motion stored as phase shifts proportional to velocity at the time of imaging. With a series of images made at different delays after tagging, the velocity can be integrated or the displacements differentiated to give the three-dimensional motions on a common scale. With both approaches to three-dimensional motion, the images are recorded at fixed planes in space, yielding an Eulerian description of the motion. Such a description is convenient for transforming to an equivalent Lagrangian description of the motion in material coordinates relative to an initial configuration at the time of tagging.

Thus, the striped intersections in an image in accordance with the invention need not define unique points, thereby allowing the possibility of a component of motion through the image plane. A triangle of intersection points in one image plane and a suitable intersection point in an adjacent image plane may also define a tetrahedron. Displacement of the vertices of the initial tetrahedron can be estimated for subsequent cardiac phases from suitable interpolation of the MR motion images. This data can then be analyzed as for the corresponding data from tetrahedral quadruplets of implanted metallic members in accordance with the techniques used by Fenton et al. in an article entitled "Transmural Myocardial Deformation in the Canine Left Ventricular Wall", American Journal of Physiology. Vol. 235, 1978, pp. H523-H530, and by Waldman et al. in an article entitled "Transmural Myocardial Deformation in the Canine Left Ventricle", Circulation Research. Vol. 57, 1985, pp. 152-163, for example.

Because of the higher dimensionality of three-dimensional data and the resulting variables characterizing the motion, it is more difficult to display such data in a readily analyzable form. Possibilities for such three-dimensional display include stereo or interactively rotatable displays such as wire frame displays of the tetrahedrons, shaded surface displays of the heart reflecting the underlying motion variables, and interactive "cutaway" displays of the heart wall. Of course, other display techniques in accordance with the invention will become apparent to those skilled in the art. Representative images of the intermediate steps and the resulting final displays are shown in FIGURES 9-16 for the preferred embodiment.

For example, FIG. 9 illustrates two-dimensional tagged MR short axis images of the heart of a human immediately after tagging at the end of diastole (left) and at the same level in late systole (right). FIG. 10 illustrates the points specified for stripe intersections for the images in FIG. 9, while FIG. 11 illustrates the displacement of specified points between the times of the taking of the two images.

FIG. 12 illustrates the corresponding triangulation for the points shown in FIG. 9, while FIG. 13 illustrates the magnitude of mean displacement of triangles between the imaging times of FIG. 9, where different grey scale magnitudes represent differing displacement magnitudes for different triangular regions.

FIG. 14 illustrates eigenvectors of the transformation of triangles between the imaging times of FIG. 9, represented as a symbol with an initial unit circle and superimposed major and minor axes of a corresponding subsequent ellipse into which it would be transformed. FIG. 15 illustrates the angle of the principal eigenvectors for the triangles of FIG. 14, displayed as grey levels.

FIG. 16 illustrates a combined display of point displacements (vectors) and the magnitude of mean displacement of triangles derived from long axis views of the same subject and at the same times as in FIG. 9.

The analysis technique of the invention allows the diagnostician to interactively pick corresponding locations of tag intersections on consecutive images. The resulting set of points can be interactively used to specify a triangular tiling of the wall from which an eigenvalue analysis of the regional deformation may be carried out. Displays of the derived motions that may be found useful by those skilled in the art include graphic overlays on the initial images of an "optical flow" display of the serial displacements of each tagged point, a display of the initial and deformed grids defined by the tiling, and a set of crossed line segments superimposed on each triangle, whose lengths and orientations correspond to the regional principal strains. Of course, other display techniques may be used by those of ordinary skill in the art.

Thus, the approach to analysis of tagged MR images of heart wall motion in accordance with the invention can provide quantitative measures of regional heart wall motion. The values calculated for the variables characterizing the action do not in principle depend on a specific orientation of the tagging grid used. Instead, these values can be used to generate functional images that directly display the regional distribution of wall motion. Moreover, the technique of the invention is also useful for studying the flows of blood or cerebrospinal fluid (CSF) and for distinguishing slow flowing blood from thrombus. In addition, by using the techniques of the invention, an infarct can be localized and identified and then displayed on a display monitor for analysis by a diagnostician. Furthermore, not only may the technique of the invention be used to separately image motion within the heart wall, but also the invention may be generally used to image elastic strain deformation in all soft tissue structures.

However, a limitation of the motion analysis from tagged images in accordance with the invention is that the spatial sampling is limited by the stripe spacing, which is, in turn, limited by the imaging resolution. The duration of stripe persistence is thus limited both by longitudinal relaxation of the tissue and blurring caused by respiratory motion. The spacing of the stripes can be easily adjusted by changing the strength and duration of the gradient pulse in the pre-imaging sequence in accordance with the techniques described above. By so adjusting the gradient pulse, stripes as fine as tens of microns apart have been created, thereby enabling very small motions involved in diffusion and perfusion to be monitored. Thus, the SPAMM imaging technique in accordance with the invention provides a simple and effective way to directly demonstrate cardiac motion. It is also useful for studying the flows of blood or CSF and for distinguishing slow flowing blood from thrombus. Adaptation of the methods described in accordance with this invention permit a wide range of other applications, including demonstrating the distribution of magnetic field and RF field inhomogeneities and gradient non- linearities, calculating gradients and demonstrating chemical shift differences. Moreover, the MR imaging technique of the invention is much more flexible than any previously contemplated for motion analysis.

Although an exemplary embodiment of this invention has been described in detail above, those skilled in the art will readily appreciate that many additional modifications are possible in the exemplary embodiment without materially departing from the novel teachings and advantages of this invention. For example, the radio frequency pulses need not be applied along the same axis. They can instead be incremented in value to create a net phase shift. Also, since the pulses herein described are designed to the specifications of the observer, minor variations in the amplitude and phase of the optimal pulses are believed to be within the scope of the invention as herein described and claimed. In addition, the imaging described herein can be applied to a patient's blood or the patient's cerebrospinal fluid as well as the patient's heart wall to determine motion in accordance with the techniques of the invention. Furthermore, since the most time-consuming part of the analysis process of the invention is the interactive picking of locations of stripe intersections, image processing approaches may be used to speed up or automate this step so as to greatly facilitate the motion analysis. Accordingly, all such modifications are intended to be included within the scope of this invention as defined in the following claims. APPENDIX A

0001 IMPLICIT REAL*8(A-H,O-Z)

0002 COMMON

PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID

0003 DIMENSION IEXT(512),AD(512),ALPHA(512),X(512),Y(512) 0004 DIMENSION H(512)

0005 DIMENSION DES(8196),GRID(8196),WT(8196)

0006 DIMENSION EDGE(20),FX(10),WTX(10),DEVIATE(10)

0007 DIMENSION XINT(8196)

0008

0009 PI2=6.283185307179586

0010 PI=3.141592653589793

0011 NFMAX = 256 ¦ MAXIMUM LENGTH OF EXPONENTIAL

SERIES

0012 NNBBAANNDDSS=2 ¦ HAS 2 BANDS

0013 JB=2*NBANDS

0014 LGRID=6 ¦ GRID DENSITY

0015 C

0016 C

0017 C

0018 C PROGRAM INPUT SECTION

0019 C HERE WILL SPECIFY INPUT

0020

0021 PRINT *,' PULSE SYNTHESIS PROGRAM '

0022 PRINT *,' ************************************' 0023 PRINT *,' PROGRAM WILL SYNTHESIZE A HARD PULSE

SEQUENCE YIELDING

0024 PRINT *,' ANY DESIRED FREQUENCY SPECIFICATIONS 0025 PRINT *,' NUMBER OF PULSES ? '

0026 READ *,NFCNSC

0027 NFILT=2*NFCNS-1

0028 IF (NFILT .GT. NFMAX .OR. NFILT .LT. 3) CALL ERROR 0029

0030

0031 C THERE ARE TWO BANDS

0032 C THE FIRST BAND RUNS FROM THE MIDDLE OF THE STRIPE TO

THE END OF T

0033 C STRIPE

0034 C THE SECOND BAND (WHICH IS THE PART BETWEEN THE

STRIPES) RUNS FROM

0035 THE END OF THE TRANSITION ZONE TO THE

MIDDLE OF THE

0036 C INTERSTRIPE AREA

0037 C THE TRANSITION ZONE IS THE PART BETWEEN THE TWO

BANDS

0038 C THE WEIGHTS DETERMINE THE RATIO OF THE RIPPLE

AMPLITUDES IN THE

0039 C BANDS - MEANINGFUL ONLY IF THE STRIPE IS

WIDER THAN ONE

0040 C RIPPLE (A FUNCTION OF THE NUMBER OF PULSES

AND THE WIDTH

0041 C OF THE STRIPE)

0042 PRINT *,' NOW ENTER DESIRED VALUES. '

0043 PRINT *,' BAND EDGES UNITS ARE FRACTION OF

INTERSTRIPE DISTANCE' 0044 EDGE(1)=0

0045 BAND # 1 STARTS AT 0 AND ENDS AT

0046 READ *,EDGE(2)

0047 PRINT *,' AND SHOULD HAVE VALUE ? '

0048 READ *,FX(1)

0049 PRINT *,' AND SHOULD HAVE WEIGHT ? '

0050 READ *,WTX(1)

0051 PRINT *,' BAND 2 ENDS AT .5 BUT STARTS AT??'

0052 READ *,EDGE(3)

0053 FX(2)=1

0054 PRINT *,' AND SHOULD HAVE WEIGHT ? '

0055 READ *,WTX(2)

0056 EDGE (4) =.5

0057 C

0058 C ************************************

0059 C WILL NOW USE PARK - MCCLELL*AN ALGORITHM TO GENERATE

FOURIER SERIE

0060 C MAGNETIZATION

0061 C

0062 GRID (1)=EDGE(1)

0063 DELF=LGRID*NFCNS

0064 DELF=0.5D0/DELF

0065 J=1

0066 L=1

0067 LBAND=1

0068 140 FUP=EDGE(L+1)

0069 145 TEMP=GRID(J)

0070 C

0071 C CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE

0072 C WEIGHT FUNCTION ON THE GRID

0073 C

0074 DE5 (J) =EFF (TEMP, FX,WTX, LBAND)

0075 WT (J) =WATE (TEMP, FX,WTX, LBAND)

0076 J=J+1

0077 GRID (J) =TEMP+DELF

0078 IF(GRID(J). GT. FUP)GO TO 150

0079 GO TO 145

0080 150 GRID(J-1)=FUP

0081 DES (J-1) =EFF (PUP,FX,WTX,LBAND)

0082 WT (J-1) =WATE (FUP, FX,WTX, LBAND)

0083 LBAND=LBAND+1

0084 L=L+2

0085 IF (LBAND .GT. NBANDS) GO TO 160

0086 GRID(J)=EDGE(L)

0087 GO TO 140

0088 160 NGRID=J-1

0089

0090

0091 C

0092 C SET UP A NEW APPROXIMATION PROBLEM WHICH IS

0093 C EQUIVALENT TO THE ORIGINAL PROBLEM

0094 C

0095 C INITIAL GUESS FOR THE EXTERMAL FREQUENCIES—EQUALLY

0096 C SPACED ALONG THE GRID

0097 C

0098 200 TEMP=FLOAT (NGRID-1)/FLOAO (NFCNS) 0099 DO 210 J=1,NFCNS

0100 210 IEXT(J)=(J-1)*TEMP+1

0101 IEXT (NFCNS+1)=NGRID

0102 NH1=NFCNS-1

0103 NZ=NFCNS+1

0104

0105 C CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE C

APPROXIMATION PROBL

0106 C

0107 CALL REMEZ (EDGE,NBANDS)

0108 C

0109 DO J=1,2048

0110 XINT(J)=0.D0

0111 ENDDO

0112 dev=dabs (dev)

0113 aconl=1/ (1+dev)

0114 C PRINT *,ACON1

0115 DO J=1,NFCNS

0116 XINT(J)=ALPHA(J)*acon1

0117 ENDDO

0118 PRINT *,XINT(1)

0119 C WILL NOW NORMALIZE DO FFT OF XINT

0120 CALL FFT(xint)

0121 amax=0.0

0122 do j=1,4096

0123 x1=dabs (XINT(j)/2.d0)

0124 if(x1 .gt. amax) amax=x1

0125 enddo

0126 C CALCULATE THE IMPULSE RESPONSE

0127 C

0128 acon-aconl/amax

0129 PRINT *,' AMAX=',AMAX

0130

0131 OPEN(UNIT=21,NAME= 'poly. DAT' ,TYPE='NEW' )

0132 C HERE ARE STORED COEEFICIENTS OF FOURIER SERIES

(COMPLEX EXPONENTI

0133 C FORM - NEED TO CONVERT FROM COSINE FORM 0134 write (21,*) 2*nfcns-1

0135 do j=1,nfcns-1

0136 h(j)=alpha(nfcns+1-j)*acon/2.do

0137 write(21,*)h(j)

0138 enddo

0139 write(21,*) alpha (1) *acon

0140 do j=1,nfcns-1

0141 write (21,*)alpha(j+1)*acon/2.do

0142 enddo

0143 close (unit=21)

0144

0145 CALL PULS (NBANDS)

0146 STOP

0147 C

0148 END 0001

0002 REAL*8 FUNCTION EFF (TEMP, FX,WTX, LBAND)

0003 C

0004 C FUNCTION TO CALCULATE THE DESIRED

MAGNITUDE RESPONSE

0005 C AS A FUNCTION OF FREQUENCY,

0006 C

0007 IMPLICIT REAL*8(A-H,O-Z)

0008 DIMENSION FX(5),WTX(5)

0009 EFF=FX (LBAND)

0010 RETURN

0011 END

0001

0002

0003 REAL*8 FUNCTION WATE (TEMP, FX,WTX, LBAND)

0004 C FUNCTION TO CALCULATE THE WEIGHT FUNCTION

0005 C OF FREQUENCY,

0006 C

0007 IMPLICIT REAL*8(A-H,O-Z)

0008 DIMENSION FX(5),WTX (5)

0009 WATE=WTX(LBAND)

0010 RETURN

0011 END

0001

0002 SUBROUTINE ERROR

0003 PRINT 1

0004 1 FORMAT ( '********** ERROR IN INPUT DATA **********') 0005 STOP

00006 END

0001

0002

0003 SUBROUTINE REME2 (EDGE,NBANDS)

0004 C

0005 C THIS SUBROUTINE IMPLEMENTS THE REME2 EXCHANGE

0006 C ALGORITHM

0007 C FOR THE WEIGHTED CHEBYCHEV APPROXIMATION OF A

CONTINUOUS

0008 C FUNCTION WITH A SUM OF COSINES. INPUTS TO THE

SUBROUTINE

0009 C ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS,

THE

0010 C DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION

ON THE

0011 C GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF

THE

0012 C EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE

CHEBYCHEV

0013 C ERROR BY DETERMINING THE BEST LOCATION OF THE

EXTREMAL

0014 C FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN

CALCULATES

0015 C THE COEFFICIENTS OF THE BEST APPROXIMATION.

0016 C

0017 IMPLICIT REAL*8(A-H,0-Z)

0018 COMMON

PI2,AD,DEV,X,Y,GRID, DES ,WT,ALPHA,TEXT,NFCNS ,NGRID

0019 DIMENSION EDGE (20)

0020 DIMENSION IEXT(512),AD(512),ALPHA(512),X(512),Y(512) 0021 DIMENSION DES (8196),GRID(8196),WT(8196)

0022 DIMENSION A(256),P(256),Q(256)

0023

0024 C THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF

25

0025 C

0026 ITRMAX=25

0027 DEVL=-1.0

0028 NZ=NFCN5+1

0029 NZZ=NFCNS+2

0030 NITER-0

0031 100 CONTINUE

0032 IEXT (NZZ) =NGRID+1

0033 NITER=NITER+1

0034 IF(NITER .GT. ITRMAX) GO TO 400

0035 DO 110 J=1,NZ

0036 DTEMP=GRID (IEXT (J))

0037 DTEMP=DCOS (DTEMP*PI2)

0038 110 X(J)=DTEMP

0039 JET=(NFCNS-1)/15+1

0040 DO 120 J=1,NZ

0041 120 AD(J)=D(J,NZ,JET)

0042 DNUM=0.0

0043 DDEN=0.0

0044 K*=l

0045 DO 130 J=1,NZ

0046 L=IEXT(J) 0047 DTEMP=AD (J) *DES (L)

0048 DNUM=DNUM+DTEMP

0049 DTEMP=K*AD(J)/WT(L)

0050 DDEN=DDEN+DTEMP

0051 130 K=-K

0052 DEV=DNUM/DDEN

0053 NU=1

0054 IF(DEV .GT. 0.0) NU=-1

0055 DEV=-NU*DEV

0056 K=NU

0057 DO 140 J=1,NZ

0058 L=IEXT(J)

0059 DTEMP=K*DEV/WT (L)

0060 Y(J) =DES (L) +DTEMP

0061 140 K=-K

0062 IF(DEV .GE. DEVL) GO TO 150

0063 CALL OUCH

0064 GO TO 400

0065 150 DEVL=DEV

0066 JCHNGE=0

0067 K1=IEXT(1)

0068 KNZ=IEXT(NZ)

0069 KLOW=O

0070 NUT=-NU

0071 J=1

0072 C

0073 C SEARCHFOR THE EXTREMAL FREQUENCIES OF THE BEST

0074 C APPROXIMATION

0075 C

0076 200 IF(J .EQ. NZZ) YNZ=COMP

0077 IF(J .GE. NZZ) GO TO 300

0078 KUP=IEXT(J+1)

0079 L=IEXT(J)+1

0080 NUT=-NUT

0081 IF(J EQ. 2) Y1=COMP

0082 COMP=DEV

0083 IF(L .GE. KUP) GO TO 220

0084 ERR=GEE(L,NZ)

0085 ERR= (ERR-DES (L) ) *WT (L)

0086 DTEMP=NUT*ERR-COMP

0087 IF(DTEMP .LE. 0.0) GO TO 220

0088 COMP=NUT*ERR

0089 210 L=L+1

0090 IF(L .GE. KUP) GO TO 215

0091 ERR=GEE(L,NZ)

0092 ERR=(ERR-DES (L))*WT(L)

0093 DTEMP=NUT*ERR-COMP

0094 IF(DTEMP .LE. 0.0) GO TO 215

0095 COMP=NUT*ERR

0096 GO TO 210

0097 215 IEXT(J)=L-1

0098 J=J+1

0099 KLOW=L-1

0100 JCHNGE=JCHNGE+1

0101 GO TO 200

0102 220 L=L-1 0103 225 L=L-1

0104 IF(L .LE. KLOW) GO TO 250

0105 ERR=GEE(L,NZ)

0106 ERR=(ERR-DES (L))*WT(L)

0107 DTEMP=NUT*ERR-COMP

0108 IF(DTEMP .GT. 0.0) GO TO 230

0109 IF(JCHNGE .LE. 0) GO TO 225

0110 GO TO 260

0111 230 COMP=NUT*ERR

0112 235 L=L-1

0113 IF(L .LE. KLOW) GO TO 240

0114 ERR=GEE(L,NZ)

0115 ERR=(ERR-DES (L))*WT(L)

0116 DTEMP=NUT*ERR-COMP

0117 IF(DTEMP .LE. 0.0) GO TO 240

0118 COMP=NUT*ERR

0119 GO TO 235

0120 240 KLOW=IEXT(J)

0121 IEXT(J)-L+1

0122 J=J+1

0123 JCHNGE=JCHNGE+1

0124 GO TO 200

0125 250 L=IEXT(J)+1

0126 IF(JCHNGE .GT. 0) GO TO 215

0127 255 L=L+1

0128 IF(L .GE. KUP) GO TO 260

0129 ERR=GEE(L,NZ)

0130 ERR= (ERR-DES (L))*WT(L)

0131 DTEMP=NUT*ERR-COMP

0132 IF(DTEMP .LE. 0.0) GO TO 255

0133 COMP=NUT*ERR

0134 GO TO 210

0135 260 KLOW=IEXT(J)

0136 J=J+1

0137 GO TO 200

0138 300 IF (J .GT. NZZ) GO TO 320

0139 IF (Kl .GT. IEXT(l)) K1=IEXT(1)

0140 IF(KNZ .LT. IEXT(NZ)) KNZ=IEXT(NZ)

0141 NUT1=NUT

0142 NUT=-NU

0143 L=O

0144 KUP=K1

0145 COMP=YNZ* (1.00001)

0146

0147 LUCK=1

0148 310 L=L+1

0149 IF(L .GE. KUP) GO TO 315

0150 ERR=GEE(L,NZ)

0151 ERR= (ERR-DES (L) ) *WT(L)

0152 DTEMP=NUT*ERR-COMP

0153 IF(DTEMP .LE. 0.0) GO TO 310

0154 COMP=NUT*ERR

0155 J=NZZ

0156 GO TO 210

0157 315 LUCK=6

0158 GO TO 325 0159 320 IF(LUCK .GT. 9) GO TO 350

0160 IF(COMP .GT. Yl) Y1=COMP

0161 K1=IEXT(NZZ)

0162 325 L=NGRID+1

0163 KLOW=KNZ

0164 NUT=-NUT1

0165 C0MP=Y1* (1.00001)

0166 330 L=L-1

0167 IF(L .LE. KLOW) GO TO 340

0168 ERR=GEE(L,NZ)

0169 ERR= (ERR-DES (L))*WT(L)

0170 DTEMP=NUT*ERR-COMP

0171 IF(DTEMP .LE. 0.0) GO TO 330

0172 J=NZZ

0173 COMP=NUT*ERR

0174 LUCK=LUCK+10

0175 GO TO 235

0176 340 IF (LUCK .EQ. 6) GO TO 370

0177 DO 345 J=1,NFCNS

0178 345 IEXT (NZZ-J) =IEXT(NZ-J)

0179 IEXT(1)=K1

0180 GO TO 200

0181 350 KN=IEXT(NZZ)

0182 DO 360 J-1,NFCNS

0183 360 IEXT (J) =IEXT (J+1)

0184 IEXT(NZ)=KN

0185 GO TO 100

0186 370 IF(JCHNGE .GT. 0) GO TO 100

0187

0188 C CALCULATION OF THE COEFFICIENTS OF THE BEST

APPROXIATION

0189 C USING THE INVERSE DISCRETE FOURIER TRANSFORM

0190 C

0191 400 CONTINUE

0192 NM1=NFCNS-1

0193 FSH=1.0E-06

0194 GTEMP=GRID(1)

0195 X(NZZ)=-2.0

0196 CN=2*NFCNS-1

0197 DELF=1.0/CN

0198 L=1

0199 KKK=0

0200 IF(EDGE(1) .EQ. O.O .AND. EDGE (2*NBANDS) .EQ. 0.5)

KKK=1

0201 IF(NFCNS .LE. 3) KKK=1

0202 IF(KKK .EQ.1) GO TO 405

0203 DTEMP=DCOS (II2*GRIDS (1) )

0204 DNUM=DCOS (PI2*GRID(NGRID))

0205 AA=2. O/DTEMP-DNUM)

0206 BB=-(DTEMP+DNUM)/ (DTEMP-DNUM)

0207 405 CONTINUE

0208 DO 430 J=1,NFCNS

0209 FT=(J-1)*DELF

0210 XT=DCOS(PI2*FT)

0211 IF(KKK .EQ. 1) GO TO 410

0212 XT=(XT-BB)/AA 0213 FT=DACOS (XT) /PI2

0214 410 XE=X(L)

0215 IF (XT .GT. XE) GO TO 420

0216 IF((XE-XT) .T. FSH) GO TO 415

0217 L=L+1

0218 GO TO 410

0219 415 A(J)=Y(L)

0220 GO TO 425

0221 420 IF((XT-XE) .LT FSH) GO TO 415

0222 GRID(1)=FT

0223 A(J)=GEE(1,NZ)

0224 425 CONTINUE

0225 IF(L .GT. 1) L=L-1

0226 430 CONTINUE

0227 GRID(1)= GTEMP

0228 DDEN=PI2/CN

0229 DO 510 J=1,NFCNS

0230 DTEMP=0.0

0231 DNUM=(J-1)*DDEN

0232 IF (NMl .LT. 1) GO TO 505

0233 DO 500 K=1,NM1

0234 500 DTEMP=DTEMP+A(K+1) *DC05 (DNUM*K)

0235 505 DTEMP=2.0*DTEMP+A(1)

0236 510 ALPHA(J)=DTEMP

0237 DO 550 J+2,NFCNS

0238 550 ALPHA(J) =*ALPHA(J)/CN

0239 ALPHA(1) =ALPHA(1)/CN

0240 IF(KKK .EQ. 1) GO TO 545

0241 P(1) =2.0*ALPHA(NFCNS) *BB+ALPHA(NMl)

0242 P (2) =2.0*AA*ALPHA(NFCNS)

0243 Q (1) =ALPHA(NFCNS-2) -ALPHA(NFCNS)

0244 DO 540 J=2,NM1

0245 IF(J .LT. NMl) GO TO 515

0246 AA=0.5*AA

0247 BB=0.5*BB

0248 515 CONTINUE

0249 P(J+1)=0.0

0250 DO 520 K=1,J

0251 A(K)=P(K)

0252 520 P(K)=2.0*BB*A(K)

0253 P(2)=P(2)+A(1)*2.0*AA

0254 JM1=J-1

0255 DO 525 K=1,JM1

0256 525 P (K) =P(K) +Q (K) +AA*A(K+1)

0257 JP1=J+1

0258 DO 530 K=3,JP1

0259 530 P (K) =P (K) +AA*A(K-l)

0260 IF(J .EQ. NMl) GO TO 540

0261 DO 535 K=1,J

0262 535 Q(K)=-A(K)

0263 Q(1)=Q(1)+ALPHA(NFCNS-1-J)

0264 540 CONTINUE

0265 DO 543 J=1,NFCNS

0266 543 ALPHA(J)=P(J)

0267 545 CONTINUE

0268 IF (NFCNS .GT. 3) RETURN 0269 ALPHA(NFCNS+1) =0.0

0270 ALPHA(NFCNS+2) =0.0

0271 RETURN

0272 END

0003 REAL*8 FUNCTION D(K,N,M)

0004 C

0005 C FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION 0006 C COEFFICIENTS FOR USE IN THE FUNCTION GEE

0007 C

0008 INPLICIT REAL*8 (A-H,O-Z)

0009 COMMON PI2 ,AD, DEV,X,Y, GRID, DES,WT,ALPHA,

IEXT,NFCNS ,NGRID

0010 DIMENSION IEXT(512),AD(512),ALPHA(512),X(512),Y(512) 0011 DIMENSION DES (8196),GRID(8196),WT(8196)

0012

0013 D-1.0

0014 Q=X(K)

0015 DO 3 L=1,M

0016 DO 2 J=L,N,M

0017 IF (J-K) 1,2,1

0018 1 D=2.0*D*(Q-X(J))

0019 2 CONTINUE

0020 3 CONTINUE

0021 D=1.0/D

0022 RETURN

0023 END

0002 REAL*8 FUNCTION GEE(K,N)

0003 C

0004 C FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING

THE

0005 C LAGRANGE INTERPOLUTION FORMULA IN THE BARYCENTRIC

FORM

0006 C

0007 IMPLICIT REAL*8(A-H,0-Z)

0008 COMMON PI2 ,AD, DEV,X,Y,GRID, DES ,WT,ALPHA, IEXT,

NFCNS ,NGRID

0009 DIMENSION IEXT(512),AD(512),ALPHA(512),X(512),Y(512) 0010 DIMENSION DES (8196),GRID(8196),WT(8196)

0011

0012 P=0.0

0013 XF=GRID(K)

0014 XF=DCOS(PI2*XF)

0015 D=0.0

0016 DO 1 J=1,N

0017 C=XF-X(J)

0018 C=AD(J)/C

0019 D=D+C

0020 1 P=P+C*Y(J)

0021 GEE=P/D

0022 RETURN

0023 END

0003 SUBROUTINE OUCH

0004 PRINT1

0005 1 FORMAT (' ************ FAILURE TO CONVERGE

**********'/

0006 l'OPROBABLE CAUSE IS MACHINE ROUNDING ERROR'/

0007 2'OTHE IMPULSE RESPONSE MAY BE CORRECT'/

0008 3'OCHECK WITH A FREQUENT RESPONSE') 0009 RETURN

0010 END

0001

0002 SUBROUTINE PULS (NPOLY)

0003 C PROGRAM WILL TAKE INPUT FROM PARK MCCLELLAN

ALGORITHM AND USE IT

0004 C DO CEPSTRAL DECONVOLUTION

0005 C IT WILL THEN FIGURE OUT FOUR PULSE SEQUENCES

0006 C AND GIVE MINIMAN MAXIMUM DIFFERENCE

0007 IMPLICIT REAL*8(A-H,O-Z)

0008 COMPLEX*16 COEFF(256)

0009 REAL*8 POLY (256),POLY1(256),POLY2(256)

0010 REAL*8 CEP1(256),CEP2(256),CEP3(256),CEP4(256) 0011 REAL*8 PULS1(256),PULS2(256),PULS3(256),PULS4(256) 0012 REAL*8 PULS5(256),PULS6(256),PULS7(256),PULS8(256) 0013 REAL*8 XINT1(256),PULSE(256)

0014 OPEN(UNIT=21,NAME='POLY.DAT',TYPE='OLD')

0015 READ (21,*) NFILT

0016 DO J=1,NFILT

0017 READ (21,*) POLY (J)

0018 POLY (J) =POLY (J)*.999999

0019 POLY1(J)=POLY(J)

0020 POLY2 (J) =-POLY(J)

0021 CEP1(J)=0.0

0022 CEP2(J)=0.0

0023 CEP3(J)=0.0

0024 CEP4(J)=0.0

0025 PULS1(J)=0.0

0026 PULS2(J)=0.0

0027 PULS3(J)=0.0

0028 PULS4(J)=0.0

0029 PULS5(J)=0.0

0030 PULS6(J)=0.0

0031 PULS7(J)=0.0

0032 PULS8(J)=0.0

0033 ENDDO

0034 CLOSE (UNIT=21)

0035 NFCNS=NFILT/2+1

0036 POLY1 (NFCNS) =POLY1 (NFCNS) +1

0037 POLY2 (NFCNS) =POLY2 (NFCNS)+1

0038

0039 DO J=1,NFILT

0040 COEFF (J) =POLY1(J)/2

0041 ENDDO

0042 CALL CEPSTRUM(COEFF,CEP1,CEP2,NFILT)

0043 do j=1,nfcns

0044 cep2(j)=-cepl(j)

0045 enddo

0046 DO J=1,NFILT

0047 COEFF(J) =POLY2 (J)/2

0048 ENDDO

0049 CALL CEPSTRUM(COEFF, CEP3 , CEP4 ,NFILT)

0050 do j=1,nfcns

0051 cep4(j)=-cep3(j)

0052 enddo

0053

0054 CALL PULSR1 (CEP1, CEP3 ,NFCNS , PULS1, DIFMAX1)

0055 DIFMAX=DIFMAX1 0056 DO J=1,NFCNS

0057 PULSE(J)=PULSE1(J)

0058 ENDDO

0059 CALL PULSR1 (CEP2 , CEP3 ,NFCNS , PULS2 , DIFMAX2)

0060

0061 IF(DIFMAX2 .1T. DIFMAX) THEN

0062 DIFMAX=DIFMAX2

0063 DO J=1,NFCNS

0064 PULSE (J) =PULS2 (J)

0065 ENDDO

0066

0067 ENDIF

0068 DO J=1,NFCNS

0069 CEP1(J)=-CEP1(J)

0070 CEP2(J)=-CEP2(J)

0071 ENDDO

0072 CALL PULSR1 (CEP1,CEP3,NFCNS,PULS3,DIFMAX3)

0073

0074 IF(DIFMAX3 .1T. DIFMAX) THEN

0075 DIFMAX=DIFMAX3

0076 DO J=1,NFCNS

0077 PULSE (J) =PULS3 (J)

0078 ENDDO

0079

0080 ENDIF

0081 CALL PULSRl (CEP2 EP4,NFCNS,PULS4,DIFMAX4)

0082

0083 IF(DIFMAX4 .LT. DIFMAX) THEN

0084 DIFMAX=DIFMAX4

0085 DO J=1,NFCNS

0086 PULSE (J)=PULS4(J)

0087 ENDDO

0088 ENDIF

0089 C WRITE RESULTS TO FILE

0090 OPEN(UNIT=21, TYPE= 'NEW' ,NAME= ' PULSE. OUT ' )

0091 WRITE (21,*) NFCNS

0092 DO J=1,NFCNS

0093 WRITE (21,*) PULSE (J)

0094 ENDDO

0095 CLOSE (UNIT=21)

0096 RETURN

0097 END

0001

0002 SUBROUTINE CEPSTRUM(COEFF, CEPMAX, CEPMIN,N)

0003 IMPLICIT COMPLEX*16 (A-H, O-Z)

0004 COMPLEX*16 COEFF (1), CINT (4096), CINT1 (4096)

0005 REAL*8 XINT1,XINT2, CEPMAX(N), CEPMIN(N)

0006 M=N/2+1

0007 DO J=1,4096

0008 CINT(J)=0 ¦ INITIALIZE

0009 ENDDO

0010 DO J=1,N

0011 CINT(J)=COEFF(J)

0012 ENDDO

0013 CALL FFT2C(CINT)

0014 DO J=1,4096

0015 XINT1=CDABS (CINT (J))

0016 IF(XINTl .LT. l.D-20) XINT1=1.D-20

0017 CINT (J)=DCMPLX(DLOG(XINTl)/2.) ¦ computes

logofsqrt

0018 ENDDO

0019 c now do fft

0020 CALL FFT(CINT)

0021

0022 do j=1,2047

0023 CINT1 (J+1) =DCONJG(CINT(J+1)/2048)

0024 cint(j+1)=0

0025 CINT1(4097-J)=O

0026 cint(4097-j)=dconjg(cint(4097-j)/2048

0027 enddo

0028 cint(2049)=0

0029 cint(1)=dconjg(cint(1)/4096)

0030 cint1(2049)=0

0031 cint1(1)=dconjg(cint(1)/4096)

0032

0033 c now have cepstrum; need to compute minimal

phase

0034 call fft (cint)

0035 do j=1,4096

0036 cint(j)=dconjg(cdexp(cint(j) ) )

0037 enddo

0038 call fft(cint)

0039 cepMIN (1) =dREAL(cint(1)/4096)

0040 DO J=2,M

0041 CEPMIN(J)=DREAL(CINT(4098-j))/4096

0042 ENDDO

0043 xint1=dsqrt(dabs(dreal(coeff(1))/(cepmin(1)*cepmin(m

))))

0044 do j=1,m

0045 cepmin(j)=cepmin(j)*xint1

0046 enddo

0047 C WILL NOW COMPUTE MAXIMAL PHASE

0048 call fft(cint1)

0049 do j=1,4096

0050 cint1(j)=dconjg(cdexp(cint1(j)))

0051 enddo

0052 call fft(cint1)

0053 DO J=1,M 0054 CEPMAX(J)=DREAL(CINT1(j))/4096

0055 ENDDO

0056 xintl=dsqrt(dabs(dreal(coeff(1))/(cepmax(1)*cepmax(m

))))

0057 do j=1,m

0058 cepmax(j)=cepmax(j)*xint1

0059 enddo

0060 RETURN

0061 END

0001

0002

0003 SUBROUTINE PULSRl (POL1, POL2 ,N, PULSES , DIFMAX) 0004 C WILL INTAKE THE SPINOR COMPONENTS OF THE

MAGNETIZATION AND OUTPUT

0005 C FLIP ANGLES

0006 C THIS ONE ASSUMES REAL INPUTS

0007 c spinor is (poll,pol2)

0008 c pol1 is a1 + a2s +a2s -2+... as is pol2

0009 IMPLICIT REAL*8(A-H,O-Z)

0010 REAL*8 POLl (1), POL2 (1), PULSES (1)

0011 REAL*8 XINTl (200) ,XINT2(200) ,XINT3(200) ,XINT4(200) 0012 C FIRST INITIALIZE XINTl,XINT2

0013 DO J=1,N

0014 XINT1 (J)=POL1(J)

0015 XINT2(J)=POL2(J)

0016 ENDDO

0017 DIFMAX=0

0018 DIFO=0

0019 DIFN=0

0020 C WILL NOW LOOP OVER

0021 NF1=(N-1)/2

0022 DO J=1,NF1

0023 J1=N+1-2*J

0024 FLIP=DATANZD(XINT2 (1) ,XINTl (1))

0025 IF((FLIP) .GT. 90) FLIP=FLIP-180

0026 IF((FLIP) .LT. -90) FLIP=FLIP+180

0027 Xl=DCOSD(FLIP)

0028 Y1=DSIND(FLIP)

0029 DIF1=2*FLIP

0030 PULSES (N+1-J) =DIF1

0031 DO K=1,J1

0032 XINT1 (K) =X1*XINT1 (K) +Y1*XINT2 (K)

0033 XINT2 (K) =-Y1*XINT1 (K+1) +X1*XINT2 (K+1) 0034 ENDDO

0035 DIF=ABS (DIF1-DIF0)

0036 IF(DIF .GT. DIFMAX) DIFMAX=DIF

0037 DIF0=DIF1

0038 C NOW FROM OTHER END

0039 FLIP=DATAN2D (XINTl (Jl) , -XINT2 (Jl) )

0040 IF((FLIP) .FT. 90) FLIP=FLIP-180

0041 IF (FLIP .LT. -90 FLIP=FLIP+180

0042 Xl=DCOSD(FLIP)

0043 Y1=DSIND(FLIP)

0044 DIF1=2*FLIP

0045 PULSES (J)=DIF1

0046 DO K=1,J1-1

0047 XINT2 (K) =X1*XINT1 (K) -Y1*XINT2 (Jl+l-K) 0048 XINT4 (K) =-Yl*XINTl (Jl=1-K)+X1*XINT2 (K) 0049 ENDDO

0050 DO K=1,J1-1

0051 XINT1 (K)=XINT3(K)

0052 XINT2(K)=XINT4(K)

0053 ENDDO

0054 DIF=ABS (DIF1-DIFON)

0055 IF(DIF .GT. DIFMAX) DIFMAX=DIF 0056 DIFN=DIF1

0057

0058 ENDDO

0059 FLIP=DATAN2D(XINT2 (1), XINTl (1)) 0060 IF (FLIP .GT. 90) FLIP=FLIP-180 0061 IF (FLIP .LT. -90) FLIP=FLIP+180 0062 DIF1-2*FLIP

0063 PULSES (1)*-DIF1

0064 DIF-ABS (DIF1-DIF0)

0065 IF(DIF .GT. DIFMAX) DIFMAX=DIF 0066 DIF=ABS (DIF1-DIFN)

0067 IF(DIF .GT. DIFMAX) DIFMAX=DIF 0068 PRINT *, DIFMAX

0069 RETURN

0070 END

Claims

WE CLAIM:

1. A method of nonintrusively determining displacements of an image slice of a body portion of a patient during a given time interval using a magnetic resonance imaging device, comprising the steps of:

inputting parameters specifying the desired visual characteristics of said image slice;

synthesizing a radio frequency pulse sequence which, when applied to said body portion prior to imaging by said magnetic resonance imaging device, will cause the creation of an image slice having said desired visual characteristics;

applying to the body portion an external magnetic field so as to produce a resultant magnetization;

applying a pre-imaging sequence comprising said synthesized radio frequency pulse sequence and magnetic field gradient pulses to the body portion so as to divide said body portion into a grid of intersecting stripes having altered magnetization, said stripes intersecting at a plurality of intersecting points;

applying an imaging sequence to the body portion to make visible said grid of intersecting stripes;

quantitatively measuring physical displacements of predetermined ones of said plurality of intersecting points during said given time interval; and

providing an image representing said displacements.

2. A method of detecting motion of an image slice of a body portion of a patient using a magnetic resonance imaging device, comprising the steps of:

inputting parameters specifying the desired visual characteristics of said image slice;

synthesizing a radio frequency pulse sequence which, when applied to said body portion prior to imaging by said magnetic resonance imaging device, will cause the creation of an image slice having said desired visual characteristics;

applying to the body portion an external magnetic field so as to produce a resultant magnetization; applying a pre-imaging sequence to the body portion so as to spatially modulate transverse and longitudinal components of the resultant magnetization to thereby produce a spatial modulation pattern, said pre-imaging sequence comprising said synthesized radio frequency pulse sequence and magnetic field gradient pulses;

applying an imaging sequence to the body portion to make visible said spatial modulation pattern of intensity; and detecting motion of said body portion during a time interval between application of said pre-imaging and imaging sequences by examining displacements of the spatial modulation pattern of interest occurring during said time interval.

3. The method as in claim 2, wherein said parameters include stripe thickness in said spatial modulation pattern of intensity.

4. The method as in claim 2, wherein said parameters include weights for relating sharpness and flatness of respective stripes in said spatial modulation pattern of intensity.

5. The method of claim 2, wherein said parameters include a number of radio frequency pulses to be included in said radio frequency pulse sequence.

6. The method of claim 2, wherein said parameters include a size of a transition zone between a striped and a non-striped region in said spatial modulation pattern of intensity.

7. The method of claim 2 wherein said pre-imaging sequence comprises said synthesized radio frequency pulse sequence with said magnetic field gradient pulses interspersed between respective radio frequency pulses of said synthesized radio frequency pulse sequence.

8. The method of claim 2, wherein said pre-imaging sequence comprises said synthesized radio frequency pulse sequence and overlapping magnetic field gradient pulses.

9. The method of claim 2, wherein said synthesized radio frequency pulse sequence causes a selective periodic excitation of said resultant magnetization.

10. The method of claim 9, wherein said synthesized radio frequency pulse sequence comprises respective radio frequency pulses having amplitudes which are not related in accordance with a binomial amplitude function.

11. A method of detecting motion of an image slice of a body portion of a patient using a magnetic resonance imaging device, comprising the steps of:

applying to the body portion an external magnetic field so as to produce a resultant magnetization;

applying a pre-imaging sequence to the body portion so as to spatially modulate transverse and longitudinal components of the resultant magnetization to thereby produce a spatial modulation pattern, said pre-imaging sequence comprising magnetic field gradient pulses and a radio frequency pulse sequence which has been synthesized from parameters specifying the desired visual characteristics of said image slice, whereby when said radio frequency pulse sequence is applied to said body portion prior to imaging by said magnetic resonance imaging device said radio frequency pulse sequence will cause the creation of an image slice having said desired visual characteristics;

applying an imaging sequence to the body portion to make visible said spatial modulation pattern of intensity; and detecting motion of said body portion during a time interval between application of said pre-imaging and imaging sequences by examining displacements of the spatial modulation pattern of interest occurring during said time interval.

12. The method of claim 11, wherein said radio frequency pulse sequence is synthesized in accordance with the following steps:

inputting parameters specifying the desired visual characteristics of said image slice;

calculating a Fourier cosine series in M_z which optimally fits said parameters, where M_z is a magnetization response of said body portion to said radio frequency pulse sequence;

converting said Fourier cosine series into a complex exponential Fourier series;

generating a Fourier series for spinor components of said complex exponential Fourier series of the form: (1+M₂)/2=P² and (1-M_z)/2=Q²;

deconvoluting said Fourier series for spinor components to find P and Q; and

calculating a radio frequency pulse sequence which yields a spinor with P and Q as its spinor components.

13. A method of nonintrusively determining displacements of a body portion of a patient during a given time interval, comprising the steps of:

applying to the body portion an external magnetic field so as to produce a resultant magnetization;

applying a pre-imaging sequence to the body portion so as to divide said body portion into a grid of intersecting stripes having altered magnetization, said stripes intersecting at a plurality of intersecting points;

applying an imaging sequence to the body portion; quantitatively measuring physical displacements of predetermined ones of said plurality of intersecting points during said given time interval; and

providing an image representing said displacements.

14. The method of claim 13, comprising the further step of specifying a triangular grid of said body portion wherein said predetermined ones of said plurality of intersecting points are respective vertices of respective triangles of said triangular grid.

15. The method of claim 14, wherein said displacements measuring step comprises the step of determining which of said intersecting points in said triangular grid at the beginning of said given time interval correspond to said predetermined ones of said plurality of intersecting points in a deformed triangular grid at the end of said given time interval.

16. The method of claim 15, wherein said displacements measuring step comprises the step of comparing said triangular grid at the beginning of said given time interval with said deformed triangular grid at the end of said given time interval.

17. The method of claim 16, wherein said image representing said displacements comprises a comparison of an image representing said triangular grid at the beginning of said given time interval with an image representing said deformed triangular grid at the end of said given time interval.

18. The method of claim 13, wherein said body portion is the patient's heart wall.

19. An apparatus for noninvasively determining displacements of a body portion of a patient during a given time interval, comprising:

means for applying to the body portion an external magnetic field so as to produce a resultant magnetization;

imaging means for applying a pre-imaging sequence to the body portion so as to divide said body portion into a grid of intersecting stripes having altered magnetization, said stripes intersecting at a plurality of intersecting points, and for applying an imaging sequence to the body portion;

means for quantitatively measuring physical displacements of predetermined ones of said plurality of intersecting points during said given time interval; and a display for displaying an image representing said displacements.

20. The apparatus of claim 19, further comprising means for specifying a triangular grid of said body portion wherein said predetermined ones of said plurality of intersecting points are respective vertices of respective triangles of said triangular grid.

21. The apparatus of claim 20, wherein said displacements measuring means comprises means for determining which of said intersecting points in said triangular grid at the beginning of said given time interval correspond to said predetermined ones of said plurality of intersecting points in a deformed triangular grid at the end of said given time interval.

22. The apparatus of claim 21, wherein said displacements measuring means comprises means for comparing said triangular grid at the beginning of said given time interval with said deformed triangular grid at the end of said given time interval.