Acoustic Imaging

Warning

This treatment is provided by an external module and is not part of the core system. If you’re interested in accessing or utilizing this feature, please contact us for further details and support.

Description

This treatment perform an acoustic imaging algorith to estimate the source location given the signal pressure of a set of observers. The treatment computes an acoustic imaging function \(f\) that should, ideally, have the following properties:

  • The function should be maximun at the source location

\[f(x = x_s) > f(x \neq x_s )\]

  • The function should measure the source strength

\[f(x = x_s) \propto q q^*\]

Where \(q\) is the source strength

(Treatment added in v2.6.0)

Parameters

  • base: Base

    Input base where each instant corresponds to different time step.

  • fwh_base: Base

    Input base formatted as the output base of the FWH treatment.

  • pressure_variable: str, default= 'pressure'

    The name of the pressure variable in the input base.

  • pref: float, default= 101325.0

    Reference pressure.

  • source_mesh Base

    The mesh in which compute the source strength.

  • obs_base: Base, default= None

    Base with the observers, if fwh_base or base contain already the observer positions, this keyword can be ignored.

  • build_obs_connectivity: bool, default= False

    If True and the observers data does not have a connectivity table, use a Delaunay triangulation to build one. The triangulation is only valid if all points are in X-Y plane.

  • nb_block: int

    The number of blocks used to compute the cross-spectral matrix.

  • overlap: float, default= 0.66

    Fraction of the signal length to overlap between consecutive blocks (e.g., 0.66 = 66%)

  • freq_min_max: list(float), default= None

    Truncate frequencies in CSM computation. The input is a list of 2 floats [min, max]. If either value is None, then there is no bound for that value. For example [100, None] will keep all frequencies above 100Hz, and [None 1000] will keep all frequecies below 1000Hz.

  • dt: float, default= 1.0

    The timestep of the input base.

  • time_t0: float, default= None

    If not None, clip the time series so it starts at time_t0.

  • time_tf: float, default= None

    If not None, clip the time series so it ends at time_tf.

  • time_duration: float, default= None

    If not None, this argument trims the time series so that its duration matches the specified value. Clipping is done from the end, so the resulting series ends at the same point as the original time series or time_tf.

  • csm_algorithm: str, default= 'welch'

    The algorithm to compute the CSM matrix. Possible values are:

    • 'welch': Uses the Welch method to estimate the CSM

    • 'multitaper': Uses the multitaper method to estimate the CSM.

    • 'precomputed': Uses an already precomputed CSM (see precomputed_csm).

  • precomputed_csm: str, default= None

    The path to the file contained the precomputed CSM.

  • params_multitaper: dict, default= {}

    Arguments to pass to the multitaper function of the form: {'arg_name': arg_value}

    • max_order: Maximum order of the Slepian sequence.

    • normalized_bandwith: Standardized half bandwidth of the window.

  • imaging_functions: python:callable or list(callable), default= beamforming.ConventionalFreqDomainBeamforming()

    The functions that compute the value of \(f\) at each source location. The function arguments must be the following:

    • A matrix of size \(N_{obs} \times N_{obs}\) containg the CSM matrix at a given frequency

    • The output of the transfer_function

    • The dtype in which the computation must be performed (either numpy.float32 or numpy.float64)

    • A flag to only return a string describing the name of the function

    If the last argument is True, the function must immediately return a string with its description. This string is used to name the variables in the output base.

    The function must return a numpy vector of size \(N_{s}\) containing the value of the function \(f\) at each source location.

    def img_func(CSM, tf, dtype, description):

    if description:
        return "Imaging description"

    # Compute f
    return f

