# 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 |