UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 EXPERIMENTAL STUDY OF THERMOACOUSTIC IMAGING SYSTEM By GUAN XIN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSO PHY UNIVERSITY OF FLORIDA 2012 PAGE 2 2 2012 Xin Guan PAGE 3 3 To my parents and my wife PAGE 4 4 ACKNOWLEDGMENTS Thanks to my advisor, Professor Henry Zmuda, for his kind help and excellent guidance for my research work. Without him, I could not have the achieve ment today Thanks to Professor Jian Li, for her kind help to assist me in understanding the Robust Capon Beamforming algorithm, which is the major algorithm used and implemented in my thesis. Thanks to Professor Mark Sheplak for his helpful advices on aco ustic measurement and theoretical analysis for thermoacoustic imaging. Thanks to Professor Louis Cattafestta for his kind help on data measurement and all the discussions we have had PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 ABSTRACT ................................ ................................ ................................ ................... 11 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 13 1.1 Thermoacoustic Imaging ................................ ................................ ................... 13 1.2 Research Background ................................ ................................ ...................... 16 ................................ ......... 16 ................................ .................... 22 1.3 Thesis Organization ................................ ................................ .......................... 25 2 THERMOACOUSTIC SIGNAL GENERATION ................................ ....................... 28 2.1 CW Microwave Stimulation ................................ ................................ ............... 28 2.1.1 Diebold Solution ................................ ................................ ...................... 28 2.1.2 Resonant Thermoacoustic Signal Generation ................................ ......... 32 2.1.3 RF Amplitude Mod ulation ................................ ................................ ........ 35 2.1.4 CW Stimulation Experiment ................................ ................................ ..... 37 2.2 Pulsed Microwave Excitation ................................ ................................ ............ 40 3 IMAGING ALGORITHM AND FPGA IMPLEMENTATION ................................ ...... 43 3.1 Digital Beamforming ................................ ................................ .......................... 43 3.2 Robust Capon Beamforming Algorithm ................................ ............................. 46 3.2.1 Standard RCB For Narrowband Signal ................................ .................... 46 3.2.1.1 Capon Beamforming ................................ ................................ ...... 46 3.2.1.2 Ro bust Capon Beamforming ................................ .......................... 48 3.2.1.3 Numerical simulation ................................ ................................ ...... 49 3.2.2 Time Delayed RCB for Wideband Signal ................................ ................. 53 3.2.3 Computational Complexity for Different Algorithm Implementation ......... 55 3.3 FPGA Implementation of RCB ................................ ................................ .......... 58 3.3.1 Implementation of Standard RCB ................................ ............................ 60 3.3.1.1 System block diagram ................................ ................................ .... 60 3.3.1.2 Complex Multiply and Accumulate Unit ................................ .......... 61 3.3.1.3 Extended Cyclic Jacobi Method ................................ ..................... 62 3.3.1.4 Lagrange Multiplier Solver ................................ ............................. 65 3.3. 1.5 Adaptive Weight Calculation ................................ .......................... 66 PAGE 6 6 3.3.1.6 FPGA implementation ................................ ................................ .... 67 3.3.2 Implementation of Wideband RCB ................................ .......................... 68 3.3.2.1 Cyclic Jacobi Method for Eigendecomposition ................................ ..... 69 3.3.2.2 FPGA implementation of wideband RCB ................................ ....... 72 4 EXPERIMENTAL RESULTS OF THERMOACOUSTIC IMAGING SYSTEM .......... 73 4.1 Experimental Step up of TAI System ................................ ................................ 73 4.1.1 Introduction ................................ ................................ .............................. 73 4.1.1.1 High Power Subsystem ................................ ................................ .. 74 4.1.1.2 Transducer array subsystem ................................ .......................... 77 4.1.1.3 DAQ s ubsystem ................................ ................................ ............. 79 4.2 Thermoacoustic Signal Generation Experiment ................................ ................ 83 4.3 Thermoacoustic Imaging Experiment ................................ ................................ 8 7 4.4 Imaging Results ................................ ................................ ................................ 90 5 CONCLUSION ................................ ................................ ................................ ........ 94 5.1 Contribution ................................ ................................ ................................ ...... 94 5.2 Future Work ................................ ................................ ................................ ...... 95 LIST OF REFERENCES ................................ ................................ ............................... 97 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 101 PAGE 7 7 LIST OF TABLES Ta ble page 1 1 Typical Dielectric Properties of Breast Tissue ................................ .................... 15 2 1 Peak Power Levels for different experiment setup ................................ ............. 39 3 1 Resource Utilization for RCB Implementation ................................ .................... 67 3 2 Resource Utilization for Wideband RCB Implementation ................................ ... 72 PAGE 8 8 LIST OF FIGURES Figure page 1 1 Schematic of Thermoelastic effects ................................ ................................ .... 15 1 2 Schematic of Experimental Setup for Professor L ................. 16 1 3 Schematic of the circular measurement configuration ................................ ........ 20 1 4 Schematic of the new combined imaging system ................................ ............... 22 1 5 ........................ 22 1 6 ................ 25 2 1 Amplitude frequency response with density ratio = 1 ................................ ......... 33 2 2 Frequency spectrum with density ratio scanned from 1 to 5 ............................... 34 2 3 Frequency spectrum with density ratio scanned from 1 to 0.6 ............................ 34 2 4 Experiment Schematic for CW Excitation ................................ ........................... 37 2 5 Output Signal on Spectrum Analyzer with 150 kHz Input Signal ........................ 38 2 6 Numerical simulation for pulsed microwave stimulation ................................ ..... 42 3 1 Adaptive Antenna Array ................................ ................................ ...................... 43 3 2 Beamforming Network ................................ ................................ ........................ 45 3 3 System Diagram of Adaptive Beamformer ................................ ......................... 45 3 4 Signal Comparison ................................ ................................ ............................. 51 3 5 SCB Filtered Signal ................................ ................................ ............................ 51 3 6 RCB Filtered Signal ................................ ................................ ............................ 52 3 7 Beam Pattern Plot ................................ ................................ .............................. 52 3 8 Power Estimates Plot ................................ ................................ ......................... 53 3 9 Illustration of time delayed RCB ................................ ................................ ......... 54 3 10 Schematic for Computational Complexity Analysis for Matrix Multiplication ....... 56 3 11 Run time Performance Com parison ................................ ................................ ... 57 PAGE 9 9 3 12 System Block Diagram ................................ ................................ ....................... 60 3 13 Datapath of Complex MAC Module ................................ ................................ .... 61 3 14 Datapath for Eigenvalue Decomposition Module ................................ ................ 63 3 15 Datapath for CORDIC Processor ................................ ................................ ........ 64 3 16 Datapath for Diagonal Processing Element ................................ ........................ 65 3 17 Datapath of Lagrange Multiplier Solver ................................ .............................. 65 3 18 Datapath for Eigenvalue Decomposition Module ................................ ................ 70 3 19 Datapath for CORDIC Processor ................................ ................................ ........ 71 3 20 Datapath for Diagonal Processing Element ................................ ........................ 71 3 21 Dat apath of Off Diagonal Processing Element ................................ ................... 71 4 1 TAI System Block Diagram ................................ ................................ ................. 73 4 2 High Power Subsystem ................................ ................................ ...................... 74 4 3 Control Panel of Pulse Generator ................................ ................................ ....... 75 4 4 Summary of Operational Conditions ................................ ................................ ... 76 4 5 Transducer ................................ ................................ ................................ ......... 77 4 6 Motor and Transducer System ................................ ................................ ........... 78 4 7 Scanning Schematic ................................ ................................ ........................... 79 4 8 DAQ Subsyste m Schematic ................................ ................................ ............... 80 4 9 Pulser/Receiver Operating Parameters ................................ .............................. 80 4 10 Configuration for NI 5122 Digitizer ................................ ................................ ...... 81 4 11 Trigger Setting for NI 5122 Digitizer ................................ ................................ ... 82 4 12 Control Panel ................................ ................................ ................................ ...... 83 4 13 Tank without phantom (right), T ank with phantom (left) ................................ ...... 84 4 14 Relative Positions between Transducer and Phantom ................................ ....... 84 4 15 Thermoacoustic Signal in the First Experimen t ................................ .................. 85 PAGE 10 10 4 16 Zoom in View of the Thermoacoustic Signal ................................ ...................... 86 4 17 Thermoacoustic Signal in the Second Experiment ................................ ............. 87 4 18 Thermoacoustic Signal in the Third Experiment ................................ ................. 87 4 19 New Experimental Tank ................................ ................................ ..................... 88 4 20 Schemati c of the phantom ................................ ................................ .................. 89 4 21 Thermoacoustic signals at 0, 120, 240 degree ................................ ................... 89 4 22 Imaging Result with 4mm 15% salt phantom using RCB ................................ .... 90 4 23 Imaging Result with 4mm 15% salt phantom using SCB ................................ .... 90 4 24 Imaging Result with 4mm 15% salt phantom using DAS ................................ .... 91 4 25 Imaging Result with 2mm phantom using RCB ................................ .................. 91 4 26 Imaging Result with 4mm phantom using RCB ................................ .................. 92 4 27 Imaging Result with 6mm phantom using RCB ................................ .................. 92 PAGE 11 11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of D octor of Philosophy EXPERIMENTAL STUDY OF THERMOACOUSTIC IMAGING SYSTEM B y X in G uan M ay 201 2 Chair: Henry Zmuda Major: Electrical Engineering Thermoacoustic Imaging is a relatively new medical imaging modality, which combines the benefits of pure mic rowave imaging and conve ntional ultrasound imaging, achieving both their high image contrast and high spatial resolution respectively In operation, the soft tissues in the imaging area are exposed to microwave pulses with a high peak power. U pon irradiati on the tissues will absorb the electromagnetic energy, expand in volume and gen erate a thermoacoustic signal propagating in all directions. A transducer system picks up the generated thermoacoustic signal, which is then amplified, digitized and transferr ed to the computer. An image reconstruction program ru n ning in the computer solves an inverse problem to obtain the image of the microwave absorption distribution inside the ti ssue. C ancerous tissue has a relatively strong electromagnetic absorption due to its higher water content [1], and distinguishes itself from other tissues in the image. An experimental study of the thermoacoustic imaging system is performed in this thesis, which includes building the prototype system, conducting the initial experime nts for generating thermoacoustic signals and obtaining images of different phantoms with this system. At the beginning, CW (Continuous Wave) microwave stimulation was PAGE 12 12 employed to find the resonance of the phantom, and finally fail ed for various reasons. Then pulsed microwave radiation was used and successfully gene rated the thermoacoustic signal A single transducer, coupled with a motor system, was used to simulate a multiple transducer array for data acquisition Based on the time delayed Robust Capon Beamforming (RCB) algorithm, several images were acquired from different phantoms. A Field Programmable Gate Arrays (FPGA) device was used to investigate a hardware implementation of the RCB algorithm. Two implementations were achieved, with one realizin g the narrowband RCB used in smart antenna system s and another realizing the wideband RCB used in medical imaging. The possibility of making the thermoacoustic imaging system run in real time is studied in this work. PAGE 13 13 CHAPTER 1 INTRODUCTION 1.1 Ther moacoustic Imaging Breast cancer is considered as the most common type of non skin cancer in women, and the fifth most usual cause of cancer death worldwide among others The successful treatment of breast cancer greatly depends on the early detection of t he malignant tissue The current primary screening method for breast cancer is X ray mammography. X ray is a kind of ionizing radiation and could be dangerous for human body in frequent usage During the diagnosis, the patient s need to co mpress their brea sts and suffer a lot of pain. In addition, the dense glandular tissue, connectiv e tissue, and the region near chest wall are difficult to image due to their high density [2]. Ultrasound imaging is commonly used as a complementary method to mamm ography for breast cancer diagnosis. However the poor soft tissue contrast between benign and malignant tumor in breast limits the application of ultrasound imaging in breast cancer diagnosis [2]. Magnetic Resonance Imaging (MRI), which is also used as a complementary imaging technology for breast imaging, could de liver improved tissue contrast Although t he sensitivity of MRI in detecting malignant tissue has been excellent, the most exp ensive cost of MRI imaging among all other breast imaging modalities hinders the a doption of MRI as a routine screening method for b reast cancer [2]. Recently, thermoacoustic imaging has gained a lot of attention as a promising breast cancer screening method Firstly microwave radiation is non ionizing and much healthier for the human body than other ionizing radiatio n, such as X ray. Secondly t hermoacoustic imaging can deliver t he high soft tissue contrast together with high PAGE 14 14 spa tial resolution In addition, thermoacoustic imaging system is more cost effective than MRI. Thermoaco ustic imaging is a relatively new imaging modality [3], which irradiates the soft tissue by electromagnetic radiation, such as optical or microwave radiation, records the ultra sound signal generated by soft tissue, and obtains the image about distribution of electromagnetic absorption inside the soft tissue through various image reconstruction algorithms. The imaging technology using optical radiation is usually called photoacoustic imaging to show differences from microwave induced thermoacoustic imaging. Conventional microwave imaging depicts the distribution of dielectric property inside soft tissue by impinging microwave beam onto the target and measuring the back scattered microwave signal Cancerous tissue has higher permittivity and conductivity (5 to 10 times) than normal tissue due to higher water content in malignant cells as seen in Table 1 1 [4], and generates larger scattered signal which c ould be exploited to figure out tumor locations. Microwave imaging has better image contrast than ultrasou nd imaging, but poor resolution because of long wavelength of microwave. On the other side, ultrasound imaging is able to achieve the resolution in millimeter range because of the short wavelength but has lower image contrast due to the small difference i n acoustic property between tumor and normal tissue. Thermoacoustic imaging combines the benefits of both microwave imaging and ultrasound imaging, and therefore has high image contrast due to the large difference in conductivity between tumor and normal t issue, together with high spatial resolution due to the short wavelength of thermoacoustic signal [5]. PAGE 15 15 Table 1 1. Typical Dielectric Properties of Breast Tissue (Source : PhD Thesis University of Florida, Gainesville, FL, 2007. ) TISSUES DIELECTRIC PROPERTIE S PERMITTIVITY CONDUCTIVITY IMMERSION LIQUID 9 0 CHEST WALL 50 7 SKIN 36 4 FATTY BREAST TISSUE 9 0.4 NIPPLE 45 5 GLANDULAR TISSUE 11 15 0.4 0.5 TUMOR 50 4 Figure 1 1. Schematic of Thermoelastic effects In microwave induced th ermoacoustic imaging system, pulsed microwave radiation is used to heat the tissue in sho rt duration. After being heated the tumor expands and generates the acoustic pre ssure wave as illustrated in Fig ure 1 1 The acoustic signal is picked up by the ultrasonic transducer or transducer array surrounding the imaging area. After amplific ation and filtering, the signal is digitized by Analog to Digital (AD) circuitry and proc essed by image reconstruction algorithm to obtain the PAGE 16 16 image of conductivity distribution inside the tissue. The peak power of the microwave pulse shoul d be high enough to generate thermoacousti c signal with sufficient signal to noise ratio (SNR) in order t o be observable. The recorded signal from AD card is usually avera ged multiple times to reduce noise level and increase SNR 1.2 Research Background The earliest research about thermoacoustic imaging was published in 1981 by Theodore Bowen [6]. In his p aper, issues about t hermoacoustic signal generation and detection were investigated, an d the possibility of imaging tissue characteristics by thermoacoustic signal was confirmed James C. Lin published another paper in 1988 [7], in which he discussed some system design issues for thermoacoustic imaging such as signal detection, signal conditioning and data conversion techniques. 1.2.1 Resear ch i In the first paper published in 1999 [8], the experimental setup and initial re sults were illustrated in det ails. The schematic is as below Figure 1 (Source: induced acoustic Review O f Scientific Instruments vol. 70, no. 9, pp. 3744 3748, Sep 1999 ) PAGE 17 17 In this first experimental setup, they used a 9 GHz microwave source with 10 kW peak power. The thermoacoustic signal was clearly recorded with different power levels and different phantoms used in experiments Two focused transducers were utilized to receive acoustic signal with one centered at 1 MHz bandwidth 0.6 5 MHz and focal length 2.5 cm at 1MHz and t he other centered at 3.5 MHz, bandwidth 2.5 MHz and focal length 1.8 cm at 3.5 MHz In this paper, the initial image was acquired for a phantom. T he time resolved therm oacoustic signal recorded by transducer was directly converted to the one dimensional depth resolved image As a result, n o image recon struction algorithm was involved for imaging In the following papers published in 2000 [1][9], the microwave frequency was changed to 3 GHz because of the less attenuation in soft tissue. As a result, the image obtained is much better than the original image acquired from 9 GHz microwave sou rce As same as the initial experiment two focused transducers with different center frequency and bandwid th were utilized to receive thermoacoustic signal. In the first case, the transducer with center frequency 1 MHz and bandwidth 0.6 MHz was deployed f or reception The axial resolution was estimated to b e 3.3 mm which coincide with measurement result [1]. In the second case, the transducer with center frequency 3.5 MHz and bandwidth 2.5 MHz was utilized to receive acoustic signal, and the axial resoluti on was estimated to b e 1.4 mm which is also close to experimental result [1]. In experiments, the axial resolutio n is primarily determined by the pulse w idth of microwave source and frequency spectrum of the transducer. Shorter microwave pulse and wide ban d transd ucer should result in better axial resolution. The lateral resolutio n is primarily determined by numer ical aperture of the transducer, which is a m easure of the capability for transducer to collect ultrasonic PAGE 18 18 signal and expose the de tails of soun d source at fixed distance In experiments, the lateral resolution was e stimated to 2.1 mm which matches well with measurement result [1]. In p aper [10] published in 2001, advanced signal processing t echnique was used to improve axial resolution and latera l resolution in the original image. The axial resolution is determined by the frequency band of the aco ustic signal, while based on analysis about numerical aperture of transducer, the lateral resolution is inversely proportional to the frequency component of the acoustic signal [10]. If we simply use a band pass filter to get rid o f the low frequency noise in acoustic signal, the lateral resolution will be improved, but the axial resolution is sacrificed due to the reduce d frequency band of acoustic signal As a result, a reshaping filter was used for filtering in order to broaden the frequency band and emphasize the high frequency components in acoustic signal. Therefore the axial resolution and lateral resolution were improved at the same time with higher computational complexity. Until then focused trans ducer was deployed to detect ultrasound signal and the imag e was obtained directly from one dimensional time resolved thermoacoustic signal, meaning that no image reconstruction method was utilized. In p aper [11] published in 2001, synthetic aperture method was used for image reconstruction. In experiment, unfocused transducer with center frequency 2.25 MHz was used U nfocused transducer is able to produce unfocused ultrasound beam in near field region, a nd also receive ultrasound signals from all directions with a receiving directivity pattern. Therefore it can be used to record ultrasound signal s generated by different parts of phantom and propagating at different directions The simple and robust synthe tic aperture method PAGE 19 19 (Delay and Sum) was used to reconstru ct the image accordingly. In planar measurement configurat ion with focused transducer, boundaries nearly perpendicular to the acoustic axis of transducer can be imaged clearly, because large part of the thermoacoustic signal generated propagates along the acoustic a xis of transducer and is received by tr ansducer. On the other hand, boundaries paralleling to the acoustic axis canno t be imaged, because most of thermoacoustic signal generated propagates perpendicular to the acoustic ax is and cannot be received by transducer. This problem was solved by usage of unfocused transducer and synthetic aperture method for image reconstruction. S ynthetic aperture method could perform even better in circular or sph erical measurement configuration. In p aper [12] published in 2001, focused transducer was used in multi sector scanning to scan the illuminated sample at all directions on a circular trajectory. The image reconstruction method was still the direct method a s before and the image was also successfully acquired. In paper [13] and paper [14] published in 2002, another image reconstruction algorithm filtered back projection algorithm, was proposed and used to obtain the image. Filtered back projection algorith m is the algorithm which tries to reconstruct the original image from in computed tomography technology. In thermoacoustic imaging research, this algorithm is usually consulted for image reconstruction. At first, the detailed analytical derivation was illustrated [13], which leads to an approximated solution under two assumptions The first assumption is that the detection radius is usually much larger than wavele ngth of thermoacoustic sig nal and the s econd assumption states that low frequency component of thermoacoustic signal does not significa ntly contribute to spatial PAGE 20 20 resolution [13]. As a result, the modified filtered back proj ection algorithm was used as image rec onstruction method t o depict electromagnetic ab sorption distribution inside phantom in [14] The experiment was p erformed in circular measurement configuration, and unfocused transducer with center frequency 2.25 MHz was deployed to receive thermoacoustic signal. The schemati c is showed in Figure 1 3. M icr owave pulse was used to illuminate the tissue In image reconstruction, the filtered back projection algorithm was used In order to get spatial resolution, point spread function (PSF) was calculated. For cutoff frequency of 4 MHz, 2 MHz and 1 MHz, the full width half maximum (FWHM) of PSF was obtained as 0.4 mm, 0 .9 mm and 1.5 mm In p aper [17] published in 2003, time domain reconstruction algorithm for various measure configuration, including planar, cylindrical and spherica l geometry, was unified in one analytical formulation. Figure 1 3. Schematic of the circular measurement configuration (Source : M. Xu and L. microwave induced thermoacoustic tomography: Filtered backprojection in a circular measurement Med. Phys ., vol. 29, no. 8, pp. 1661 1669, Aug 2002. ) PAGE 21 21 I mage reconstru ction algorithms mentioned are all based on the assumption that surrounding tissue is acoustically homogeneous, which means that the sound speed and acoustic property sho uld be constant across the entire tissue. Ho wever, in real applications, sound speed may vary from 1430 m/s to 1570 m/s around the assumed 1500 m/s [18], and as a result the phase and amplitude of therm oacoustic signal received at transducer is usually di storted due to heterogeneity of acoustic property across the tissue. In paper [18] published in 2008, a new algorithm, Adaptive and Robust Methods of Reconstruction (ARMOR), was used to perform image reconstruction in thermoacoustic im aging. This algorithm allows amplitude and ph ase of thermoacoustic signal have some uncertainty, and uses the Robust Capon Beamforming algorithm and peak searching technique to cope with thi s uncertainty. As a result, image reconstruction became more robust. This ARMOR algorit hm was further extended to exploit one more dimension microwave stimulation frequency to reduce more noise and interferences in paper [1]. In paper [19] published in 2008, a new device combining thermoacoustic imaging and photo acoustic imaging was pro posed and developed. The schematic of the new device is showed in Figure 1 4. T hermoacoustic imaging can expose the tissue part having high electromagnetic absorption, such as cancerous tissue due to high water content, while photo acoustic imaging ca n dep ict the tissue part showing stronger optical absorptio n such as blood vessel. The combination of both medical imaging modalities could enhance early stage can cer detection, and also have several benefits. Firstly, it will reduce the imaging time and is co st effective S econdly acquiring two PAGE 22 22 images on the same setup avoids moving and realigning patient all over again [19]. Imaging results with good quality were achieved on this combined system. Figure 1 4. Schematic of the new combined imaging system (So urce : M. Pramanik, G. Detection System Combining Both Thermoacoustic and Photoacoustic Med. Phys. vol. 35, no. 6, pp. 2218 2223, Jun 2008. ) 1.2.2 Research i n Professor K The first thermoacoustic image was obtained by Dr Robert A. group. In paper [3] published in 1999, the prototype system was presented wit h schematic showed below Figure 1 5 Schematic of Experimental Setup for Professor (Source : R. Radiology vol. 121, no. 1, pp. 275 278, Apr 1999. ) PAGE 23 23 In this prototype system, 434 MHz microwave pulse with pulse wid th 0.5 us and peak power 25 kW was used to illumi nate the tissue. Repetition rate for microwave pulse was set to 4000 Hz, provi ding 50 W of average power. M icrowave power was fed into the imaging tank by a helical RF antenna, and thermoacou stic signal was picked up by transducer array mounted on the rotational hemisphere bowl, which rotates into many different circular locations to hav e transducer array receiving acoustic signals in different positions during imaging process After received by transducer, t hermoacoustic signal was amplified and averaged for multiple times in order to increase signal to noise ratio (SNR). The averaged signal was then digitized and transferred to computer. F iltered ba ck projection algorithm was utilized to perfor m three dimen s ional image reconstruction. Reconstructed image of exc ised porcine kidney was showed [3] In paper [20] published in 1999, technical consideration s were reported regarding design and implementation of thermoacoustic computed tomography. The approximated Radon transform was used to derive filtered back projection algorithm for image reconstruction, which lends itself into easy implementation in computer program or special hardware. According to [20], 434 MHz RF radiation was chosen as stimulation frequenc y based on several reasons. Firstly, electromagnetic absorption for soft tissue is peaked around this fre quency range, and 434 MHz is assigned by FCC for medical usage. S econdly, relative SNR of thermoacoustic signal was simulated to peak in the range of 2 00 600 MHz. At this low frequency, penetration depth of microwave is adequate for medical imaging applicat ion [20]. The positioning of 64 element transducer arr ay was explained in details Based on the mounting positio ns of transducer relative to imaging bowl, the rotational movement of the bowl was determined to achieve PAGE 24 24 maximum coverage for ther moacoustic signal reception. RF heating effects for tissue was also experimented to show that less than 1 degree temperature rise was observed during experiment w hich suggests the safety of human being for performing thermoacoustic imaging. The data acquisition scheme and transducer characteristics were also introduced in that paper. In p aper [21] published in 2000, imaging principles of thermoacoustic CT was illu strated in details. R econstruction methodology based on approximate Radon transform was derived and demonstr ated in computer simulation. S patial resolution was studied extensively using computer simulation. The best spatial resolution near the center axis was showed to be 1.2 mm under 0 .5 us microwave pulse stimulation The limiting factors for spatial resolution were pointed out as the duration of microwave pulse, frequency response of transducer together with ass ociated electronics, and physical size of t he face of transducer [21]. D ifferent effects on spatial resolution from different factors were demonstrated by computer simulation. After that, pilot study of the application of thermoacoustic imaging system on breast cancer detection was conducted, and r esults were published in the paper in 2001 [22]. Later on, new systems were proposed, built and experimented [23][24]. The schematic is shown below in Fig ure 1 6 Wa veguide array was used to feed microwave pulse into imaging tank in order to make the elect romagnetic field distribution more uniform than previous situation Planar transduc er array was used to replace discrete transducer with improved sensitivity, frequency response and immunity to electromagnetic induction (EMI). Another imaging system was de scribed in 2010 paper by Kruger, et. al [25], which combined thermoacoustic imaging, photoacoustic imaging PAGE 25 25 and ultrasound imaging modalities together. The co registered image by different imaging technologies was presented in that paper. Figure 1 6 Schem atic of New Experimental (Source : Medical Imaging 2002: Physics of Medical Imaging pp. 521 525, May 2002. ) There are also many other research g roups working on this topic group proposed and develo ped fast thermoacoustic imaging system with large transducer array and parallel data acquisition system to reduce data collection time, which is described in their paper published in 2006 [26] They performed research for enhancing thermoacousti c image by microwave contr ast agent in [27]. In [28], the feasibility of detecting renal calculi using thermoacoustic imaging technique was investigated 1.3 Thesis Organization In Chapter 2, theore tical analysis of thermoacoustic signal generation is presented. Firstly, the analysis for thermoacoustic signal generation under CW microwave stimulation is presented. Diebold solution is showed to depict frequency PAGE 26 26 response of a spherical phantom upon hea ting which points out the possibility of resonance of phantom in time varying heating condition. A mplitude modulation technique for RF signal is then presented, in order to provide baseband ti me varying heating function. After introducing negative experim ental results for thermoacoustic signal generation under CW microwave stimulation, r easons for the discrepancy between theoretical analysis and experimental results are d iscussed. At last, t heoretica l analysis for pulsed microwave induced thermoacoustic si gnal generation from Professor search is quoted to complete the analysis for thermoacoustic signal generation In Chapter 3, the major imaging algorithm used in our thermoacoustic imaging appl ication is discussed. First, original narrowband Robust Capon Beamforming algorithm [29] is outlined, together with simulation results. Second, extended wideband Robust Capon Beamforming algorithm used in our medical imaging application is intro duced [18]. After discussing software implementation of the algorithm, proposed FPGA implementation of the algorithm is presented. Two implementations ar e developed respectively for narrowband and wideband Robust Capon Beamformer. The current implementation is still an initial version, which needs to be improved t o satisfy the real time requirements. In Chapter 4, current experimental setup is introduced as three subsystems, high power subsystem, transducer array subsystem and DAQ subsystem. E xperimental results of our thermoacoustic imaging system are presented, including typical thermoacoustic signal waveform and initial imaging results. Different sets of phantom PAGE 27 27 are experimented, and resu lts are compared The current imaging functionality and limitations are pointed out in this Chapter. In Chapter 5, contributio ns of our thermoacoustic imaging research are identified In addition, future works are presented. PAGE 28 28 CHAPTER 2 THERMOACOUSTIC SIGNA L GENERATION When electromagnetic or o ptical radiation impinges on soft tissue, energy is absorbed leading to heating effects and subsequent expansion of the tissue gives rise to acoustic wave which can be picked up by transducer This is called thermoacoustic or phot oacoustic effects, and has been the research topic for a long time in order to utilize this phenomenon for imagi ng purpose or characterization of different particulate matter. In experiments, usually short pulse of electromagnetic or optical radiation with high enough energy density is used to trigger observable thermoacoustic or photoac oustic signal which is above noise floor. 2.1 CW Microwave Stimulation 2.1.1 Diebold Solution G J. Diebold, etc, performed research on photoacoustic effect s gener ated by a spherical droplet in fluid. In their p aper published in 1988 [30], the analytical solution of the generated acoustic wave was presented in frequency domain and time domain M ajor derivati ons are listed as follow. If electromagnetic radiation is modulated with baseband frequency and the travel time of electromag netic wave can be neglected, energy absorbed by phantom can be formulated in equation 2.1 as below: = 0 (2 1) w here 0 is the illumination intensity, is the absorption coefficien t of phantom material and is theamount of energy absorbed per unit volume per unit time. Under thermal confinement condition ( the duration of heating pul se sh ould be short enough so that heat conduction to surround ing environment could be neglected) the inhomogeneous wave equation inside the sphere is presented below: PAGE 29 29 2 1 2 2 2 = (2 2) w here is th e pressure wave generated by heating is the sound speed in phantom, is thermal expansion coefficient and is specific heat of phantom material Substituting equation (2 1) into equation (2 2), we can get the following equation (2 3): 2 1 2 2 2 = 0 (2 3) Due to symmetrical condition, we can further assume spherically symmetric solution as = ( ) and plugging into equation (2 3) yields: 2 2 + 2 + 2 = 0 (2 4) w here = is the wave number inside phantom. With particular solution = 0 2 and homogeneous solution = sin the full solution for pressure wave inside phantom body is listed below, = 0 2 1 + sin (2 5) and the full solution for pressure wave outside phantom body is showed as following, = ( ) (2 6) w here = is the wave number outside phantom, and is sound speed of surrounding liquid. T erms in the equations above are unknown coefficients which need to be determined by apply ing boundary condition Until now, we have got the analytical expression of pressure wave inside the spherical phantom a nd outside the phantom. Then pressure and normal velocity are matched at interface. Regarding acceleration, Eul = should be PAGE 30 30 consulted, where is density, is axial velocity and is pressure. After setting pressure and normal velo city from inside and outside phantom to be equal at interface two unknown coefficients and can be s olved as below: = ( 1 ) cos + 1 sin sin (2 7) = 0 2 1 + ( 1 ) cos + 1 sin sin sin (2 8) w here is the density of phantom material, is the density of surrounding liquid, is radius of the spherical phantom. Substituting equations (2 7) and (2 8) into equations (2 5) and (2 6) yields the pr essure wave solution inside and outside phantom as below: = 0 1 + ( 1 ) 1 sin cos + sin sin (2 9) = 0 sin cos 1 sin cos + sin ( ( ) ) (2 10) In order to further simplify the expression, we introduce following dimensionless terms: = dimensionless time = ( ) retarded time = wave vector = relative distance from sphere center Plugging in the terms above, inner and outer solution of pressure wave can be formulated as follow: PAGE 31 31 = 0 1 + 1 sin 1 sin cos + sin (2 11) = 0 sin cos 1 2 1 sin cos + sin (2 12) By using time harmonic function for heating, thermoacoustic pressure wave is solved anal ytically in frequency domain. Impulse response of thermoacoustic process in time domain could be found by performing Inverse Fourier Transform on frequency domain solution. T hermoacoustic pressure wave under radiation with arbitrary temporal profile could be found by convolving impulse response with temporal radiation profile. D etailed derivation is presented in [30]. 2.1.1.1 Harmonic Solution in Cylindrical Coordinate System The frequency domain solution in previous section is achieved in spherical coordi nate system while their experiment was conducted for the spherical droplet in a fluid. However our thermoacoustic imaging experiment was performed in cylindrical coordinate system in which the phantom was made into tall cylinder. Therefore we parallel the theoretical analysis above for cylindrical coordinate system to make the derivation complete. After we substituting the heating function in equation (2 1) into the wave equation in (2 2), we end up with the differential equation (2 3). If we assume symmetr ical condition and infinite cylindrical source, equation (2 3) could be simplified as: 2 2 + + 2 = 0 (2 13 ) where = is the wave number inside phantom. With particular solution = PAGE 32 32 0 2 and homogeneous solution = 0 ( ) the full solution for pressure wave inside phantom body is listed below, = 0 2 1 + 0 ( ) (2 14 ) and the full solution for pressure wave outside phantom body is showed as fol lowing, = 0 2 ( ) (2 15 ) w here = is the wave number outside phantom, and is sound speed of s urrounding liquid. Terms in the equations above are unknown coefficients which need to be determined by applying boundary condition. After matching the pressure and acceleration at the boundary, we could solve the unknown coefficients as foll ow, = 1 ( 2 ) ( ) 1 0 ( 2 ) 0 1 ( 2 ) (2 16) = 0 2 1 ( ) 1 0 ( 2 ) 0 1 ( 2 ) (2 17) Substituting equation (2 17) into equation (2 15), we could find the full solution for pressure wave outside phantom as equation (2 18). = 0 2 1 1 0 ( 2 ) 0 1 ( 2 ) 0 2 ( ) (2 18) The impulse response of the thermoacoustic process in cylindrical coordinate system could be obtained by applying the Inverse Fourier Transform on the frequen cy domain solution in equation (2 18). 2.1.2 Resonant Thermoacoustic Signal Generation Accor ding to the analytical solution obtained in previous section, computer simulation has been conducted to study the frequency response of thermoacoustic pressure wa ve generation process. D ifferent frequency responses are showed in the PAGE 33 33 plot under varying density ratio which is obtained through dividing the density of phantom by the density of surrounding liquid. The case when density ratio equals one is used as refe r ence for other situations. Frequency response of the referential situation is illustrated below in Fig ure 2 1. Multiple peaks appear at a series of resonant frequencies, which give large amplitude amplification, and are relatively b roa der compared to peaks in larger density ratio cases In Fig ure 2 2, a scan of density ratio from 1 to 5 with step size 1 has been performed, and frequency responses are illustrated. R esonant pe aks are becoming narrower as density ratio becoming larger, and at the same time res onant frequencies have been shifted a little bit. In Fig ure 2 3, a scan of density ratio from 0. 6 to 1 with step size 0.1 has been performed, and frequency responses are illustrated. Figure 2 1. Amplitude frequency response with density ratio = 1 PAGE 34 34 Figu re 2 2. Frequency spectrum with density ratio scanned from 1 to 5 Figure 2 3. Frequency spectrum with density ratio scanned from 1 to 0. 6 Based on frequen cy response presented above, the resonant thermoacoustic imaging idea is developed. L ow frequency he ating function sho uld be utilized to PAGE 35 35 stimulate phantom, and the frequency should be swept in certain frequency band to locate the actual resonant frequency. If the res onant frequency was located, lower input power CW stimulation should be applied at this r esonant frequency, and larger output signal shou ld be expected due to resonant stimulation. In next section, the technique used to generate low frequency heating fu nction is illustrated, while CW stimulation experiment is described in following section. 2 .1.3 RF Amplitude Modulation In ord er to get low frequency time varying heating function, we need to apply ampli tude modulation technique on high frequency microwave carrier signal, because baseband RF signal (around 1 2 MHz) will be reflected back from t he interface between waveguide and experimental tank due to unmatched electromagnetic impedance s on both sides. The relationship between the electric field and the energy absorbed by phantom indicated by equation (2 1 9 ) will generate the low frequency heat ing fu nction, which will serve as the stimulation for phantom. Higher frequency components of heating function should have no effects, because the mechanical expansion and contraction of phantom will not be able to res pond at such high frequency. D etailed deri vation is illustrated below. SAR (Specific Absorption Rate) can be formulated as equation (2 1 9 ). = 2 ( ) ( ) (2 1 9 ) w here is the conductivity of phantom material, ( ) is the intensity of electric field imposed on phantom and ( ) is the density of phantom. represents the ener gy absorbed per unit mass by phantom. Temperature in crease can be obtained from thermodynamics principle as follow: PAGE 36 36 ( ) = ( ) = 2 ( ) ( ) = ( ) 2 ( ) (2 20 ) The temporal formulation of electromagnetic field can be expressed as below assuming uniform heating across phantom body, = 0 1 + (2 21 ) = 0 1 + cos cos where m(t) is the low frequency modulation function. Substituting equation (2 21 ) into equa tion (2 20 ) yields: ( ) = 0 2 ( 1 + cos ) cos 2 in which ( 1 + A cos ( m t ) ) cos ( c t ) 2 (2 22 ) = 2 ( ) + 2 cos ( ) 2 ( ) + 2 2 2 ( ) = 1 2 + 2 4 + cos ( 2 ) + 1 2 + 2 4 cos ( 2 ) + 2 cos 2 + + cos 2 + 2 8 cos 2 + + cos 2 If we suppress carrier signal in modulated signal, meaning th at the 1 in the amplitude of electric field is removed we end up with ( ) = 0 2 cos cos 2 in which cos ( ) cos ( ) 2 (2 23 ) = 2 cos + + cos 2 2 = 2 4 + 2 4 cos ( 2 ) + 2 4 cos ( 2 ) + 2 8 cos 2 + + cos 2 From equation (2 22 ) and (2 23 ), we can see that there is a low frequency heating component at 2 which can stimulate phantom properly. PAGE 37 37 2.1.4 CW Stimulation Experiment Using RF amplitude modulation technique illustrat ed in previous section, we could generate the low frequency heating function at twice of modulation frequency In experiment, it is realized by RF Signal Generator (RFSG) card from National Instruments (N I), which could be programmed by LabVIEW software to generate AM m odulated RF signal. Based on the analysis about resonant t hermoacoustic imaging, if we apply ampl itude modulated RF signal on phantom with an appropriate modulating frequency, we can st and a chance for stimulating phantom at one of its resonant frequencies. Therefore we might be able to leverage this res onant phe nomenon to generate thermoacoustic signal with much lower input power. This will save the cost of thermoacoustic imaging system significantly. In e xperiment we tried to sweep the modulation frequency in low frequency band with fine resolution and also uti lized a RF a mplifier to boost the power of RF signal to 100 W for stimulating the spherical phantom. However no thermoacoustic signal was observed but EMI signal, which was believed to be generated by amplifie r nonlinearity. The experiment schematic is ill ustrated in Figure 2 4 below. Figure 2 4. Experiment Schematic for CW Excitation PAGE 38 38 The little sphere insi de the tank was a phantom, which was made of distilled water, salt and gelling agent with combinations of different concentr ation levels, such as 5% salt, 5% gelling agent and distilled water, and this was all wrap ped with a piece of rubber The salt was first dissolved in distilled water. Then gelling agent was added to make the solution become a kind of gel. After that, a fin ger of rubber glove was cut off and used to wrap the phantom. Transducer was positioned on the tank wall and at the same height with phantom. The transducer used in our experiment was the unfocused immersion transducer from Paname trics NDT, model number V3 03 SU C enter fr equency is around 1 MHz, and 6 dB bandwidth is from 0.687 MHz to 1.25 MHz. RF signal was input from the coaxial to waveguide adaptor at the bottom of tank which primarily produced TE01 mode electromagnetic field while the output signal fro m transducer wa s amplified and displayed on spectrum analyzer T ypical outp ut signal is illustrated in Figure 2 5 below. Figure 2 5. Output Signal on Spectrum Analyzer Photo courtesy of Xin Guan. In picture above, input signal was set to 150 kHz, and r emarkable output signals were located at 300 kHz and 600 kHz, namely at 2 and 4 sometimes even at PAGE 39 39 higher harmonic frequencies. We swept the modulating frequency from 50 kHz to 2 000 kHz with 10 kHz step size, and the 2 and 4 har monic components kept showing up with no significant changes. The most convincing results should come from comparative study We did the same frequency sweep with and without phantom in the tank, compared s. It is believed that this output signal was generated by direct coupling from electromagnetic wave p revailing inside the tank to metal structure inside transducer, and the RF signal on low order harmonic frequencies was generated by nonlinearity of RF am plifier rather than heating effects. Based on the c urrent experimental results, we can draw the following conclusions about the absence of the thermoacoustic signal The p eak powers of s s t ypical Photoacoustic experiment s and our CW stimulation experiment s are listed and compared in Table 2 1 below It is showed that our energy density of CW stimulation is far less than all of the other experiments, which make s SNR of thermoacoustic signal too small to be observable. In CW stimulation, heat conduction might also hurt the resonance through damping effects In addition, the surface tension of the r ubber wrapping may suppre ss the acoustic vibration of phantom under heating. Table 2 1. Peak Power Levels for different experiment setup PEAK POWER ( ) ENERGY DENSITY( / 3 ) Lihong Wang 20 0.45 Kruger 25 0.56 Photoacoustic Imaging ~100 ~6 CW Stimulation < 0.1 <0.0023 PAGE 40 40 Another big concern is related to safety issue. According to IEEE safety standard [31], und er a controlled environment, m aximum permission exposure (MPE) allows the power density up to 10 2 for human safety. The power density limit is relaxed to 20 2 for partial body exposure. Given dimension s of outlet for our microwave generator ( 72 34 ), maximum input power sh ould be 0.49 W. Considering reflection du e to mismatched impedance at waveguide bounda ry and microwave attenuation which soft tissue is experiencing, input power could be increased to 10 W or 20 W at most. This power level is far below the threshold power to trigger any observable thermoacoustic signal. C urrent thermoacoustic imaging systems are all based on pulsed microwave stimulation, and therefore there is a lack of mature image reconstruction algorithm for CW microwave stimul ation Based on all considerations above, we switched to high peak power pulsed microwave stimulation to construct our thermoacoustic imaging system. 2.2 Pulsed Microwave Excitation C omplete theoretical analysis and numerical simulation for thermoaco ustic signal generation is presented in [1], while main ideas are include d here for illustration purpose. Under thermal confinement condition, the inhomoge neous wave equation relating deposited microwave energy and generated thermoacoustic pressure wave is presented below. 2 1 2 2 2 = (2 24 ) where is th e pressure wave generated by heating is the sound speed in phantom, PAGE 41 41 is the thermal expansion coefficien t and is specific heat of phantom material per unit mass. The solution for this differential equation can be formulated as below. = 4 1 ( ) (2 25 ) The position of pressure measuring point is r while the source point for thermoacoustic signal is The integral is calculated inside the sphere centered at r with radius is cal culated as below. = It is difficult to get an analytic solution for equation (2 25 ), but for simple geometry, such as a thin slab and under delta heating, pressure wave generated could be formulated as below. = 2 2 ( ) (2 26 ) In which ( ) is defined as unity between the range 0 ( ) and zero elsewhere. is the thi ckness of slab. If heating function can be expressed temporally by ( ) pressure response from this specific heating functio n can be calculated based on convolution theory as below. = ( ) (2 2 7 ) T hermoacoustic pressure signal impinges on transducer, and generates piezoelectric signal observed on oscilloscope If the impulse response of the transducer can be expressed as ( ) the piezo electric signal generated can be formulated as below based on convolution theory. = ( ) (2 2 8 ) Numerical simulati on for thin slab phantom in [1] is listed below for illustration. Figure 2 6 (a) is the thermoacoustic signal from a n infinite thin slab under delta heating, and Figure 2 6 (b) is the temporal profile of microwave pulse. Figure 2 6 (c) is the PAGE 42 42 convolution r esult of (a) and (b). Figure 2 6 (d) is the impulse response of transducer. Figure 2 6 (e) is the convolution result of (c) and (d). Figure 2 6 Numerical simulation for pulsed microwave stimulation (Source : G. Ku and Med. Phys ., vol. 27, no. 5, pp. 1195 1202, May. 2000. ) PAGE 43 43 CHAPTER 3 IMAGING ALGORITHM AN D FPGA IMPLEMENTATIO N 3.1 Digital Beamforming B eamforming algorithm is a widely used algorithm in array signal processing, with applic ations ranging from Radar/Sonar systems, wireless communication and medical imaging. B eamformer is basically a spatia l filter which operates on output of an array of sensors/antennas in order to enhance the signal of interest (SOI) relative to background n oise and directional interference [32][33][34]. Conventional beamformer, such as Delay and Sum beamformer, is data independent and robust. However th e performance is inferior to adaptive beamformer, which can automatically adapt the sp atial filter coeffici ents to changing signal characteristics as illustrated by Fig ure 3 1. A daptive antenna array is showed below which can adaptively form the beampattern according to moving locations of targeted user and interfering user. Figu re 3 1. Adaptive Antenna Array Usually beamformi ng system is working between analog front end circuitry and digital baseband p rocessing as illustrated in Figure 3 2. Firstly, the analog signal is PAGE 44 44 picked up by ante nna array in RF application or transducer a rray in ultrasound application. Then the a nalog signal is processed by receiver circuitry which performs signal conditioning, such as filtering and amplification. The first stage is commonly a low noise amplifier (LNA) which can amplify the analog signal w ith as little increase of noise power as possible. Then an analog bandpass filter or lowpas s filter is used to wipe out out of band noise and interference signals. In RF application, the signal usually needs to be down converted to intermediate frequency ( IF) or to baseband frequency by analog mixer and local oscillator. After these preprocessing, the analog signal is converted to digital signal by ADC (Analog Digital Converter). This series of operation is p erformed on every channel in antenna array. D igit al signals collected from multiple channels are then fed into digital beamforming network, whic h tries to optimally combine multiple signal streams into one signal stream which has lower no ise and interference level. Beamforming output could then be passed down to digital baseband processing stage such as demodulation in wireless communication application. A daptive beamforming network consists of a spatial filter and an adaptive weight update module. The spatial filter actuall y forms the beam pointing to the desired signal dire ction and nullifies the strong interference signal using the calculated weight coefficients The system diagram is showed in Fig ure 3 3. The adaptive weight calculation module takes channel data information as input, such as desired signal direction, signal covariance matrix or reference signa l, calculates adaptive weight coefficients and feeds th em into the spatial filter. W eight coefficients are updated a s frequent as possible given hardware lim itations in order to capture fast chan ging signal characteristics. PAGE 45 45 Figure 3 2. Beamforming Network Figure 3 3. System Diagram of Adaptive Beamformer In our thermoacoustic imaging application, a modified Robust Capon Beamfor ming algorithm is adopted as ima ging algorithm [4][18], because it has supe rior performance over conventional delay and sum beamforming algorithm. The derivation and hardware implementation details of our imaging algorithm will be presented in following sections [ 52 ] [53] PAGE 46 46 3.2 Robust Capo n Beamforming Algorithm Adaptive Beamforming is a widely used technology in array signal processing. Its various applications include Radar and Sonar systems, Smart Antenna systems for wireless communications, and medical imaging methods. Capon Beamforming is one of the earliest adaptive beamforming techniques devised. It has the potential to provide significantly better resolution and interference reject ion capability compared to conventional delay and sum beamformer [32][33][34]. However, it is sensitive to model errors, such as array steering vector errors including array calibration errors and imprecise knowledge of d irection of arrival (DOA) of SOI. Many algorithms have been proposed to mitigate this problem, and one of them being the Robust Capon Beamf orming (RCB) [29]. This algorithm belongs to diagonal loading class ( algorithm s which adds a diagonal loading factor in covariance matrix to make it robust for beamforming) and its diagonal loading level can be calculated exactly according to the uncertai nty set of array steering vector. Moreover, its computational complex ity is comparable to that of Capon beamformer [29][34]. 3.2.1 Standard RCB f or Narrowband Signal 3.2.1.1 Capon Beamforming Conventional Be amforming methods such as delay and sum beamform ing have gained wide use in many applications including medical imaging. It is simple, robust and easy to implement in software or hardware. However its drawbacks include poor resolution and high side lobe levels when compared with other adaptive beamformi ng approaches. Capon Beamforming is capable of adaptively calculating the req uired weight vector based on sample cov ariance matrix obtained from array output data, which passes the SOI in an undistorted manner while placing deep nulls in directions of PAGE 47 47 inte rferers [32][33][34], resulting in strong interference suppression. The spatial filtering formulation of Capon Beamforming is summarized below (a detailed derivation of Capon Beamforming can be found in [29]). The data model used here is: ( n ) = 0 s 0 n + ( n ) (3 1) w here y n = y 1 n y 2 n y M ( n ) T is the received signal vector at array output, with n represent ing sampling time and M being the number of antenna elements. 0 is the true steering vector of SOI, s 0 ( n ) is the waveform of SOI, and the term n represents interference plus noise. Under the framework of Capon Bea mforming, the calculation of wei ght vector can be formulated as the linearly constrained optimization problem given by: min 0 = 1 (3 2) w here is the t heor etical covariance matrix of received signal and is the output power of Capon Beamformer. E quation (3 2) illustrates the optimization problem that the total output power is minimized while the SOI is passed undistorted, which indicates that noi se and interferences are minimized. The solution to this optimization problem leads to the optimal weight vector of Capon Beamformer: = 1 0 0 1 0 (3 3) In practical applications, is difficult to obtain, and thus is usually replaced by sample covariance matrix which is calculated from received signal vectors as : = 1 ( ) ( ) = 1 (3 4) where the integer N is the number of snapshots used in beamforming. Also, the true steering vector 0 may not be available in practice, and the assumed steering vector 0 is used instead in Equation (3 3). Unfortunately, a slight deviation between the assumed PAGE 48 48 steering vector and the true steering vector of SOI may cause a significant degradation of performance for Capon Be amformer, because in this case the beamformer treats SOI as an interfering signal. Therefore, numerous algorithms have been derived to improve the robustness of Capon beamforming [29][34]. 3.2.1.2 Robust Capon Beamforming The RCB algorithm implemented in o ur thermoacoustic imaging project was propos ed in [29], which deals with narrow band complex valued array signals. The main idea of the algorithm is presented in this section. The algorithm assumes that the true steering vector of SOI can be estimated with in an uncertainty set based around th e assumed steering vector of SOI. The formulation of RCB is summarized as: min = 1 0 (3 5) w here is a user defined parameter for the size of the uncertainty set of the true steering vector around an assumed steering vector 0 By using a Lagrange multiplier methodology, the constrained optimization problem of RCB ab ove can be solved in four steps: Step 1. Compute e igendecomposition (EVD) of = (3 6) In equation (3 6), represents ei gen vector matrix, while is a diagonal matrix containing eigen value s of as diagonal elements. Step 2. An itera tive method is used to solve Lagrange multiplier from equation (3 7). 2 1 + 2 = 1 = (3 7) where is the number of antenna elements, is the element of vector = 0 PAGE 49 49 is the eigen value of (al l eigen va lues are arranged in the increasing order along the diagonal elements of ), is the Lagrange multiplier. Step 3. After obtaining the Lagrange multiplier an estimate of the true steering vector for SOI can be calculated from equation (3 8). 0 = 1 + 1 0 (3 8) Step 4. Finally substituting 0 into equation (3 3) yields the optimal weight vector for RCB algorithm: = 1 0 0 1 0 (3 9) A detailed derivation of RCB algorithm can be found in [29]. 3.2.1.3 Numerical s imulati on Based on the introduced adaptive beamforming algorithm above, numerical simulation is performed in order to compare the performance of standard Capon Beamforming (SCB) and Robust Capon Beamforming (RCB). In t he simulation configuration, uniform linear antenna array with 10 an tenna elements was assumed. N oise power was set to 0 dB, and the SOI signal power was set to 20 dB at 10 degree ( 90 degree to 90 degree), while four strong interferers with 40 d B power were assumed at 75 degree, 60 degree, 10 de gree and 25 degree respectively. The assumed DOA of SOI was deliberately set at 11 degree with 1 degree deviation from the true DOA of SOI in order to test the robustness of RCB algorithm. A ntenna elements were positioned at half wavelength separation dis tance between each other. The waveforms of SOI signal, interferer signals and noise were generat ed by random sequence generator in order to make them irrelevant to each other. Both the standard Capon Beamforming and Robust Capon Beamforming were PAGE 50 50 applied to the array data, and tried to pick out the SOI signal. The sim ulation results are listed in following figures In Fi gure 3 4, the blue line is the total received signal at channel 1. The red line is the SOI signal expected to be received at channel 1, whil e the green line is the interferer signal at channel 1. The black line represents noise which has the same magnitude at all channels It can be found that the received signal is dominated by the interferer signal because of its high er power level compared to SOI signal In Figure 3 5, the filtered signal by SCB is compared agai nst SOI waveform. It can be noticed that the SOI is almost filtered out due to the 1 degree deviation between the assumed DOA of SOI and the true DOA of SOI, which vividly demonstrate s the la ck of robustness of SCB. In Figure 3 6, the filtered signal by RCB is compared a gainst SOI waveform. It can be verified that the SOI is successfully extracted from severely corrupted totally received signal under the 1 degree imprecision of knowled ge for DOA of SOI, which confirm s the robustness of RCB. In Figure 3 7, the beampatterns of SCB and RCB are plotted together which tells us why the SCB algorithm filters out the SOI signal while the RCB algorithm maintains the SOI signal It can be found that at 10 degree where the true SOI comes from, the SCB beampattern forms a deep null killing the SOI signal because it thought the true SOI signal comes from 11 degree and treats the SOI at 10 degree as an interferer signal, while the RCB beampattern fo rms a peak which passes through the SOI signal und istorted. In Figure 3 8, power estimates from SCB and RCB are plotted against the true SOI power. It can be found that the power estimate from SCB is far less than the true SOI power, while the power estima te from RCB is around the true SOI power with acceptable error. PAGE 51 51 Figure 3 4. Signal Comparison Figure 3 5. SCB Filtered Signal PAGE 52 52 Figure 3 6. RCB Filtered Signal Figure 3 7. Beam Pattern Plot PAGE 53 53 Figure 3 8. Power Estimates Plot 3.2.2 Time Delayed RCB f or Wideband Signal Robust Capon Beamforming algorithm prese nted above mainly deals with narrowband complex valued input dat a. In narrowband situation, the steering vector is a complex valued vector representing the phase a nd amplitude relationship of r eceived signals among multiple antennas in the array. However, in our thermoacoustic imaging project, the signal picked up by transducer is a broadband (around 1 MHz at baseband) and real signal, which is not able to be processed directly by the original R obust Capon Beamforming algorithm. Another unique issue in our medi cal imaging application is that the ultrasonic signals received by different transducer channels are experiencing different path losses and need to be compensated before performing beamfor ming. Therefore, the original Robust Capon Beamforming algorithm was modified in [4][18], and applied in our thermoacoustic imaging project. PAGE 54 54 Figure 3 9. Illustration of time delayed RCB In the modified time delaye d RCB presen ted in [4][18], acoustic signals picked up by transducers firstly are compensated in amplitude according to different path loss at different channel Then the compen sated signals are aligned in time domain based on different delay on every channel, which i s determined by the distance between the imaging point and the transducer. As a result, the steering vector for this specific imaging point ha s become an all one vector. S ubsequent steps follow the standard RCB algorithm with only real valued data in all c omputations. After the adaptive weight vector for this imaging point has be en calculated, waveforms from multiple tr ansducers are combined using beamforming coefficients just obtained. E nergy or peak to peak value of the combined waveform is considered to be the image intensity at this specific imaging point The operation of time de layed RCB is illustrated in Figure 3 9, in which a typical imaging point is showed. In real operation, the entire imaging area is scanned PAGE 55 55 point by point with acceptable resoluti on. In doing so, we can get an image of the intensity of acoustic sources inside the imaging area, which corresponds to the amount of energy absorbed during electromagnetic illumination and directly re lates to the conductivity of phantom at that location. 3.2.3 Computational Complexity f or Different Algorithm Implementation M ajor computational steps for RCB are summarized as ma trix multiplication (include calculation of covariance matrix and estimated steering vector, which will be explained in follow ing sections), eigen value decomposition for symmetric matrix, iterative method for solving Lagrange multiplier and FIR filtering. The most computational intensive tasks are matrix multiplication and eigen value decomposition for symmetric matrix, which ha ve computational complexity of ( 3 ) [29]. The well known numerical algorithm for eigen value decomposition is the Cyclic Jacobi Method which is simple, robust and amenable to fixed point implementation. The computational complexity of Cyclic Jacobi M ethod is analyzed and presented in [51], which gives also ( 3 ) complexity for every sweep in calculation. In conceptual Microprocessor or DSP processor model, there is only one multiplier and one adder, or more concisely, one MAC operational unit available for computation. Although in modern Microprocessor or DSP processor, there are multiple integer and floating point functional units available, we can first assume the conceptual model intr oduced above for simplicity In our FPGA implementation we realized linear computational unit array with dimensions to exploit parallelism in matrix multiplication and eigen value decomp osition in order to speed up computation from the architectural perspective. The schematic for matrix multiplication is listed in Figure 3 10. PAGE 56 56 Figure 3 10. Schematic for Computational Complexity Analysis for Matrix Multiplication The l eft pa rt of Figure 3 10 represents Microprocessor or DSP processor implementation, which utilize th e single MAC unit to perform computations for every single elem ent in result matrix yielding ( 3 ) computational complexity. The r ight pa rt of Figure 3 10 represents FPGA implementation, which utilize the line ar MAC unit array to perform computations for every element in single column from result matrix yielding ( 2 ) computational compl exity. The same situation occurs in eigen value decomposition. In Microprocessor or DSP processor implementation, the only computational unit is used to perform computations for every row e lement in Cyclic Jacobi Method which results in ( 3 ) comput ational compl exity. On the other hand, in FPGA implementation, the linear computational u nit array is used to perform computations for entire row at the sa me time, which will speed up computation by times and yield ( 2 ) computational complexity. T his wi ll be further illustrated in following sections. C omparison of run time performance between software implementation and FPGA implementation are presented as follow. The application is an 8 channel PAGE 57 57 wideband Robust Cap on Beamforming algorithm. T empor ary delay and amplitude compensation have been done before testing. As a result, only the computation time is compared. For software implementation, the platform is a Pentium(R) dual core CPU running at 2GHz in a laptop with 3GB memory, while the platform for FPGA implementation is a Xi linx SX 35 platform device. C omputation time scaling with FPGA operating frequency is showed in Figure 3 11. From the chart, we can notice that from the frequency larger than 35 MHz, FPGA implementation shou ld run faster than software implementation. C urrent runn ing frequency is 80 MHz, which gives FPGA implementation around 2.3 times speedup over software implementation compared with the expected 8 times speedup. It should be noticed that only computation time is considere d i n comparison while data transfer time from PC to FPGA board is not taken into consideration, which could slow down overall FPGA computation further. Figure 3 11. Run time Performance Comparison 0 200 400 600 800 1000 1200 10 20 30 40 50 60 70 80 90 100 Computation Time (us) FPGA Operating Frequency (MHz) Run time Performance Comparison Software Computation Time FPGA Computation Time PAGE 58 58 3.3 FPGA Implementation of RCB In previous sections, the pr inciple and computational steps of the RCB algorithm are illustrated in details. Th is section mainly deals with hardware implementation of the algorithm. For any digital signal processing algorithm, there might be multiple alternative implementations which have different pros and cons. A digital signal processor based implementation is the most flexible and cheap approach to adopt however it may not provide the best performance. An ASIC (Application Specific Integrated Circuit) spe cially designed for runni ng beamforming algorithm usually is too expensive to design changes any more. A FPGA (Field Programmable Gate Arrays) provides a good platform for implementing digital s ignal processing algorithm due to its high processing power and gr eat flexibility when compared with its DSP processor and ASIC counterparts. It is also an ideal platform for algorithm design and verification. Therefore, in our th ermoacoustic imaging proje ct, FPGA base d hardware implementation of RCB algorithm is investigated and realized on Xilinx Virtex4 SX35 platform. As illustrated in previous section, compared with Microprocessor or D SP processor implementation, FPG A implementation can exploit paralle lism in computation more extensively by utilizing the linear array of computational unit to get a n fold s peed up (n is the dimension of matrix). In reality, this speed up can be compromised by two factors Fi rst, modern Microprocessor or DSP processor als o have multiple computational units such as in Superscalar architect ure, but the architecture of computational unit pool is for general purpose usage, which can not be specialized for any specific application. Second, the operating frequency for Microproce ssor or DSP p rocessor is much higher than FPGA implementation. As a result, FPGA PAGE 59 59 implement ation may not be faster than Microprocessor or DSP processor implementation in absolute time. However, the research and development of FPGA implementation for RCB is still meaningful in that FPGA implementation of RCB could work as the hardware a ccelerator, which could save Micropro cessor or DSP processor from computational intensive task and let them focus on other essential jobs, such as running O perating System, con trolling display. FPGA implementation could also work as the verification platform for ASIC module implementati on, which might be included in SOC design. Due to the fast developmen t of Smart Antenna Systems, many researc h effort s have been devoted to hard ware implementation of adaptive beamformer. Because of its strong interference rejection capability and simplicity, the Capon beamformer, also known as MVDR (Minimum Variance Distortionless Response) beamformer, has been the target for hardware implementat ion in several works. Dick and Harris performed a real time QR decomposition based Capon Beamformer implementation [35]. They provided a flexible and compact architecture for a real time beamformer, together with their Matlab/Simulink design flow for heter ogeneous architecture designs. Edman extensively investigated many alternative realizations for the Capon Beamformer and use d their FPGA implementation in Channel Sounder application [36]. Liu et al, presented the Capon Beamformer implementation based on d ichotomous coordinate descent iterations [41], instead of the popular QR Decomposition method, resulting in an area efficient realization. In addition, there exist many work s on implementing an MVDR beamformer with RLS (Recursive Least Squares) adaptive al gorithm [37]. PAGE 60 60 As previously mentioned, Capon Beamformer suffers from the lack of robustness. Due to this drawback any model errors could cause significant performance degradation, which limits its practical applications. Wang et al, presented an implement Beamformer [38], which has better performance than Capon Beamfor mer. It has been proven that the Robust Capon Beamformer implemented in our thermoacoustic imaging project is equivalent to Ger Beamformer [29], i.e., both algorithms produce the same weight vector. However the approach we adopted gives a simple way to eliminate the scaling ambiguity when estimating the SOI power and is also more computationa lly efficient [29]. Besides, har Beamformer is not mentioned in [38]. On the other hand, major building blocks and architectures of our RCB implementation are illustrated in this thesis. 3.3.1 Implementation of Standard RCB 3.3.1.1 System block d iagram Figure 3 12 System Block Diagram PAGE 61 61 The system block diagram for hardware implementation of the Standard Narrowban d RCB algorithm is shown in Figure 3 1 2 A rray output data is fed into the Block Ram Array, and every two Block Ram store the real part and imaginary part of one channel of d ata respectively. D ata is initially read from Block Ram Array into the Matrix Multiplication Modu le which computes sample covariance m atrix. This matrix is used as input for Eige ndecomposition Module. After e igen values and eigen vectors are computed the Lagrange Multiplier Solver module uses these terms to calculate the Lagrange multiplier based on an iterative method. The Lagrange multiplier i s then utilized to determine adaptive weight vector. This weight v ector generated is then applied to a rray output data to complete beamforming process. 3.3.1.2 Complex Multiply and Accumulate Unit The pipelined architecture is adopted for realizing multiplication and accum ulation of complex numbers. Guard bits are added in intermediate stages and accumulator stage to prevent overflow, while the final result is truncated to 18 bit fixed point representation giving sufficient accuracy. S everal such modules are instantiated in a linear array fas hion in order to parallelize t he c omputation. The datapath of M AC module is illustrated in Figure 3 1 3 as below. Figure 3 13 Datapath of Complex MAC Module PAGE 62 62 3.3.1.3 Extended Cyclic Jacobi Method Cyclic Jacobi Method is a well know numerical algorithm for eigenvalue decomposition of symmetric real valued m atrix [39][44]. To deal with complex valued data an extension of the Cyclic Jacobi Method is adopted to compute Eigen decomposition of a Hermitian Matrix [42] The for mulation is outlined as follow. 0 = + 1 = (3 10) A is the original matrix, and A k + 1 i s the matrix after k th iteration. J pq is transforma tion matrix which eliminates off diagonal element in the p th row and q th column of A k J pq can be decomposed into two matrix inner matrix M pq and outer matrix R pq as below: pq = pq M pq (3 11) In which M pq has the form: = (3 12) = = = 1 2 [ ] While pq has the form: = (3 12 ) = = cos = = sin = 1 2 tan 1 2 where i j is the Kronecker Delta and k p q is the complex valu ed off diagonal element in the p th row and q th column of k while arg [ k p q ] and k ( p q ) are the angle and norm of k p q respectively. PAGE 63 63 After multiplying both sides of k by inner matrix pq k p q and k q p wi ll become real numbers. Then multiplication of outer matrix pq reduces these two off diagonal elements to zero by plane rotation. During every transformation, two rows of elements of k are updated with diagonal elements being k p p and k ( q q ) and off diagonal elements being k ( p i ) and k q i ( i p q ) Both k ( p q ) and k q p are reduced to zero, as illustrated. + 1 = + 2 + sin ( 2 ) cos ( 2 ) A pp k A qq k 2 (3 13) + 1 = + 2 sin ( 2 ) + cos ( 2 ) 2 + 1 = + 1 = 0 + 1 = cos ( ) sin ( ) + 1 = sin ( ) + cos ( ) A detailed derivation of the extended Cyclic Jacobi Method can be found in [42]. Datapath of EVD Module consists of a linear array of Processor Elements (PE), and is illustrated in Figure 3 1 4 Figure 3 14 Datapath for Eigenvalue Decomposition Module PAGE 64 64 CORDIC processor [40] takes k p q fr om k calculates th e trigonometric functions of angle pq rotation angle pq and norm of k p q and feeds them to diagonal PE and off diagonal PEs. The diagonal PE and off diagonal PEs work in parallel to spee d up computation. The whole d atapath operates on the entire row of A k at once. Datapath for C ORDIC processor is shown in Figure 3 1 5 It contains two identical CORDIC units in order to evaluate the angles and calculate the trigonometric functions of pq 2 pq and pq in an e fficient and pa rallel way. Datapath for di agonal PE is illustrated in Figure 3 1 6 while datapath of off diagonal PE is not included for brevity. D atapath for calculating eigen vectors is almos t the same as that shown in Figure 3 1 4 but only with off diago nal PEs, and takes the identity m atrix as input and generates eigenvectors as output [44 ]. Based on experimental result a total of 3 sweeps achieves sufficient accuracy for eigenvalues on a 4 by 4 complex valued matrix. Due to the progressively small ampl itude of off diagonal elements, the evaluation of pq will become worse leading to poor accuracy of computations for eigen vectors. However the whole algorithm is not very sensitive to the inaccuracy of eigenvectors, and still finds the adaptive weight vector with an acceptable accuracy Figure 3 15 Datapath for CORDIC Processor PAGE 65 65 Figure 3 16 Datapath for Diagonal Processing Element 3.3.1.4 Lagrange Multiplier Solver The Lagrange Multiplier is obtained by solving e quation (3 7) iteratively. Many iterative root finding algorithms can be applied here and we adopt the binary search method due to its simplicity and robustness. The number of iterations i s a user defined parameter. D atapath of the Lagrange multipli er solv er is illustrated in Figure 3 1 7 Figure 3 1 7 Datapath of Lagrange Multiplier Solver The divider is based on CORDIC processor [40], and the mult iplier is from Xilinx primitive This datapath is timely shared by computations for M vecto r components. The PAGE 66 66 results are added together by an add er tree. The sum is taken by controller and compared to in equation (3 7). If the sum is larger than then we take the upper half region as the solution region due to the monotonically decreasing property of the left side of equation (3 7) and perform the next iteration of binary search on this chosen reg ion. Otherwise, we take the lower half region as the solution region and perform the next iteration of binary search on this chosen region. After a few number of iterations the Lagrange multiplier is found with an acceptable accuracy 3.3.1.5 Adaptive Weig ht Calculation From equation (3 9), we can se e that in order to calculate adaptive weight vector, we need to compute a 0 and R 1 a 0 Developing individual datapaths for each computation is possible, but consumes more hardware resources. We therefore reformulate the computations for a 0 and R 1 a 0 to make resource sharing possible, as shown below: 0 = 1 + 1 0 (3 14) = 1 + 1 0 = 1 1 + 1 0 0 0 0 2 2 + 1 0 0 0 0 0 0 0 0 + 1 0 = 1 0 1 0 = 1 1 + 1 0 (3 15) = 1 + 1 0 PAGE 67 67 = + 1 0 = 1 + 1 0 0 0 0 2 + 1 0 0 0 0 0 0 0 0 + 1 0 = 2 0 The matrices 1 and 2 have the same form, both being the result of a multiplication of a diagonal matrix and which si mply results in a scaling of rows of The matrices 1 and 2 are then multiplied by matrix As a result, these computations can be performed by utilizing the Matrix Multiplicat ion Module used in computing sample covariance ma trix with additional control This hardware resource sharing can save a s ignificant amount of area in FPGA with the penalty of longer latency for computation. 3.3.1.6 FPGA i mplementation The hardware implementation of the standard narrowband RCB algorithm is realize d using a Xilinx Virtex 4 SX35 FPGA. In our realization, an 18 bit fixed point representation is adopted which prov ides sufficient accuracy for present application. The array output consists of 4 channels of complex valued data, and the number of snapshots used for beamforming is 32. The resource utilization for our application is listed in T able 3 1. Table 3 1. Resource Utilization for RCB Implementation Resource Used Available Utilization Slices 14978 15360 97% FF 13420 30720 43% LUTs 23246 30720 75% B RAM 41 192 21% DSP48 192 192 100% PAGE 68 68 The latency for one sweep in EVD computation for a 4 by 4 complex valued matrix is around 540 cycles. With a total of 3 sw eeps, the overall latency of RCB algorithm is approximately 5485 cycles, which can be further impr oved by exploiting more parallelism and reducing control overhead. Since this design is running at a clock speed of 100MHz, roughly 54.85 microseconds delay i s expected for one update of adaptive weights. Compromises have been mad e in choosing between t ria ngular systolic array and linear array processor architecture for EVD core in order to save hardware resources at the expense of increased latency and control complexity. 3.3.2 Implementation of Wideband RCB In wideband R CB algorithm implementation, the input sign al is the real valued data, and the signal processing operations involved are all real valued computa tions. The signal alignment and amplitude compensation are performed on the host PC, and the time delayed, compensated data is then transferred t o the FPGA by the PCI interface and processed by the modif ied wideband RCB algorithm. Major differences between the original narrowband RCB algorithm and the modified wideband RCB algorithm lie in the eigenv alue decomposition step, although there are other minor differences in other computational steps In the standard RCB, the covariance matrix is a comple x hermitian matrix, while in the modified wideband RCB, the covariance matrix is a real sy mmetric matrix. As a result, the eigenv alue decomposition modul e in the wideband RCB implementation is simply based on the cyclic Jacobi method. In the wideband RCB implementation, the multiply and accumulate operation is based on the real number, and th erefore is much simpler than MAC operation in the narrowband RCB implementation. PAGE 69 69 3.3.2.1 Cyclic Jacobi Method for Eigendecomposition Cyclic Jacobi Method is a well know numerical algorithm used to calculate eigen value decomposition of symmetric real valued matrix [39][44]. It is robust and the convergence is guaranteed The formulation is outlined as follows: 0 = + 1 = pq (3 16) w here is the original matrix, and k + 1 i s the resulting matrix after k th iteration. pq is the plane rota tion matrix which eliminates off di agonal element in the p th row and q th column of k thus has the form in equation (3 17). = (3 17) = = cos = = sin = 1 2 tan 1 2 ( ) ( ) w here i j is the Kronecker Delta and k p q is off diagonal element in the p th row a nd q th column of k After multiplying both sides of k by rotation matrix pq two rows of elements of k are updated with diagonal elements k p p and k ( q q ) and off diagonal elements k ( p i ) and k q i ( i p q ) Bo th k ( p q ) and k q p are reduced to zero, as illustrated below. + 1 = + 2 + sin ( 2 ) cos ( 2 ) 2 (3 18) + 1 = + 2 sin ( 2 ) + cos ( 2 ) 2 + 1 = + 1 = 0 + 1 = cos ( ) sin ( ) PAGE 70 70 + 1 = sin ( ) + cos ( ) A detailed derivation of Cyclic Jacobi M ethod can be found in [39]. Datapath of EVD Module consists of a linear array of Processor Elements (PE), and is illustrated in Figure 3 1 8 Figure 3 1 8 Datapath for Eigenvalue Decomposition Module CORDIC processor [40] takes k p q from k calculates the trigonometric functions o f rotation angle pq and feeds them to diagonal PE and off diagonal PEs. The diagonal PE and o ff diagonal PEs work in parallel to speed up computation. Th e whole datapath operates on entire row of k at once. Datapath for C ORDIC processor is shown in Figure 3 1 9 It contains two identical CORDIC units in order to calculate the trigonometric func tions of 2 pq and pq in parallel. Datapath for di agonal PE is illustrated in Figure 3 20 while the datapath o f off di agonal PE is illustrated in Figure 3 2 1 D atapath for calculating eigen vectors is nearl y the same as that shown in Figure 3 1 8 but o nly with off diagonal PEs, and takes the identity m atrix as input and generates eigen vectors as output [44 ]. Based on experimental result a total of 4 sweeps achi eves sufficient accuracy for eigen values and eigen vectors. PAGE 71 71 Figure 3 1 9 Datapath for CORDIC Processor Figure 3 20 Datapath for Diagonal Processing Element Figure 3 2 1 Datapath of Off Diagonal Processing Element PAGE 72 72 3.3.2.2 FPGA implementation of w ideband RCB Th e hardware implementation o f wideband RCB algorithm is realized using a Xilinx Virtex 4 SX35 FPGA. The algorithm simulation is performed using MATLAB, and the implementation is accomplished by VHDL programming. In our realization, an 18 bit fixed point representation is adopted whic h prov ides sufficient accuracy for present application. The array output consists of 8 channels of real valued data, and the number of snapshots used for beamforming is 32. The resource utilization for our application is listed in T able 3 2. Table 3 2. Reso urce Utilization for Wideband RCB Implementation Resource Used Available Utilization Slices 10674 15360 69% FF 9377 30720 30% LUTs 17042 30720 55% BRAM 24 192 12% DSP48 83 192 43% The latency for one sweep in EVD computation for an 8 by 8 real value d matrix is around 2312 cycles, which can be further reduced by eliminating more control overhead. With a total of 4 sw eeps, the overall latency of RCB algorithm is approximately 11374 cycles, which can be further improved by exploiting more parallelism. S ince this design is running at a clock speed of 80MHz, roughly 140 microseconds delay i s expected for one update of adaptive weight. Compromises have b een made in choosing between t riangular systolic array and linear array processor architecture for EVD co re in order to save hardware resources at the expense of increased latency and control complexity. PAGE 73 73 CHAPTER 4 EXPERIMENTAL RESULTS OF THERMOACOUSTIC IM AGING SYSTEM 4.1 Experimental Step up of TAI System 4.1.1 Introduction The experiment al setup is built to conduct proof of principle research for our thermoacoustic ima ging system. It can generate high peak power, short duration microw ave pulse used to illuminate phantom body, pick up the emit ted thermoacoustic signal by transducer subsystem, digit ize and record the signal by Data Acquisition subsystem (DAQ). The system block diagram is showed in Figure 4 1. The whole system can be divided into three subsystems, including high power subsystem, transducer array subsystem and data acquisition subsystem, which are described in details as follow Figure 4 1. TAI System Block Diagram PAGE 74 74 4.1.1.1 High Power Subsystem High power su bsystem operates to generate high peak power m icrowave pulse, feed it into e xperimental tank and provide reflection protection for micro wave source. M ajor components involved are pulsed microwave source, pulse generator, circulator, dummy load, H bend and supporting table. T he real system is showed in Figure 4 2 below. Figure 4 2. High Power Subsystem Photo Courtesy of Xin Guan. HP 214 A pulse generator [45] is used to generat e the trigger signal for pulsed microwave source. The output signal should be a pulse with peak voltage 200 V. The pulse repetition rate could be adjusted according to specific application requirement in the range o f 30 Hz to 1100 Hz, or even manual mode. It can operate in several trigger PAGE 75 75 modes, including manual trigger, exterior trigger and internal trigger. In our experiments, we set it to internal trigger mode and 100 Hz pulse repetition rate. The control panel in o ur configuration looks like Figure 4 3 below. Figure 4 3. Control Panel of Pulse Generator Photo Courtesy of Xin Guan. The output signal from pulse generator is also used as the t rigger signal to synchronize data acquisition (a certain amount of de lay is inserted) after firing the electromagnetic illuminatio n. Because the amplitude of trigger signal is 200 V, a voltage attenuator with 20 dB attenuation at least is used to reduce the voltage to 2 V and feed into the external trigger port on digitizer card. The pulsed microwave source is from Radio Research Instrument Co., and the par t number is 12 1 21MOD [45]. 3 phase electricity with 120 V nominal voltage is used to power the equipment. If operate d at its maximum capability, this microwave source is capable of supplying 0.75 us pulse at repetit ion rate of 1100 Hz and peak power of 120 kW. The fundamental frequency of microwave pulse is fixed around 3GHz. P ulse repetition rate can be adjusted by setting the HP pulse generator to an appropriate valu e like 100 Hz. The output power could be controlled by the high voltage applied inside the system, which can be changed by the knob on the front control panel. P ulse PAGE 76 76 duration is fixed to 0.75 us, which can only be changed by replacing the pulse forming net work circuit inside the equipment. The operatio nal status is summarized in Figure 4 4 below, in which the red part indicates the nominal operation condition for our experiments. Figure 4 4. Summary of Operational Conditions In real operation, we set the high voltage approximately to 2.0 kV, thus generating peak power less than 25 kW in order to keep radiation level as low as necessary based on safety concern The columns named Typical Magnetron Current and Typical Rectifier C urrent are two indicators of o utput power level on the front control panel. Due to the low power op erational condition, both of the indicators are approaching zero. PAGE 77 77 Because of mismatch of impedance at the pow er feeding interface between experimental tank and H bend, part of the input power will be reflected back toward power source, which is harmful for the equipment. The ci rculator is inserted between H bend and power source, and serves as reflection protection device. The reflected power will be directed to the dummy load and absorbed locally. 4.1.1.2 Transducer a rra y s ubsystem The transducer in usage is the unfocused immersion transducer from Panametrics NDT, model number V303 SU, which is showed in Figure 4 5 below. The center frequency is around 1 MHz, and the 6 dB bandwidth is 0.6 MHz. The impulse response lasts for nearly 4 us [46]. Figure 4 5. Transducer Photo Courtesy of Xin Guan. In order to form a transducer array, a motor driven rotational stage is used together with one transducer to simulat e a circular transducer array by rotating the single transducer to different c ircular positions during data acquisition process. This method is simple and economically efficient, but needs more time to collect data PAGE 78 78 compared with real transdu cer array. The schematic is in Figure 4 6. The front view is on the left while the top view is on the right. The orange part is the step motor fixed on the top of the aluminum fielding box. The mechanical arm punches through the motor, connects the transducer and serves as the rotational sta ge. The blue part represents experimental tank f illed with mineral oil. P hantom sits in the middle of the tan k, which stays on the top of H bend. Togeth er with Figure 4 1, we can find that the relative position ing between phantom and t ransducer trajectory is not fixed during different rounds of experiment. The motor and transducer system is at tached to the sh ielding box, while the expe rimental tank and phantom is fixed to H bend. Every time we take off the sh ielding box, change phantom and put back the sh ielding box, we can not exactly recover the previous positioning of the sh ielding box. That is one of major drawback s for current transducer array subsystem The direct consequence of this alignment issue is that the position of phantom in our final images is not consistent, although the shape is the same. Figure 4 6. Motor and Transducer System PAGE 79 79 In our experiments, t he transducer is swept along a circle with resolution of 3 degree, meaning total 120 scan ning positions. At every circular position, high power system lau nches multiple microwave pulses and generates the same number of thermoacoustic respo nses, which are picked up by transduce r and averaged and stored by DAQ system. After sc anning the current position, motor rotates 3 degree to the next position and performs data collection at that position. In doing so, the simulated transducer array can collect data as if there were 120 transducers collecting data at the same time. The scan ning schematic is s howed in Figure 4 7 below. Figure 4 7. Scanning Schematic 4 .1.1.3 DAQ s ubsystem DAQ subsystem performs two functions, signal conditi oning and data acquisition. T hermoacoustic signal is deeply embedded in noisy background inclu ding amplifier noise (especially with large amplification gain like 60 dB) and strong interference signal s PAGE 80 80 such as EMI (Electromagnetic Induced) signal. F irst, the signal coming from transducer needs to be amplified and filtered, which is accomplished by t he pulser/receiver equipment. Then the outcome has to be digitized by digitizer card, and average d multiple times to increase SNR (Sign al to Noise Ratio). At last, averag ed signal can be recorded in data file, which will be utilized subsequently to perform imaging. The schematic of DAQ subsystem is showed in Figure 4 8 below. Figure 4 8. DAQ Subsystem Schematic The pulser/receiver in usage is Olympus Panametrics NDT 5800. The module operates in external pulser mode in orde r to solely utilize receiver module. The high pass filter cutoff frequency is set to 1 kHz, and the low pass filter cutoff frequency is set to 10 MHz. The amp lifier gain is set to 60 dB, maximum allowable value, to amplify thermoacoustic signal as large as po ss ible [47]. The operating parameters are displayed in Figure 4 9 below. Figure 4 9. Pulser/ Receiver Operating Parameters Photo Courtesy of Xin Guan. PAGE 81 81 The digitizer card in usage is a NI P XI 5122 digitizer hosted in NI PXI box. It can work from 1 MSample/ s to 100 MSample/s with 14 bit resolution on two ADC channels [48]. There is a software control panel embedded in the VI, which can be used to instantiate the digitizer in LabVIEW program. In the panel, you can set the operating parameters for applicat ion. The general setting for digitizer card in our experiments can be found in Figure 4 10 below. Figure 4 10. Configuration for NI 5122 Digitizer In the configuration above, channel 0 has been selected for data acquisition. Input range is set to 2 V, no off set voltage. DC coupling is chosen to directly sample the entire waveform. High input impedance is set to capture the most part of signal. In this particular setting, 20 MSample/s sampling rate is taken with record length 2000 samples, corresponding to 10 us in time and 15 mm in distance, which covers the entire phantom area. In some experiments, if the phantom body is large, you need to set record length large enough to make sure you collect the generated thermoacoustic PAGE 82 82 signal from all over the phantom bod y. The next important configuration is the trigger se tting which can be found in Figure 4 11 below. Figure 4 11. Trigger Setting for NI 5122 Digitizer I n the trigger setting above, external trigger source i s selected, which comes from HP pulse generator after certain amount of attenuation. The trigger type is set to edge trigger and the slope is set t o positive, which means that data acqui sition will start as soon as input signal level exceeds the threshold level (10 mV in this particular configuration) for the first time. Another important point is the delay specified as 3 us. The reason why we introduced a delay before the start of data acquisition is that at the moment of electromagnetic illumination, there wil l be a huge EMI signal from transducer whi ch generally will excee d the voltage range we set in previous section. If we increase the voltage range to accommodate this huge spike, we will lose resolution in subsequent measurement. As a result, a 3 us delay period is introduced to gate out this huge spike. PAGE 83 83 A LabVIEW program was developed to coordinate between digitizer, motor control module, average functionality and data file generation. The control panel o f the LabVIEW program is in Figure 4 12 below. Figure 4 12. Control Panel The angle indicat or tells the current position of transducer, while the summation counter shows current averaging ti mes. The top window displays sampled signal in real time, w hile the bottom window shows current averaged result. The data file will be generated for every tr ansducer position, and can be loaded into matlab for further processing. 4.2 Thermoacoustic Signal Generation Experiment In previous section, experimental setup of our TAI system is i ntroduced, which consists of high power subsystem, transducer array s ubs ystem and data acquisition subsystem. Using the system described earlier, we did the initial experiments to confirm the feasibility of thermoacoustic signal generation. In this initial experiment, we used the PAGE 84 84 old cylindrical tank with a bulk of phantom. Th e transducer was hung from the top surface of the tank in order to have better reception of acoustic signal. The important issue was that mineral oil should be used to fill the tank as bath material instead of distilled water [8], which has larger attenuat ion for microwave radiation. The experimental tank together with phantom and transducer is showed in Figure 4 13. Figure 4 13. Tank without phantom (left), Tank with phantom (right ) Photo Court esy of Xin Guan. C omparative studies have confirm ed the source of acoustic signal, one experiment without phantom and the o ther experiment with phantom as depicted by Figure 4 13. The relative positioning between transduc er and phantom is showed in Figure 4 14. Figure 4 14. Relative Positions between Transducer and Phantom PAGE 85 85 Distance is converted to time of flight by sound speed 1.5 mm/us [49]. The signal gen erated by the top surface of phantom should arrive at transducer around 3 2.6 us, and the signal genera ted by the bottom surface of phantom should arrive at transducer around 56.6 us. This signal should be corrupted by the signal gene rated by the table surface. A coustic signal generated i n waveguide should arrive at transducer a round 80 us. Due to the inaccuracy in measurements, estimates above are allowed to have several us deviation. L ow frequen cy signal generated inside phan tom body was filtered out by transducer. T hermoacoustic signal generated in this configuratio n is showed in Figure 4 1 5, which is the summation of recorded acoustic signal 16000 times instead of averaging Figure 4 15. Thermoacoustic Signal in the First Experiment In this initial experiment, the delay in data acquisition has not been inserted, and therefor e the EMI signal appeared at the beginning of data sequence. The signal from the top surface of phantom was located around 35 us, while the signal from the bottom surface of phantom was located around 58.6 us, which are in good agree ment with PAGE 86 86 previous time of flight analysis. To be specific, the sig nal appeared around 55 us in experiment without phantom is the thermoacoustic signal generated by the supporting table, while the sig nal appeared around 55 us in experiment with phantom is the thermoacoustic sign al generated by the bottom surface of phantom loaded with the supporting table. Th e signal generated inside waveguide arrived near 82 us as it was supposed t o be [50]. A zoom in view of above picture is showed in Figure 4 16, which mainly focuses on signal generated by the top and bottom surface of phantom. Figure 4 16. Zoom in View of the Thermoacoustic Signal It can be seen clearly that the signal existed when phantom was in the tan k, while it disappeared when phantom was outside the tank, which confirm ed that it was the thermoacoustic signal generated by the top surface of phantom. In order to prove the repeatability of the above experimental result, two more experime nts have been conducted with same operating conditions. Pictures of thermoacoustic sig n al generated are showed in Figure 4 17 and Figure 4 18 respectively. PAGE 87 87 Figure 4 17. Thermoacoustic Signal in the Second Experiment Figure 4 18. Thermoacoustic Signal in the Third Experiment As can be seen from pictures above, thermoacoustic signal genera ted from the top surface of phantom repeatedly showed up, which further confirmed that the signal studied was the thermo acoustic signal generated by phantom. 4.3 Thermoacoustic Imaging Experiment The initial experimental setup h as several drawbacks. Firs t, transducer was hung from the top of the tank, which can only be used to receive acoustic signa l in one PAGE 88 88 direction, while in real imaging operation a transduce r array is needed to collect signal from different directions. Second, the old cylinder tank is too small to allow t he rotation of mechanical arm used in transducer array subsystem described in Chapter 3. We need a larg er tank to make room for all circular positions transducer will stay. Thi rd, the phantom is a bulk of material which has higher condu ctivity. As a result, only the thermoacoustic signal generated at the boundary can be detected and recorde d. We need to build new phantom with bulk of background material whic h has lower conductivity and real target embedded inside background material whic h has higher conductivity in order to perform imaging. Based on these considerations, a new experimental tank has been made to accommodate phantom and transducer array, which is showed in Figure 4 19. Figure 4 19. New Experimental Tank Photo Courtesy of Xin Guan. The tan k can be fixed on the top of H bend, while phantom can be put at the center of the tank. New set of phantoms have been made to perform comparative studies. The background material is a cylinder with 30 mm in diameter and 40 mm in height, which consists of only distilled water and gelling agen t. The targets are made into ver y slim cylinders embedded in background with 2 mm, 4 mm and 6 mm in PAGE 89 89 diameter and same height as background, consisting of 15% salt and gelling agent. Different sizes of targets are used to test the spatial re solution of our imaging system. The schematic of new phantom is showed in Figure 4 20. Figure 4 20. Schematic of the phantom T hermoacoustic s ignals received at three circular positions dur ing experiment such as 0 degree, 120 degree and 2 40 degree, are displayed in Figure 4 21. Figure 4 21. Thermoacoustic signals at 0, 120, 240 degree PAGE 90 90 4.4 Imaging Results In order to perform comparison, different imaging algorithms, such as RCB, Standard Capon Beamforming (SCB) and Delay a nd Sum (DAS), are applied to thermoacoustic data collected from different phantom samples. The first set of images is constructed using data from phantom with 4 mm diameter, 15% salt target material. Figure 4 22. Imagi ng Result with 4mm 15% salt phantom using RCB Figure 4 23. Imaging Result with 4mm 15% salt phantom using SCB PAGE 91 9 1 Figure 4 24. Imaging Result with 4mm 15% salt phantom using DAS It is noticed that the target image obtained by SCB is smaller and isolated, w hile the image acquired from RCB is better having smoother boundary and clean background. On the other hand, the image from DAS is very noisy w ith many imaging artifacts. Actually RCB algorithm works between SCB and DAS algorithms. If the parameter in RCB is large, the result is more close to DAS. Otherwise, the result is more close to SCB. RCB results from phantom A (2mm, 15% salt), phantom B (4mm, 15% salt) and phantom C (6mm, 15% salt) are listed below for comparison with =0.3 Figure 4 25. Imaging Result with 2mm phantom using RCB PAGE 92 92 Figure 4 26. Imaging Result with 4mm phantom using RCB Figure 4 27. Imaging Result with 6mm phantom using RCB Regarding spatial resolution, m any factors could come into play. Th e bandwidth of transducer can af fect pulse width of thermoacoustic signal, which can blur t he image boundary and reduce spatial resolution. Large pulse width of microwave stimulation also PAGE 93 93 reduces spatial resolution. A lignment errors between phan tom and transducer array could hurt the con sistency of images for phantom during different experimental runs For rough estimates, we could get an approximate spatial resolution based on the simple reasoning below. Assume that pulse width of thermoacoustic pressure wave under delta heating is pulse width of microwave stimulation is and the frequency bandwidth of transducer is MHz. Roughly speaking, pulse width of the impulse response for transducer is = 1 According to convolution equation (2 21) and (2 22), pulse width of the piezo electric signal received by transducer could be estimated as below: = + + (4 1) and then spatial resolution could be estimated as: = (4 2) w here is spee d of sound. If we substitute microw ave pulse width 0.75 us and transducer frequency bandwidth 0.65 MHz into equations above we could end up with the spatial resolution at least 3.4 mm. PAGE 94 94 CHAPTER 5 CONCLUSION 5.1 Contribution In this thesis, a theoretic al and experimental study of thermo acoustic imaging system is presented, together wi th some research findings in hardware implementation of imaging algorithm, namely Robust Capon Beamforming algorithm. The initial imaging r esults are obtained based on experiment al data and dedicated imaging algorithm, which show the proof of principle study. Further research of the imag ing results is restricted by limitations of current experimental setup. T heoretical analysis about thermoa coustic signal generation by CW microwave stimulation is presented, including AM modulation technique and Diebold analysis [30] for thermoacous tic signal generation. After theoretical analysis experiment al results are presented with no observable thermoacoustic signal recorded. In reasoning for the absence of thermoacous tic signal, it is pointed out that peak power of CW microwave is too low to generate any observable signal compar ed with successful thermoacoustic signal generation under pulse d microwave stimulation In CW stimulation, heat conduction might hurt the reson ance of phantom through damping effects. RF safety concern also limits the highest CW RF power which can be used in medical application. As a result, we switched to pulsed microwave stimulation. Two FPGA implementations for narrow band and wideband RC B are performed in Xilinx Virtex4 SF35 platform device [52][53] which is specialized in signal processing application. Current implementation can successfully fit an 8 channel wideband RCB beamformer into the chip, while actual working algorithm used to g et initial imaging results is based on a 24 channel RCB beamformer In order to entirely put current PAGE 95 95 imaging algorithm in real time implementation, a larger FPGA device should be selected for implementation or some modifications should be done to reduce th e hardware utilization rate with penalty of longer latency. T he series of experiments for thermoacoustic imaging under pulsed microwave stimulation h as been performed including initial signa l generation experiments and later imaging experiments. The prel iminary images are obtained to serve as a successful proof of principle stud y. Several major limitations of current experimental platform are listed as follow. First, alignment between experimental tank (therefore phantom) and transducer ar ray subsystem is an issue, which can affect the imaging result in various aspects. Second, the bandwidth of transd ucer could be increased, and microwave pulse width could be reduced to achieve higher spatial resolution according to equation (4 1) N on uniform distribution of electromagnetic field inside experiment tank can be partially improved by some techniques, such as waveguide array power feeding, although it is expensive. 5.2 Future Work In previous section, major accomplishments and limitations of our thermoac oustic imaging research have been sum marized. Generally speaking, theoretical analysis, experimental investigation and practical implementation of imaging algorithm are involved. Based on limitations o f current system, fu ture works are laid out as follow. N ew acoustic sensor or MEMS sensor array with wider frequency response and smaller size will be deployed for signal reception. Shorter microwave pulse will be used, together with newer acoustic sensor, to improve spatial resolution of the image. New microw ave power system might be designed to improve the hom ogeneity of electromagnetic field dist ribution, which could reduce imaging artifacts. More aggressiv e PAGE 96 96 hardware implementation of RCB algorithm co uld be investigated to boost computation for image reconst ruction PAGE 97 97 LIST OF REFERENCES [1] Med. Phys ., vol. 27, no. 5, pp. 1195 1202, May. 2000. [2] brea st cancer detection system combining both thermoacoustic (TA) and Med. Phys ., vol. 35, no. 6, pp. 2218 2223, Jun. 2008. [3] Radio Waves: A Medical Im Radiology vol. 121, no. 1, pp. 275 278, Apr 1999. [4] PhD Thesis University of Florida, Gainesville, FL, 2007. [5] ermoacoustic and Photoacoustic Technology in Cancer Research & Treatment vol. 4, no. 5, pp. 559 565, Oct 2005. [6] Ultrasonics Symposium pp. 817 822, 1981. [7] IEEE Transactions on Microwave Theory and Techniques vol. 32, no. 8, pp. 854 860, Aug 1984. [8] L. V. Wang, X. Zhao, H. Sun and G. Ku, induced acoustic imaging of Review Of Scientific Instruments vol. 70, no. 9, pp. 3744 3748, Sep 1999. [9] induced thermoacoustic Med. Phys ., vol. 28, no. 1, pp. 4 10, Jan 2001. [10] Med. Phys ., vol. 28, no. 7, pp. 1519 1524, Jul 2001. [11] e induced thermoacoustic Med. Phys ., vol. 28, no. 12, pp. 2427 2431, Dec 2001. [12] induced thermoacoustic tomography using multi Med. Phys ., vol. 28, no. 9, pp. 1958 1963, Sep 2001. PAGE 98 98 [13] IEEE Trans. Medical Imaging vol. 21, no. 7, pp. 814 822, Jul 2002. [14] micr owave induced thermoacoustic tomography: Med. Phys ., vol. 29, no. 8, pp. 1661 1669, Aug 2002. [15] Thermoacoustic To IEEE Trans. Medical Imaging vol. 21, no. 7, pp. 823 828, Jul 2002. [16] IEEE Trans. Medical Imaging vol. 21, no. 7, pp. 829 833, Jul 2002. [17] IEEE Trans. Biomedical Engineering vol. 50, no. 9, pp. 1086 1 099, Sep 2003. [18] IEEE Trans. Biomedical Engineering vol. 55, no. 12, pp. 2741 2752, Dec 2008. [19] M. Pramanik, G. Ku and L. V. Wa Breast Cancer Detection System Combining Both Thermoacoustic and Med. Phys. vol. 35, no. 6, pp. 2218 2223, Jun 2008. [20] computed tomography Med. Phys. vol. 26, no. 9, pp. 1832 1837, Sep 1999. [21] Proc. SPIE. vol. 3916, no. 150, Jan 2000. [2 2] Proc. SPIE ., vol. 4256, Jan 2001. [23] Medical Imaging 2002: Physics of Medical Imagin g pp. 521 525, May 2002. [24] Proc. SPIE., vol. 3982, pp. 354 359, 2000. [25] reg istered Photoacoustic, Proc. SPIE., vol. 7564, 2010. PAGE 99 99 [26] induced Thermoacoustic Proc. SPIE., vol. 6047, 2006. [27] induced Thermoacoustic Imaging IEEE/ICME International Conference on Complex Medical Engineering pp. 989 993, 2007. [28] he Feasibility of Using Microwave induced Physics in Medicine and Biology vol. 55, pp. 5203 5212, 2010. [29] IEEE Trans. Signal Processing vol. 51, no. 7, pp. 1702 1715, Jul 2003. [30] Journal of the Acoustical Society of America vol. 84, no. 6 pp. 2245 2251, Dec 1988. [31] IEEE Standard for Safety Levels with Respect to Human Exposure to Radio Frequency Electromagnetic Fields, 3 kHz to 300 GHz pp. 6 15, 1999. [32] Artech House Publishers 1996. [33] Wiley 2006. [34] IEEE ASSP Magazine pp. 4 24, 1988. [35] C. Dick and F Time QRD Based Beamforming on an FPGA IEEE Fortieth Asilomar Conference on Signals, Systems and Computers 2006. [36] PhD Thesis Lund University, 2006. [37] A. Ir Calculation Core Using QRD UCSD Technical Report CS2009 0937D. [38] Beamformers Used for Smart A EuCAP Nov 2007. [39] like Algorithm for Parallel IEEE Transactions on Computers Sep 1993. PAGE 100 100 [40] Proceeding of ACM/SIGDA Sixth International Symposium on Field Programmable Gate Arrays pp. 192 200, Feb 1998. [41] Based MVDR Beamformer Using IEEE International Con ference on Communications 2007. [42] Princi pal Values of A Complex Matrix, Transactions of the American Mathematical Society Jan 1960. [43] PhD Thesis Rice University, Apr 1994. [44] The 13 th IEEE International Symposium on PIMRC 2002 vol. 1, pp. 120 124, Sep 2002. [45] Radio Research Instrument Co, Operational Manual of 70 kW RF Package 2010. [46] Olympus Panametrics NDT, Transducer Test Report Mar 2005. [47] Olympus Panametrics NDT, Ma nual 2005. [48] National Instruments, NI PXI 5122 High Resolution Digitizer Datasheet [49] Ultrasonic Symposium pp. 633 638, 1984 [50] US Patent File 6216025 Apr 2001. [51] Symmetric Matrix in Fixed Tenth Annual Symposium on Multimedia Comm unications and Signal Processing Nov 23, 2001. [52] domain 2011 IEEE International Symposium on Antennas and Propagation pp. 2849 2852, Jul 2011 [53] SPIE 2012 Medical Imaging Conference Feb 2011 PAGE 101 101 BIOGRAPHICAL SKETCH Xin Guan was born in Baoji City, Shaan Xi provinc e, P. R. China on 1982. Xin Guan received Bachelor of Science in Tsinghua University, Beijing, P. R. China, majoring in Engineer ing Physics. After that Xin Guan got Master of Science in Tsinghua Accelerator Lab in D epartment of E ngineering P hysics. During that time, Xin Guan studied and researched on the S band photo injector for microw ave linear accelerator. After the ma ster stud ies Xin Guan came to University of Florida, Gainesville, FL for PhD in the Elect rical and Computer Engineering D epartment. During his PhD stud ies Xin Guan was primarily involved in the exp erimental research on the Th ermoacoustic Imaging System. research interests lie in the medical imaging, microwave engineering, digital signal processing and hardware implementation of DSP algorithm. 