k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
- Getting Started
- Examples
- Initial Value Problems
- Example: Homogenous Propagation Medium
- Example: Using A Binary Sensor Mask
- Example: Defining A Sensor Mask By Opposing Corners
- Example: Loading External Image Maps
- Example: Heterogeneous Propagation Medium
- Example: Saving Movie Files
- Example: Recording The Particle Velocity
- Example: Defining A Gaussian Sensor Frequency Response
- Example: Comparison Of Modelling Functions
- Example: Setting An Initial Pressure Gradient
- Example: Simulations In One Dimension
- Example: Simulations In Three Dimensions
- Example: Photoacoustic Waveforms in 1D, 2D and 3D
- Time Varying Source Problems
- Example: Monopole Point Source In A Homogeneous Propagation Medium
- Example: Dipole Point Source In A Homogeneous Propagation Medium
- Example: Simulating Transducer Field Patterns
- Example: Steering A Linear Array
- Example: Snell's Law And Critical Angle Reflection
- Example: The Doppler Effect
- Example: Diffraction Through A Slit
- Example: Simulations In Three-Dimensions
- Sensor Directivity
- Example: Focussed Detector in 2D
- Example: Focussed Detector in 3D
- Example: Modelling Sensor Directivity in 2D
- Example: Modelling Sensor Directivity in 3D
- Example: Sensor Element Directivity in 2D
- Example: Focussed 2D Array with Directional Elements
- Photoacoustic Image Reconstruction
- Example: 2D FFT Reconstruction For A Line Sensor
- Example: 3D FFT Reconstruction For A Planar Sensor
- Example: 2D Time Reversal For A Line Sensor
- Example: 2D Time Reversal For A Circular Sensor
- Example: 3D Time Reversal For A Planar Sensor
- Example: 3D Time Reversal For A Spherical Sensor
- Example: Image Reconstruction With Directional Sensors
- Example: Image Reconstruction With Bandlimited Sensors
- Example: Iterative Image Improvement Using Time Reversal
- Example: Attenuation Compensation Using Time Reversal
- Example: Attenuation Compensation Using Time Variant Filtering
- Example: Automatic Sound Speed Selection
- Diagnostic Ultrasound Simulation
- Example: Defining An Ultrasound Transducer
- Example: Simulating Ultrasound Beam Patterns
- Example: Using An Ultrasound Transducer As A Sensor
- Example: Simulating B-mode Ultrasound Images
- Example: Simulating B-mode Images Using A Phased Array
- Numerical Analysis
- Example: Controlling The Absorbing Boundary Layer
- Example: Source Smoothing
- Example: Filtering A Delta Function Input Signal
- Example: Modelling Power Law Absorption
- Example: Modelling Nonlinear Wave Propagation
- Example: Optimising k-Wave Performance
- Using The C++ Code
- Elastic Wave Propagation
- Example: Explosive Source In A Layered Medium
- Example: Plane Wave Absorption
- Example: Shear Waves And Critical Angle Reflection
- Example: Simulations In Three Dimensions
- Functions - By Category
- Functions - Alphabetical List
- Release Notes
- License
k-Wave Toolbox |
Attenuation Compensation Using Time Reversal Example
On this page… |
---|
Running the forward simulation |
Overview
This example demonstrates how the acoustic attenuation present in photoacoustic forward problem can be compensated for using time reversal image reconstruction. It builds on the 2D Time Reversal Reconstruction For A Circular Sensor Example.
For a more detailed discussion of this example and the underlying techniques, see B. E. Treeby, E. Z. Zhang, and B. T. Cox, "Photoacoustic tomography in absorbing acoustic media using time reversal," Inverse Problems, vol. 26, no. 11, p. 115003, 2010.
Running the forward simulation
The sensor data is simulated using kspaceFirstOrder2D
in a similar manner to previous examples. The initial pressure is set to the Shepp Logan phantom and a circular Cartesian sensor mask with 200 sensor points is used to record the data. After simulation, random Gaussian noise is added using addNoise
to give a signal to noise ratio of 40dB (based on the peak of the recorded signal and the rms noise level).
% run the forward simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % add noise to the recorded sensor data signal_to_noise_ratio = 40; % [dB] sensor_data = addNoise(sensor_data, signal_to_noise_ratio, 'peak');
Selecting the regularisation parameters
To compensate for the acoustic attenuation in the forward problem, the absorption parameter must be reversed in sign (the frequency content must grow rather than decay). This is achieved by setting medium.alpha_sign
to [-1, 1]
. The two inputs control the sign of the absorption and dispersion parameters, respectively. The dispersion parameter (which controls the dependence of the sound speed on frequency) should not be reversed because if the high frequency components of the wave field have travelled to the detector faster than the low frequency components (as is the case for photoacoustic signals in biological tissue), they again need to travel faster than the low frequency components in the time reversal reconstruction to regain their initial position within the medium.
If a photoacoustic image is reconstructed from real or noisy data with the absorption parameter reversed, the high frequency content (where the signal to noise ratio is typically quite low) can quickly grow to mask the low frequency content, effectively obscuring the desired features within the reconstructed image. This is because, as the waves propagate, the high frequencies are increased at a much faster rate (recall the acoustic attenuation in soft biological tissue follows a frequency power law). To avoid this, the absorption and dispersion parameters are filtered using a frequency domain Tukey window. This regularises the reconstruction (in effect, stops it from 'blowing up') by restricting the range of frequencies that are allowed to grow. The filter is applied to the absorption parameters by assigning it to medium.alpha_filter
. The filter must have the same number of dimensions and grid points as kgrid.k
.
To choose an appropriate filter cutoff frequency, the average power spectrum of the simulated sensor data is computed. This is plotted below (shown in black), along with the power spectrum of the signal recorded at the first sensor point (shown in red). In this example, the filter cutoff frequency (shown as the left dashed line in the figure) is chosen based on the noise floor observable in the power spectrum. The maximum frequency supported by the grid is also shown (right dashed line in the figure).
Here, a symmetric Tukey window (or tapered cosine window) scaled to a particular frequency cutoff is used so that the correct power law absorption and dispersion characteristics are maintained within the filter pass band. The filter is created using getAlphaFilter
which creates the filter using getWin
via rotation.
% define the cutoff frequency for the filter f_cutoff = 3e6; % create the filter to regularise the absorption parameters medium.alpha_filter = getAlphaFilter(kgrid_recon, medium, f_cutoff); % reverse the sign of the absorption proportionality coefficient medium.alpha_sign = [-1, 1]; % [absorption, dispersion];
A visualisation of the filter is given below.
Time reversal image reconstruction
The initial pressure distribution, and the reconstructions with and without compensation for acoustic attenuation are shown below. Profiles through x = 0 are also shown for comparison. Without correction for attenuation, the edges of the reconstructed pressure become blurred, and the overall magnitude is reduced. When compensation for acoustic attenuation is included into the reconstruction, the sharpness and magnitude of the reconstructed pressure is considerably improved. Due to the band-limited frequency response of the sensors (high frequencies are not detected because their magnitudes are below the noise floor of the added noise), the reconstruction is unable to completely recover the frequency content of the initial pressure distribution. This results in Gibbs phenomenon appearing in the reconstruction where there are sharp changes in gradient.
Iterative Image Improvement Using Time Reversal | Attenuation Compensation Using Time Variant Filtering |
© 2009-2014 Bradley Treeby and Ben Cox.