t = antares.Treatment('acousticimaging')
t['imaging_functions'] = img_func

    This keyword can accept multiple imaging functions at the same time. See pairing keyword for a more detailed explanation.

  • transfer_functions: python:callable or list(callable), default= monopole_transfer_func()

    The functions that compute the transfer functions. The function arguments must be the following:

    • A matrix \(3 \times N_{src}\) with the coordinates of the source base.

    • A matrix \(3 \times N_{obs}\) with the coordinate of the observers.

    • A float value with the frequency of interest.

    • A flag to only return a string describing the name of the function

    If the last argument is True, the function must immediately return a string with its description. This string is used to name the variables in the output base.

    The function must return a numpy matrix of size \(N_{obs} \times N_{src}\) containing the value of the transfer function between the source and the observers location.

    def tf_func(src_coords, obs_coords, f, description):

    if description:
        return "Transfer function description"

    # Compute a
    return a

t = antares.Treatment('acousticimaging')
t['transfer_functions'] = tf_func

    This keyword can accept multiple imaging functions at the same time. See pairing keyword for a more detailed explanation.

  • pairing: str, default= 'combination'

    Defines how imaging and transfer functions are combined:

    • 'combination': each imaging function is paired with every transfer function (forming all possible pairs).

    • 'one_to_one': each imaging function is paired with the transfer function at the same position (first with first, second with second, etc.).

  • instant_name_format: str, default= '{freq:.0f}_Hz'

    A string used to name the instants in the output. You can include the {freq} placheholder in the string that will be replaced with the frequency value. For example "{freq:.2f}_Hz" will generate instant names with the frequency value with a two decimal values followed by an underscore and Hz (for example: '123.45_Hz').

    For more detailed information about string formatting you can consult: https://docs.python.org/3/tutorial/inputoutput.html#tut-f-strings

  • output_dir: str, default= '.'

    The output directory

  • output_psd_obs: bool, default= True

    If True, write in the output_dir the computed PSD for each observer. The base is saved in hdf_antares format with the name psd_obs

  • output_one_third_octave_bands: bool, default= False

    True if the treatment should compute the one third octave bands for the source and backpropaged base. The output base filenames are:

    • 'sources_one_third_octaves' for the source data

    • 'obs_one_third_octaves' for the backpropaged observers data

    Both files are saved in the output directory.

  • output_csm: bool, default= False

    If True, write in the output_dir the computed CSM to be reused later. The base is saved in hdf_antares format with the name CSM.h5

  • back_propagate: bool, default= False

    If True, the treatment uses the computed source strength to propagate the signal back to the observers. The backpropagated observers PSD is saved in the output directory. The files are saved in format hdf_antares and they are named 'psd_back_propag_*.h5.

  • dtype: str, default= 'float64'

    Float precision for numerical operations. If 'float64' all operations are done in double precision. If 'float32' are done in single precision.

  • verbose: bool, default= False

    If True, print the progress messages.

Preconditions

The treatment accepts two forms of data input:

  1. fwh_base: A base with the same format as the output of the FWH treatment.

  2. base: A more generic base, formatted as:

    • A single-zone base.

    • Containing as many instants as timesteps.

    • Each instant includes a pressure_variable with the pressure values for all observers.

    • Each instant contains an attribute called ‘Time’ with the physical time of the instant.

In addition, obs_base is not required if:

  • If base is used.

  • If fwh_base already includes observer positions.

The obs_base must follow the format:

  • A single-zone base.

  • Containing coordinates variables (e.g., x, y, z) that specifies the spatial position of each observer.

  • The number of observers defined in obs_base must match the number of pressure values given in fwh_base.

Postconditions

The treatment returns one base with the beamforming algorithm applied to the source_mesh. The base contains as many instants as frequencies computed. For each instant there are as many variables as combinations between transfer functions and imaging functions. Each variable name is structured as: "imaging_func_description - transfer_funtion_description".

The treatment can also output the PSD and backpropagated PSD for the observers studied (see output_psd_mics and back_propagate keywords). In each case, the output base contains as many instants as frequencies computed. Each instant contains the PSD value for each microphone. Additionally, for visualization purposes, if the microphones base does not have a connectivity dictionary, a Delaunay triangulation is performed to created a connected list of points.

