Comparison of intensity-modulated continuous-wave lasers with a chirped modulation frequency to pulsed lasers for photoacoustic imaging applications

Adam Petschke; Patrick J. La Rivière

doi:10.1364/BOE.1.001188

1. Introduction

The photoacoustic effect (which is also called the optoacoustic effect) is the generation of an acoustic wave by thermoelastic expansion caused by heating from optical absorption. Photoacoustic imaging is a hybrid imaging system that combines the advantage of optical absorption contrast with ultrasound resolution. In addition, since the absorption properties of tissue are influenced by biological activity, photoacoustic systems can perform functional imaging. In particular, by using multiple excitation wavelengths it is possible to extract quantitative information [1]. Most photoacoustic imaging systems use high-energy pulsed lasers, which tend to be expensive and bulky, for excitation. The pulse duration in these systems is usually about 10 ns. However, it is also possible to use intensity-modulated continuous-wave (CW) lasers. One choice of modulation is to chirp the modulation frequency [2–7]. The pulse duration in these systems is usually about 1 ms. One potential advantage of chirped photoacoustic imaging systems is that they can use compact and inexpensive CW diode lasers. The aim of this paper is to compare these chirped systems to pulsed systems.

2. Pulse compression using matched filters

To produce a chirp signal, the CW laser is modulated to create an irradiance (or fluence rate) that has a time-swept frequency. The irradiance is given by

I_{chirp} (t) = I_{0} [1 + \cos (ω_{0} t + π b t^{2})] rect (\frac{t}{T}),

where I ₀ is the average irradiance of the pulse (2I ₀ is the peak irradiance) with units W/m², b is the frequency sweep rate with units s⁻², and T is the length of the pulse with units s. The fluence of the pulse is I ₀ T with units J/m². The angular frequency is given by d/dt(ω ₀ t + πbt ²) = ω ₀ + 2πbt, which gives the frequency, f, as f ₀−bT/2 ≤ f ≤ f ₀ + bT/2, the bandwidth as bT, and the time-bandwidth product as bT ². For example, if the chirp goes from 1 MHz to 5 MHz and T is 1 ms, then f ₀ is 3 MHz and b is 4×10⁹ s⁻².

Eq. (1) can be split up into a low-frequency term (a pulse with constant amplitude I ₀ and length T) and a high-frequency term (a pulse with amplitude I ₀ cos(ω ₀ t + πbt ²) and length T). In most practical situations, the low-frequency term is outside of the bandwidth of the acoustic transducer since T is on the order of 1 ms and the transducer has a bandpass frequency response with a center frequency in the MHz range. The Fourier transform of the high frequency term can be split up into positive and negative frequency components. By completing the square, the positive frequency components (using the unitary Fourier transform) are given by

