# A Method to Obtain the Transducers Impulse Response (TIR) in Photoacoustic Imaging
## 2 Basic Principles of Solving TIR  
Light absorbance results in transient local heating and a pressure increase in proportion to the absorbed light energy. This leads to the creation of pressure waves. The propagation of photoacoustic pressure waves may be given through a generic wave equation [30]. \[(\
abla^{2}-\frac{1}{v^{2}}\frac{\partial^{2}}{\partial t^{2}})x(\mathbf{r},t )=-\frac{\beta}{C_{p}}\frac{\partial}{\partial t}\delta(\mathbf{r},t), \tag{1}\]  
where \(\mathbf{r}\) denotes the coordinates of any position in 3D space, \(t\) is time, \(v\) is the speed of sound, \(\beta\) is the coefficient for expansion of isobaric volume, and \(C_{p}\) is constant pressure-specific heat. \(\delta\) is the heat source function. We know that TIR is affected not only by the temporal impulse response, but also by the spatial impulse response [22; 23]. Because the transducer used in the simulation in this paper is a point-like transducer [33; 34; 35], its size is close to a point and can be ignored. Moreover, the simulation in this paper is based on far-field detection, which means that the distance between the support and the transducer is much longer than the wavelength corresponding to the center frequency of the detection device. In this article, we only consider temporal TIR. The PAPS at time \(t\) can be expressed as  
\[x(t)=\frac{\beta}{4\pi C_{p}}\iiint\frac{d\mathbf{r}}{\mathbf{r}}\frac{ \partial\delta(\mathbf{r},t^{\prime})}{\partial t^{\prime}}\bigg{|}_{t^{ \prime}=t-\frac{t}{v}}, \tag{2}\]  
where \(r=\left\lvert\mathbf{r}\right\rvert.\) If the transducer is limited by the bandwidth, then the OPAS is affected by the TIR and the obtained OPAS signal at the transducer becomes [36]  
\[y(t)=h(t)*x(t),\]where \(y(t)\) is the OPAS, \(x(t)\) is the PAPS, and \(h(t)\) is the TIR in this PAI system. To simplify the convolution operation in the time domain, it is transformed to the frequency domain through Fourier transform. So, when the number of transducers is \(J\) and by assuming that all the transducers have the same TIR (in this article, the transducers we use are made by the same manufacturer and the same model; all TIRs are considered to be the same), the following equation can be obtained  
\[H(\omega)=\frac{Y_{1}(\omega)}{X_{1}(\omega)}=\frac{Y_{2}(\omega)}{X_{2}( \omega)}=\cdots=\frac{Y_{J}(\omega)}{X_{J}(\omega)}. \tag{4}\]  
Here, \(H(\omega)\), \(Y_{j}(\omega)\), and \(X_{j}(\omega)\) are the TIR, the OPAS, and the PAPS after the Fourier transform, respectively, \(j=1\), \(2\), \(\ldots\), \(J\). Suppose there are only two transducers and \(j=1\), \(2\); then, the equation is  
\[Y_{2}(\omega)X_{1}(\omega)-Y_{1}(\omega)X_{2}(\omega)=0. \tag{5}\]  
Perform the inverse Fourier transform of Equation (5) to obtain the equation for the time domain:  
\[y_{2}(t)*x_{1}(t)-y_{1}(t)*x_{2}(t)=0. \tag{6}\]  
Discretize the continuous function of the time domain:  
\[y_{2}(k)*x_{1}(k)-y_{1}(k)*x_{2}(k)=0. \tag{7}\]  
Then, write the above equation in matrix form:  
\[Y_{2}\mathbf{x}_{1}-Y_{1}\mathbf{x}_{2}=0, \tag{8}\]  
where matrices \(\mathbf{Y_{1}}\) and \(\mathbf{Y_{2}}\) are the called convolutional matrices, or simply C-matrices. \[Y_{j}(0)= \begin{pmatrix}a_{1,j}&0&\ldots&0&0\\ a_{2,j}&a_{1,j}&\ldots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\ldots&a_{M,j}&a_{M-1,j}\\ 0&0&\ldots&0&a_{M,j}\end{pmatrix}  \tag{9}\]  
where \(a_{j}\) is the OPAS data from the \(j\)-th sensor. The first column of the matrix consists of \(M\) data of the OPAS, appended \(K-1\) zeros, \(M\) is the size of the OPAS vector, and \(K\) is the size of the PAPS vector. The other columns are cyclic copies of the first column. Similarly, \(L\) is the size of the TIR vector \(\mathbf{h}\).