Transfer functions

The beamforming algorithm decomposes the pressure signal for each observers as:

\[p(\vec{x}) = a(\vec{x}, \vec{x_0}, \vec{x_s}) q(\vec{x_s})\]

where \(\vec{x}\) is the observers location, \(\vec{x_s}\) is the source location, \(\vec{x_0}\) is the reference location, \(a\) is the transfer matrix and \(q\) is the source strength. This transfer matrix is key in the algorithm to compute the steering vectors of the beamforming algorithm.

The treatment allows the user to input its own transfer function to the treatment using the keyword transfer_functions

The function must be a regular python function with three parameters in the following order:

  • The observers coordinates as a \(3 \times N_{obs}\) array

  • The source coordinates as a \(3 \times N_{s}\) array

  • The current frequency in Hz.

  • A boolean description_only to return only the description of the function.

The function must return 1 value:

  • If description is set to True, it must return a string with the description of the transfer function. This string will be used to name the output variables.

  • Otherwise, the function must return the computed transfer function as an array \(N_{obs} \times N_{s}\)

Antares implements already implements a simple monopole transfer function:

\[a\left(\vec{x}_{obs}, \vec{x}_0, \mathbf{x}_s\right)=\frac{r_{s, 0}}{r_{s, obs}} e^{-j k\left(r_{s, obs}-r_{s, 0}\right)}\]

Where \(r\) is the distance between the two points

This function can be used as:

monopole = antares.plugins.acoustics.utils.monopole_transfer_func(c=340.0, mach=np.array([0, 0, 0]))

t = antares.Treatment('acousticimaging')
t['tranfer_function'] = monopole

Imaging algorithm

The imaging algorithm computes the function \(f\) on the source mesh. The value of this keyword must be a function that accepts the following parameters:

  • The cross spectral matrix as an array of \(N_{obs} \times N_s\).

  • The output of the transfer function.

  • The dtype in which the calculation are performed (either np.float32 or np.float64).

  • a boolean that if it is set to True the function must only output its description.

If the last bool parameters is False, the function must output the value of \(f\) as an array of length \(N_s\).

Antares implements the beamforming imaging algorithm (Sarradj 2012.).

\[f(x) = \vec{h}^{H} G \vec{h}\]

Where \(G\) is the cross spectral matrix and \(\vec{h}\) are the called steering vectors. Different formulation for these steering vector are available in the literature:

  • Formulation I

\[h_i^{\mathrm{I}}=\frac{1}{N} \frac{a_i\left(\vec{x}_0, \vec{x}_t\right)}{\left|a_i\left(\vec{x}_0, \vec{x}_t\right)\right|}\]

  • Formulation II

\[h_i^{\mathrm{II}}=\frac{1}{N} \frac{a_i\left(\vec{x}_0, \vec{x}_t\right)}{a_i\left(\vec{x}_0, \vec{x}_t\right) a_i^*\left(\vec{x}_0, \vec{x}_t\right)}\]

  • Formulation III

\[h_i^{\mathrm{III}}=\frac{a_i\left(\vec{x}_0, \vec{x}_t\right)}{\mathbf{a}^H\left(\vec{x}_0, \vec{x}_t\right) \mathbf{a}\left(\vec{x}_0, \vec{x}_t\right)}\]

  • Formulation IV

\[h_i^{\mathrm{IV}}=\frac{1}{\sqrt{N}} \frac{a_i\left(\vec{x}_0, \vec{x}_t\right)}{\sqrt{\mathbf{a}^H\left(\vec{x}_0, \vec{x}_t\right) \mathbf{a}\left(\vec{x}_0, \vec{x}_t\right)}}\]

To use these steering vectors:

imaging_func = antares.plugins.acoustics.utils. ConventionalFreqDomainBeamforming(formulation="III")
t = antares.Treatment('acousticimaging')
t['imaging_functions'] = imaging_func

Example