{\tilde{I}}_{pos} (ω) = \frac{I_{0}}{2 \sqrt{2 π}} e^{- i \frac{π}{b} {(f - f_{0})}^{2}} \frac{1}{\sqrt{2 b}} \int_{- T \sqrt{\frac{b}{2}} - \sqrt{\frac{2}{b}} (f - f_{0})}^{T \sqrt{\frac{b}{2}} - \sqrt{\frac{2}{b}} (f - f_{0})} e^{i \frac{π}{2} {(t')}^{2}} d t',

which can be solved in terms of the Fresnel integrals

C (x) = \int_{0}^{x} \cos (\frac{π}{2} t^{2}) d t

with the properties C(0) = S(0) = 0 and C(±∞) = S(±∞) = ±1/2 [8]. Eq. (2) then becomes

{\tilde{I}}_{pos} (ω) = \frac{I_{0}}{2 \sqrt{2 π}} e^{- i \frac{π}{b} {(f - f_{0})}^{2}} \frac{1}{\sqrt{2 b}} [C (T \sqrt{\frac{b}{2}} - \sqrt{\frac{2}{b}} (f - f_{0})) + i S (T \sqrt{\frac{b}{2}} - \sqrt{\frac{2}{b}} (f - f_{0}))

which resembles a rectangular function centered at f ₀ with bandwidth bT for large time-bandwidth products (bT ² > 30) [9]. The spectrum of a chirp signal is shown in Fig. 1(a) for f ₀ = 3 MHz and bT = 4 MHz. Eq. (4) is then approximately given by

{\tilde{I}}_{pos} (ω) \approx \frac{I_{0}}{2 \sqrt{2 π}} e^{\frac{- i {(ω - ω_{0})}^{2}}{4 π b}} \frac{1}{\sqrt{2 b}} (1 + i) rect (\frac{ω - ω_{0}}{2 π b T}) .

Similarly, the negative frequency components are approximately given by

{\tilde{I}}_{neg} (ω) \approx \frac{I_{0}}{2 \sqrt{2 π}} e^{\frac{i {(ω + ω_{0})}^{2}}{4 π b}} \frac{1}{\sqrt{2 b}} (1 - i) rect (\frac{ω + ω_{0}}{2 π b T}) .

A matched filter is defined as the complex conjugate of the spectrum of the signal that is being filtered. It can be shown that applying a matched filter to a signal will give the best signal-to-noise ratio (SNR) when the signal is corrupted by white noise [9]. In our case the unity-gain matched filter is

{\tilde{G}}_{chirp} (ω) = e^{\frac{i {(ω - ω_{0})}^{2}}{4 π b}} \frac{1 - i}{\sqrt{2}} rect (\frac{ω - ω_{0}}{2 π b T}) + e^{\frac{- i {(ω + ω_{0})}^{2}}{4 π b}} \frac{1 + i}{\sqrt{2}} rect (\frac{ω + ω_{0}}{2 π b T}),

and the filtered chirp spectrum is given by

{\tilde{I}}_{fc} (ω) = \frac{I_{0}}{2 \sqrt{2 π}} \frac{1}{\sqrt{b}} [rect (\frac{ω - ω_{0}}{2 π b T}) + rect (\frac{ω + ω_{0}}{2 π b T})] .

Taking the inverse Fourier transform gives

I_{fc} (t) = I_{0} \sqrt{{bT}^{2}} sinc (π b T t) \cos (ω_{0} t),

where sinc(πbTt) = sin(πbTt)/(πbTt). This is shown in Fig. 1(b) for f ₀ = 3 MHz and bT = 4 MHz. The filtered signal has undergone pulse compression by a factor 1/(bT ²) and the axial resolution of the photoacoustic imaging system (ignoring the frequency dependence of photoacoustic generation and the transducer bandwidth, which are discussed later) is approximately ν_s/(bT) where ν_s is the speed of sound (~ 1500 m/s in tissue). The signal is compressed by delaying different frequency components by different times. Since the chirp pulse has an instantaneous frequency that depends on time, the pulse will increase in amplitude and decrease in duration if the delay is done in the same manner as the frequency sweep, which is exactly what the matched filter accomplishes.

3. Green’s function solution to the photoacoustic equation

In order to solve the photoacoustic wave equation, we start with the heat conduction equation given by

\frac{\partial T (r, t)}{\partial t} - D_{T} \nabla^{2} T (r, t) = \frac{H (r, t)}{ρ C_{V}},

where T(r, t) is the temperature with units K, D_T is the thermal diffusivity (~ 1×10⁻⁷ m²/s for tissue), H(r, t) is the heating function with units W/m³, ρ is the mass density (~ 1000 kg/m³ for tissue), and C_V is the specific heat at constant volume (~4000 J/(kg·K) for tissue). The first two terms in Eq. (10) go like ωT and D_TT/l ², respectively, where l is the characteristic length of the absorber, which means we can ignore the second term in most practical situations [10]. Eq. (10) then becomes

Fig. 1. (a) The spectrum and (b) time signal, compressed by applying a matched filter, for a chirp signal with f ₀ = 3 MHz and bT = 4 MHz.

Download Full Size | PDF

\frac{\partial T (r, t)}{\partial t} = \frac{H (r, t)}{ρ C_{V}} .

The general photoacoustic equation [11] is

(\nabla^{2} - \frac{1}{ν_{s}^{2}} \frac{\partial^{2}}{\partial t^{2}}) p (r, t) = - \frac{β}{κ ν_{s}^{2}} \frac{\partial^{2} T (r, t)}{\partial t^{2}},

where p(r, t) is the pressure with units Pa, β is the thermal coefficient of volume expansion (~ 2×10⁻⁴ K⁻¹ for tissue), and κ is the isothermal compressibility (~ 5×10⁻¹⁰ Pa⁻¹ for tissue). We use κ = C_P/(ρC_Vν ² _s), where C_P is the specific heat at constant pressure (~ 4000 J/(kg·K) for tissue), Eq. (11), and Eq. (12) to get

(\nabla^{2} - \frac{1}{v_{s}^{2}} \frac{\partial^{2}}{\partial t^{2}}) p (r, t) = - \frac{β}{C_{P}} \frac{\partial H (r, t)}{\partial t} .

In order to solve Eq. (13), we take the Fourier transform to get

(\nabla^{2} + k^{2}) \tilde{p} (r, ω) = - \frac{i ω β}{C_{P}} \tilde{H} (r, ω),

where k = −ω/ν_s. The Green’s function solution is

\tilde{p} (r, ω) = \frac{i ω β}{4 π C_{P}} \int_{V} \frac{e^{ik ∣ r - r_{0} ∣}}{∣ r - r_{0} ∣} \tilde{H} (r_{0}, ω) d^{3} r_{0} .

In the far-field, where r ≫ r ₀ we have |r−r₀| ≈ r − r̂ · r₀. We keep both terms for the phase but only the first term for the amplitude. Eq. (15) becomes

\tilde{p} (r ≫ r_{0}, ω) = \frac{i ω β}{4 π C_{P}} \frac{e^{ikr}}{r} \int_{V} e^{- ik \hat{r} \cdot r_{0}} \tilde{H} (r_{0}, ω) d^{3} r_{0} .

We now consider a cubic absorber, and we are interested in the pressure at the point (x = 0,y = 0,z ≫ a). The heating function for a cubic absorber has the form

H (x_{0}, y_{0}, z_{0}, t) = {\begin{matrix} \begin{matrix} χ μ e^{- μ (\frac{a}{2 - z_{0}})} I (t) & - \frac{a}{2} \leq x_{0}, y_{0}, z_{0} \leq \frac{a}{2} \end{matrix} \\ \begin{matrix} 0 & otherwise, \end{matrix} \end{matrix}

where χ is the dimensionless efficiency of absorbed energy converted to heat, µ is the absorption coefficient with units m⁻¹, $e^{- μ (\frac{a}{2 - z_{0}})}$ represents the decrease in irradiance due to the absorption of a plane wave light source incident from the positive z axis, I(t) is the irradiance of the plane wave light source, and a is the length of one side of the absorber with units m. This is shown in Fig. 2. Using r̂ = ẑ = ẑ₀ and k = −ω/ν_s, Eq. (16) becomes

\tilde{p} (ω) = \frac{β χ μ v_{s} a^{2}}{4 π z C_{P}} \tilde{I} (ω) \frac{ω}{ω - i μ v_{s}} (e^{- i ω (\frac{2 z - a}{2 v_{s}})} - e^{- μ a} e^{- i ω (\frac{2 z + a}{2 v_{s}})}),

where p̃(ω) is shorthand for p̃(x = 0,y=0,z ≫ a, ω). Even though photoacoustic generation is less efficient at high frequencies [7], we notice that the pressure is maximum in the far-field for high frequencies because diffraction scales inversely with frequency. However, in this analysis we are ignoring absorption in the surrounding medium, so using the highest possible frequency will not always be the best choice.

Fig. 2. Illustration of coordinate system and geometry of absorber.

Download Full Size | PDF

4. Far-field pressure for a pulse excitation

For a pulse given by

I_{pulse} (t) = \frac{F_{0}}{τ} rect (\frac{t}{τ}),

where F ₀ is the fluence per pulse and τ is the pulse length with units s, we have

{\tilde{I}}_{pulse} (ω) = \frac{F_{0}}{\sqrt{2 π}} sinc (\frac{ω τ}{2}) .

Taking the inverse Fourier transform of Eq. (18) and using Eq. (20) we have

p_{pulse} (t) = \frac{β χ μ v_{s} a^{2} F_{0}}{4 π z τ C_{P}} [e^{- μ v_{s} t_{1}} u (t_{1}) - e^{- μ v_{s} t_{2}} u (t_{2}) - e^{- μ a} e^{- μ v_{s} t_{3}} u (t_{3}) + e^{- μ a} e^{- μ v_{s} t_{4}} u (t_{4})],

where p _pulse(t) is shorthand for p _pulse(x = 0,y = 0,z ≫ a, t), u(t) is the unit step function, t ₁ = t + τ/2 − (2z − a)/(2ν_s), t ₂ = t − τ/2 − (2z−a)/(2ν_s), t ₃ = t + τ/2 − (2z+a)/(2ν_s), and t ₄ = t − τ/2 − (2z + a)/(2ν_s). This is shown in Fig. 3 for χ = 0.25, µ = 1 cm⁻¹, a = 100 µm, F ₀ = 20 mJ/cm², z = 3 cm, and τ = 10 ns.

Fig. 3. Pressure at a distance of 3 cm for a 10 ns pulse excitation. The temporal profiles of each of the pressure pulses are almost exact replicas of the laser pulse for the parameters chosen.

Download Full Size | PDF

5. Comparison of SNRs

If we also include the frequency response of an acoustic transducer, T̃(ω), and a unity-gain matched filter, G̃(ω), in Eq. (18), we have

\tilde{S} (ω) = \tilde{A} (ω) \tilde{I} (ω) \tilde{T} (ω) \tilde{G} (ω),

where S̃(ω) is the detected signal spectrum and

\tilde{A} (ω) = \frac{β χ μ v_{s} a^{2}}{4 π z C_{P}} \frac{ω}{ω - i μ v_{s}} (e^{- i ω (\frac{2 z - a}{2 v_{s}})} - e^{- μ a} e^{- i ω (\frac{2 z + a}{2 v_{s}})}) .

Note that the sensitivity of the transducer is included in T̃(ω). If the pulse is short enough, Ĩ _pulse(ω) in Eq. (20) will be approximately constant with amplitude $\frac{F_{0}}{\sqrt{2 π}}$ over the transducer bandwidth so that ${\tilde{I}}_{pulse} (ω) \tilde{T} (ω) \approx \frac{F_{0} \tilde{T} (ω)}{\sqrt{2 π}}$ . Eq. (22) then becomes

{\tilde{S}}_{pulse} (ω) = \frac{F_{0}}{\sqrt{2 π}} \tilde{A} (ω) \tilde{T} (ω) [rect (\frac{ω - ω_{0}}{2 π Δ f}) + rect (\frac{ω + ω_{0}}{2 π Δ f})],

where ω ₀ is the transducer center frequency, Δf is the transducer bandwidth, and we have used

{\tilde{G}}_{pulse} (ω) = rect (\frac{ω - ω_{0}}{2 π Δ f}) + rect (\frac{ω + ω_{0}}{2 π Δ f})

to reject out of band noise. Note that this analysis is general and can be applied to any definition of the transducer bandwidth and any T̃(ω); however, G̃ _pulse(ω) is not matched exactly unless the transducer spectrum is a rect and Ã(ω) is constant over the transducer bandwidth.

In order to make a fair comparison, the parameters for the chirp signal should be chosen so that the chirp bandwidth matches the transducer bandwidth (using the same definition of Δf as in the short pulse case), in which case we have bT = Δf. Using Eq. (8), where Ĩ _fc(ω) = Ĩ _chirp(ω)G̃ _chirp(ω), we can write Eq. (22) as

{\tilde{S}}_{chirp} (ω) = \frac{I_{0}}{2 \sqrt{2 π}} \frac{1}{\sqrt{b}} \tilde{A} (ω) \tilde{T} (ω) [rect (\frac{ω - ω_{0}}{2 π Δ f}) + rect (\frac{ω + ω_{0}}{2 π Δ f})] .

As with the short pulse case, G̃ _chirp(ω) is not matched exactly unless the transducer spectrum is a rect and Ã(ω) is constant over the transducer bandwidth. We see that for both the short pulse and the chirp signals, the axial resolution of the imaging system is the same and is determined by the transducer and Ã(ω), which is dependent on the material properties of the object. Note that for actual objects, Ã(ω) will be more complicated than what is given in Eq. (23). If we ignore Ã(ω), then the resolution is approximately ν_s/Δf.

The noise amplitude, Ñ(ω), is filtered by |G̃(ω)|. This means the filtered noise amplitudes are the same in both cases and are given by

{\tilde{N}}_{pulse} (ω) = {\tilde{N}}_{chirp} (ω) = \tilde{N} (ω) [rect (\frac{ω - ω_{0}}{2 π Δ f}) + rect (\frac{ω + ω_{0}}{2 π Δ f})] .

Defining SNR to be the signal amplitude divided by the noise amplitude gives

\frac{{SNR}_{chirp}}{{SNR}_{pulse}} = \frac{I_{0} T}{2 F_{0} \sqrt{b T^{2}}} .

We see that the ratio of SNRs is determined by the fluence of the chirp pulse, the fluence of the short pulse, and the time-bandwidth product of the chirp pulse.

6. SNR in the context of the ANSI limits

The American National Standards Institute (ANSI) laser safety limits for skin [12] for a short pulse (1 ns ≤ τ ≤ 100 ns), as used in pulsed systems, are given by

min {\begin{matrix} 20 C_{A} \\ {\begin{matrix} \begin{matrix} \frac{1100 C_{A} t^{\frac{1}{4}}}{N_{τ}} & t \leq 10 s \end{matrix} \\ \begin{matrix} \frac{200 C_{A} t}{N_{τ}} & t \geq 10 s, \end{matrix} \end{matrix} \end{matrix}

and the ANSI limits for a longer pulse (100 ns ≤ T ≤ 10 s), as used in chirped systems, are given by

min {\begin{matrix} 1100 C_{A} T^{\frac{1}{4}} \\ {\begin{matrix} \begin{matrix} \frac{1100 C_{A} t^{\frac{1}{4}}}{N_{T}} & t \leq 10 s \end{matrix} \\ \begin{matrix} \frac{200 C_{A} t}{N_{T}} & t \geq 10 s \end{matrix} \end{matrix} \end{matrix}

in units of mJ/cm², where t is the total exposure time in seconds, N is the total number of pulses (N_τ = N _pulse = R_τt and N_T = N _chirp = R_Tt, where R_τ is the short pulse repetition rate and R_T is the chirp pulse repetition rate), and C_A is a wavelength correction factor defined in [12]. The maximum F ₀ for a short pulse is 20C_A mJ/cm² when R_τ is set to the optimum rate given by 55/t ^3/4 in units of Hz for t ≤ 10 s and 10 Hz for t ⁥ 10 s. The maximum I ₀ for a chirp pulse is 1100C_A/T ^3/4 mW/cm² when R_T is set to the optimum rate given by 1/(T ^1/4 t ^3/4) in units of Hz for t ≤ 10 s and 10/(55T ^1/4) in units of Hz for t ⁥ 10 s. The optimum rates are plotted in Fig. 4 for a short pulse and a 1 ms chirp pulse. With multiple pulses, we can perform averaging, which improves the SNR by √N. Note that the SNR is also maximized by choosing the same repetition rates that optimize the fluence and irradiance. Including averaging and assuming the optimum repetition rates are used, Eq. (28) becomes

\frac{{SNR}_{chirp}}{{SNR}_{pulse}} = \frac{\sqrt{55}}{2 T^{\frac{3}{8}} \sqrt{Δ f}} .

For typical values of the relevant parameters, we see that the SNR is about 20 dB to 30 dB larger for pulsed laser systems. In light of Eq. (9), it might seem counter-intuitive that shorter chirp pulses give a better SNR; however, the noise and ANSI limits also depend on T in such a way as to create that dependence. In the following discussion, we assume white noise for ease of explanation. Since white noise depends on $\sqrt{Δ f} = \sqrt{b T}$ , we have ${SNR}_{chirp} \propto I_{0} T \frac{\sqrt{N_{T}}}{\sqrt{T}} \propto \frac{1}{T^{\frac{3}{8}}}$ because the fluence is proportional to T ^1/4. However, for applications where the fluence is not determined by the ANSI limits, the maximum repetition rate is 1/T and so SNR_chirp ∝ I ₀√t and does not depend on T. In both cases, SNR_chirp does not depend on the bandwidth. The $\sqrt{Δ f}$ in the denominator of Eq. (31) arises because ${SNR}_{pulse} \propto F_{0} Δ f \frac{\sqrt{N_{τ}}}{\sqrt{Δ f}} \propto \sqrt{Δ f}$ .

Fig. 4. The optimum repetition rates for a short pulse and a 1 ms chirp pulse.

Download Full Size | PDF

7. Conclusion

We have shown that the SNR of photoacoustic imaging systems based on CW lasers with a chirped modulation frequency are about 20 dB to 30 dB worse than systems based on pulsed lasers if we are constrained by the ANSI safety limits. We have also shown that both systems have the same resolution. However, the chirp based systems have the advantage of being able to employ CW diode lasers. This advantage could be especially important for spectroscopic studies where it could be possible to image simultaneously at multiple wavelengths using a compact and inexpensive system.

Acknowledgments

This work was supported by ACS Research Scholar grant RSG-08-117-01-CCE to Patrick La Riviére. Adam Petschke acknowledges the support of a Paul C. Hodges Award from the Paul C. Hodges Alumni Society of the Department of Radiology at the University of Chicago.

References and links

1. C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54(19), R59–R97 (2009). [CrossRef] [PubMed]

2. Y. Fan, A. Mandelis, G. Spirou, and I. A. Vitkin, “Development of a laser photothermoacoustic frequency-swept system for subsurface imaging: theory and experiment,” J. Acoust. Soc. Am. 116(6), 3523–3533 (2004). [CrossRef] [PubMed]

3. Y. Fan, A. Mandelis, G. Spirou, I. A. Vitkin, and W. M. Whelan, “Laser photothermoacoustic heterodyned lock-in depth profilometry in turbid tissue phantoms,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 72(5), 051908 (2005). [CrossRef] [PubMed]

4. S. A. Telenkov and A. Mandelis, “Fourier-domain biophotoacoustic subsurface depth selective amplitude and phase imaging of turbid phantoms and biological tissue,” J. Biomed. Opt. 11(4), 044006 (2006). [CrossRef] [PubMed]

5. S. A. Telenkov and A. Mandelis, “Fourier-domain methodology for depth-selective photothermoacoustic imaging of tissue chromophores,” Eur. Phys. J. Spec. Top. 153(1), 443–448 (2008). [CrossRef]

6. S. Telenkov, A. Mandelis, B. Lashkari, and M. Forcht, “Frequency-domain photothermoacoustics: Alternative imaging modality of biological tissues,” J. Appl. Phys. 105(10), 102029 (2009). [CrossRef]

7. S. A. Telenkov and A. Mandelis, “Photothermoacoustic imaging of biological tissues: maximum depth characterization comparison of time and frequency-domain measurements,” J. Biomed. Opt. 14(4), 044025 (2009). [CrossRef] [PubMed]

8. H. H. Barrett, and K. J. Myers, Foundations of Image Science (Wiley, Hoboken, NJ, 2004).

9. C. E. Cook, and M. Bernfeld, Radar Signals: An Introduction to Theory and Application (Academic, New York, NY, 1967).

10. R. A. Kruger, P. Liu, Y. R. Fang, and C. R. Appledorn, “Photoacoustic ultrasound (PAUS)--reconstruction tomography,” Med. Phys. 22(10), 1605–1609 (1995). [CrossRef] [PubMed]

11. L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron. 14(1), 171–179 (2008). [CrossRef]

12. Laser Institute of America, American National Standard for Safe Use of Lasers ANSI Z136.1–2007 (American National Standards Institute, Orlando, FL, 2007).

Comparison of intensity-modulated continuous-wave lasers with a chirped modulation frequency to pulsed lasers for photoacoustic imaging applications

Abstract

1. Introduction

2. Pulse compression using matched filters

3. Green’s function solution to the photoacoustic equation

4. Far-field pressure for a pulse excitation

5. Comparison of SNRs

6. SNR in the context of the ANSI limits

7. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Equations (33)

Biomedical Optics Express