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

1. Compute the average time-frequency distribution of the input signals (the method can be defined using the optional input `'Distribution'`)
2. Threshold the time-frequency distribution to remove noise (the threshold value can be defined using the optional input `'NoiseThreshold'`)
3. Calculate the integral of the thresholded time-frequency distribution at each time point using `cumsum`
4. 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'`)
5. 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'`)
6. Smooth the variation of the cutoff frequency over time (the smoothing function can be defined using the optional input `'FitType'`)
7. 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

## Examples

- Attenuation Compensation Using Time Variant Filtering