The following example illustrates the application of the acoustic imaging algorithm to the output of the FWH treatment. The pressure field was propagated to an array of observers (shown as blue points). The source location is then estimated using the Beamforming algorithm.

In this case, the original FWH surface is a sphere surrounding a single monopole, while the source mesh is defined as a plane intersecting the center of the sphere.

#!/usr/bin/env python3

import antares
import os

from antares.plugins.acoustics.utils import monopole, beamforming

output_folder = os.path.join("OUTPUT", "TreatmentAcousticImaging")
os.makedirs(output_folder, exist_ok=True)

data_folder = os.path.join('..', 'data', 'ACOUSTIC_IMAGING')

## Read input bases
## ================
obs_file = os.path.join(data_folder, 'micros.txt')
obs_base = antares.io.Read(filename=obs_file, format='column')

fwh_file = os.path.join(data_folder, 'fwh.h5')
fwh_base = antares.io.Read(filename=fwh_file, format='hdf_antares')

source_file = os.path.join(data_folder, 'source.cgns')
source_base = antares.io.Read(filename=source_file, format='hdf_cgns')

## Run Acoustic Imaging treatment
## ==============================

monopole_transfer_func = monopole.monopole_transfer_func()
imaging_func = beamforming.ConventionalFreqDomainBeamforming(formulation="III")

img = antares.Treatment('acousticimaging')
img['fwh_base'] = fwh_base
img['source_mesh'] = source_base
img['obs_base'] = obs_base
img['pressure_variable'] = 'pressure'
img['dtype'] = 'float64'
img['pref'] = 0. # Data is already centered around reference pressure
img['freq_min_max'] = [500, 6e3] # [Hz], frequency range for the imaging analysis
img['nb_block'] = 8 # number of PSD blocks
img['overlap'] = 0.66  # fraction of overlap
img['transfer_functions'] = monopole_transfer_func
img['imaging_functions'] = imaging_func
img['instant_name_format'] = "{freq_int}_Hz"
img['output_psd_obs'] = True
img['output_dir'] = output_folder
img['verbose'] = True

base_img = img.execute()

## Plot resutls for 1st frequency
## ==============================

surf_fwh = antares.io.Read(filename=os.path.join(data_folder, 'fwh_surface.h5'), format='hdf_antares')
surf_fwh.attrs['plotbases'] = {
    'color_by': 'zones',
    'zone_colors': 'grey',
    'transparency': 0.85,
    'colorbar': False,
}

base_img.attrs['plotbases'] = {
    'color_by': 'Conventional Beamforming - Monopole (no flow)',
    'colorbar_horizontal': True,
    'colorbar_title_fontsize': 56,
    'colorbar_ticks_fontsize': 56,
    'colorbar_labelformat': '%.2f'
}

obs_base.attrs['plotbases'] = {
    'color_by': 'zones',
    'zone_colors': 'blue',
    'colorbar': False,
    'point_radius': 0.04,
}

first_freq = base_img[0].keys()[0]
title = {
    'text': f"Acoustic imaging for frequency: {first_freq}",
    'position': (0.3, 0.98),
    'fontsize': 46,
    'bold': True,
  }

plot = antares.Treatment('plotbases')
plot['bases'] = [obs_base, base_img[:, (0,), :], surf_fwh]
plot['show'] = False
plot['color_by'] = 'in_attr'
plot['figsize'] = os.path.join(data_folder, 'camera.txt')
plot['texts'] = [title]
plot['camera_save_settings'] = os.path.join(data_folder, 'camera.txt')
plot['camera_load_settings'] = os.path.join(data_folder, 'camera.txt')
plot['output_file'] = os.path.join(output_folder, 'results.jpg')

plot.execute()
../../../../_images/result.jpg

References

Sarradj E. (2012). Three-Dimensional Acoustic Source Mapping with Different Beamforming Steering Vector Formulations. Advances in Acoustics and Vibration, 2012. (link)