A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
- Getting Started
- 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
Attenuation Compensation Using Time Variant Filtering Example
|On this page…|
This example demonstrates how the acoustic attenuation present in the photoacoustic forward problem can be compensated for using time-variant filtering. It builds on the 2D Time Reversal Reconstruction For A Circular Sensor and Attenuation Compensation Using Time Reversal examples.
For a more detailed discussion of this example and the underlying techniques, see B. E. Treeby "Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering," J. Biomed. Opt., vol. 18, no. 3, p.036008, 2013.
Running the forward simulation
The recorded sensor data for an attenuating medium is simulated using
kspaceFirstOrder2D with the initial pressure distribution set to the Shepp Logan phantom, in the same way as the previous example. A circular Cartesian sensor mask with 200 sensor points is used to record the photoacoustic signals. After the simulation, random Gaussian noise is added using
addNoise to give a signal-to-noise ratio of 40 dB. The sensor data is then used to reconstruct the photoacoustic image using time reversal. The initial pressure distribution and reconstructed image are shown below. The reconstructed image is very blurred due to the loss of high frequency information, and the magnitude of small features in the image is significantly reduced. This is also visible in the profiles, where the black line is a profile through the initial pressure distribution at x = 0, and the red line is the same profile through the reconstructed image.
Attenuation compensation with a fixed cut-off frequency
In the Attenuation Compensation Using Time Reversal Example, the acoustic attenuation in the forward problem is corrected as part of the time-reversal image reconstruction procedure by reversing the sign of the absorption term in the numerical model, thereby causing the waves to grow rather than decay. For a homogeneous distribution of absorption parameters, it is also possible to compensate for acoustic attenuation directly in the recorded time series (i.e., before image reconstruction) using time-variant filtering. The general idea is that the signals are modified using a filter that depends on both time and frequency. Here this is achieved using a form of nonstationary convolution using the function
Analogous to attenuation compensation using time-reversal, if the time-variant filter is applied to experimental or noisy data, the high frequency content (where the signal-to-noise ratio is typically quite low) can become amplified and mask the desired features within the signals. To avoid this, the time-variant filter is regularised using a window with a time-variant cutoff frequency. The cutoff frequency can be calculated automatically (as discussed in the following section), or specified manually. In the example below, the cutoff frequency for the window is defined manually at 3 MHz. This matches the cutoff used in the Attenuation Compensation Using Time Reversal Example. The attenuation compensation is then applied to each time series in the sensor data, where the time-variant filter is constructed based on the specified sound speed and attenuation parameters. The attenuation is assumed to follow a frequency power law of the form
alpha_coeff is defined in units of [dB/(MHz^y cm)]. When using a fixed cutoff frequency, the compensation is extremely fast to apply.
% correct for acoustic attenuation using time variant filtering regularised % by a Tukey window with a fixed cutoff frequency sensor_data_comp = attenComp(sensor_data_lossy, kgrid.dt, medium.sound_speed,... medium.alpha_coeff, medium.alpha_power, 'FilterCutoff', [3e6 3e6]);
The sensor data after attenuation compensation is then used to reconstruct the photoacoustic image using time reversal. As shown below, the compensation makes the reconstructed image much sharper, and the magnitude of the edges of the phantom are now much closer to the initial pressure distribution. Choosing a fixed cutoff frequency for the regularising window in this way gives the same correction as using time-reversal with attenuation compensation (discussed in the previous example). However, the drawback of choosing the cutoff frequency in this way is that the frequency content of the recorded photoacoustic signals is not normally uniform over time. This results in a compromise between the level of attenuation compensation, and the level of noise in the reconstructed image. In this example, it is clear that noise has been introduced (visible as a rippled texture) that wasn't present in the initial pressure distribution. This noise is also visible in the profile.
Using a time-variant window
The use of time-variant filtering for attenuation compensation introduces the possibility of selecting a time-variant window to regularise the compensation based on the local frequency content of the signal. This is significant, because high-frequencies generated by deeper features will undergo more attenuation before reaching the detector compared to shallow features, as the waves have travelled further. This means that high-frequency components from shallow features (which appear early in the time signal) may be above the noise floor, while the same frequency components from deeper features (which appear later in the time signal) may not. The use of a time-variant window to regularise the attenuation compensation means the filter can adapt to the local frequency content of the recorded signal.
If the optional input
'FilterCutoff' is not specified, the function
attenComp automatically chooses a time-variant cutoff frequency for the window based on the average time-frequency distribution of the signals. The steps for computing the cutoff frequency are detailed in the documentation for
attenComp. An image of the average time-frequency distribution of the sensor data along with the extracted time-variant cutoff frequency before (yellow) and after (white dashed) smoothing are shown below. It is clear from this example (shown with a dynamic range of 40 dB) that the frequency content of the signals changes over time.
The reconstructed photoacoustic image after attenuation compensation using a window with a time-variant cutoff frequency is shown below. Again, compared to the uncompensated case, the reconstructed image is much sharper, and the magnitude of small features is closer to their expected values. Compared to the image after compensation using a fixed cutoff frequency of 3 MHz, the reconstructed image shown below is not as sharp. However, the noise has now been significantly suppressed.
Correcting the signals individually
Correcting the signals using the average time-frequency distribution has a distinct computational advantage in that the time-varying cutoff frequency and the time-variant filter only need to be calculated once. This makes it reasonably fast to apply, even to a large number of signals. However, the downside is that the time-varying cutoff frequency may not be optimal for each individual signal. In the reconstructed image shown above, this results in some artefacts in the reconstructed image, most noticeably at the inside bottom edge of the phantom.
As an alternative, it is also possible to apply attenuation compensation to the recorded time signals individually by calling
attenComp in a loop. This way, the cutoff frequency is calculated based on the time-frequency distribution of the individual time series. As the loop operations are completely independent, it is possible to parallelise the calculations using the parallel computing toolbox (set the flag
true in the example .m file).
% correct for acoustic attenuation using time variant filtering regularised % by a Tukey window with a time-variant cutoff frequency based on the % time-frequency distribution for each signal sensor_data_comp = zeros(size(sensor_data_lossy)); if ~USE_PARALLEL_COMPUTING_TOOLBOX % loop through signals for index = 1:size(sensor_data_lossy, 1) % correct signal sensor_data_comp(index, :) = attenComp(sensor_data_lossy(index, :), kgrid.dt,... medium.sound_speed, medium.alpha_coeff, medium.alpha_power, ... 'FrequencyMultiplier', 3); end else % start matlab pool parpool; % loop through signals using parfor parfor index = 1:size(sensor_data_lossy, 1) % correct signal sensor_data_comp(index, :) = attenComp(sensor_data_lossy(index, :), kgrid.dt,... medium.sound_speed, medium.alpha_coeff, medium.alpha_power, ... 'FrequencyMultiplier', 3, 'DisplayUpdates', false); end % delete matlab pool delete(gcp('nocreate')); end
The reconstructed image after attenuation compensation is shown below. The reconstructed image is much sharper without noticeable noise or artefacts, and the magnitude of the edges has been restored. The primary drawback of correcting the signals individually is that calculating the time-varying cutoff frequency and constructing the time-variant filter for each individual signal is computationally more expensive than using an average or fixed cutoff frequency.
|Attenuation Compensation Using Time Reversal||Automatic Sound Speed Selection|
© 2009-2014 Bradley Treeby and Ben Cox.