attenComp
Attenuation compensation using time-variant filtering.
Syntax
signal = attenComp(signal, dt, c, alpha_0, y) signal = attenComp(signal, dt, c, alpha_0, y, ...) [signal, tfd, cutoff_freq] = attenComp(signal, dt, c, alpha_0, y) [signal, tfd, cutoff_freq] = attenComp(signal, dt, c, alpha_0, y, ...)
Description
attenComp
corrects for frequency dependent acoustic attenuation in photoacoustic signals using time-variant filtering [1]. The time-variant filter is constructed to correct for acoustic attenuation and dispersion following a frequency power law of the form alpha_0*f^y
under the assumption the distribution of attenuation parameters is homogeneous. The filter is applied directly to the recorded time-domain signals using a form of non-stationary convolution. The approach is computationally efficient and can be used with any detector geometry or reconstruction algorithm.
To prevent high-frequency noise from being amplified, the compensation is regularised using a Tukey window with a time-variant cutoff frequency. The cutoff frequency can be specified manually using the optional input 'FilterCutoff'
. This is set as a two-element vector corresponding to the cutoff frequency in Hz for the first and last time points, respectively. For a fixed cutoff, these should be specified as the same value, e.g., [3e6, 3e6]
. Alternatively, if 'FilterCutoff'
is set to 'auto'
(the default), the cutoff frequency is chosen based on the local time-frequency distribution of the recorded signals using the following steps:
- Compute the average time-frequency distribution of the input signals (the method can be defined using the optional input
'Distribution'
) - Threshold the time-frequency distribution to remove noise (the threshold value can be defined using the optional input
'NoiseThreshold'
) - Calculate the integral of the thresholded time-frequency distribution at each time point using
cumsum
- Find the cutoff frequency at each time point where the integral reaches a given percentage of the maximum value (this percentage can be defined using the optional input
'EnergyThreshold'
) - Increase the filter cutoff frequency by a fixed multiplier so the cutoff corresponds to the edge of the passband for the Tukey window (the multiplier can be defined using the optional input
'FrequencyMultiplier'
) - Smooth the variation of the cutoff frequency over time (the smoothing function can be defined using the optional input
'FitType'
) - Threshold any values of the cutoff frequency below zero or above the Nyquist limit
If the input contains a matrix of signals, the cutoff frequency is based on the average time frequency distribution. To calculate the cutoff frequency for each signal individually, this function should be called in a loop. This can be parallelised, for example, using parfor
from the parallel computing toolbox. For further details about this function and attenuation compensation using time variant filtering, see the reference below.
[1] B. E. Treeby (2013) "Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering," J. Biomed. Opt., vol. 18, no. 3, p.036008.
Inputs
signal |
matrix of time series to compensate indexed as (sensor_index, time_index) |
dt |
time step [s] |
c |
sound speed [m/s] |
alpha_0 |
power law absorption prefactor [dB/(MHz^y cm)] |
y |
power law absorption exponent [0 < y < 3, y ~= 1] |
Optional Inputs
Optional 'string', value pairs that may be used to modify the default computational settings.
Input | Valid Settings | Default | Description |
---|---|---|---|
'DisplayUpdates' |
(Boolean scalar) | true |
Boolean controlling whether command line updates and compute time are printed to the command line. |
'Distribution' |
'Rihaczek' 'Wigner' |
'Rihaczek' |
Time-frequency distribution used to automatically compute the filter cutoff frequency if 'FilterCutoff' is set to 'auto' . Note, the option 'Wigner' requires the Discrete TFD toolbox from http://tfd.sourceforge.net. |
'EnergyThreshold' |
(numeric scalar) | 0.98 |
Threshold value given as a percentage of the total amplitude spectrum used to choose the filter cutoff frequency at each time point. |
'FilterCutoff' |
(numeric two element vector) or'auto' |
'auto' |
Option to manually define the cutoff frequencies for a linear variation in the filter cutoff instead of using an automatic search. |
'FitType' |
'spline' 'linear' 'mav' |
'spline' |
Fitting type used to smooth the filter cutoff frequency after an automatic search, where 'spline' fits a smoothed spline, 'linear' fits a linear line, and 'mav' computes the moving average. |
'FrequencyMultiplier' |
(numeric scalar) | 2 |
By default, the compensation is regularised using a Tukey window with a time-variant cutoff frequency. The default Tukey window has a taper ratio of 0.5, so the filter cutoff frequency found by the automatic search is increased by a frequency multiplier so that the filter cutoff frequency corresponds to the edge of the passband of the Tukey window. |
'NumSplines' |
(numeric scalar) | 40 |
Number of spline segments used in the smoothing spline if 'FitType' is set to 'spline' . |
'NoiseThreshold' |
(numeric scalar) | 0.03 |
Threshold value given as a percentage of the signal maximum used to threshold the TFD before the automatic search for the filter cutoff. |
'Plot' |
(Boolean scalar) | false |
Boolean controlling whether a plot of the time frequency distribution and filter cutoff frequency are displayed. |
'PlotRange' |
(numeric two element vector) or'auto' 'full' |
'auto' |
Option to manually set the plot range in the frequency axis when 'Plot' is set to true . This can be manually specified, or set to 'auto' (1.5 x the maximum filter cutoff frequency) or 'full' (maximum supported frequency range). |
'TaperRatio' |
(numeric scalar) | 0.5 |
Taper ratio used to construct the Tukey Windows. |
'T0' |
(numeric scalar) | 0 |
Time index of T0 in the input signals. For photoacoustic imaging, T0 corresponds to the arrival of the excitation laser pulse at the sample. |
Outputs
signal_comp |
time series after attenuation compensation |
tfd |
average time frequency distribution of the input signals |
cutoff_freq |
filter cutoff frequency for each time index |