Added room acoustics test statistics analysis scripts + related packages legal notices.

# Statistical analysis of IVAS Reverb Characterization testing

This project is used to analyze the result logs for the multiple renderer preference testing carried out as part of the IVAS characterization testing.

# Prerequisites 

## Data Structure

The code assumes that all results files are stored in same folder, and conform to the specification for ARL Step result files. That is, a tab separated list, with the header of: 

```

Lab	Listener	Session	File	Sig	Sys	Score 

```

The folder must also contain the `RoomAcoustics.asi` file, detailing for each test item the files used for each renderer condition. This conforms to the `mushra` style session file of ARL Step.

## Python

The project uses `uv` ro manage the python dependencies. `uv` can be installed with your favorite package manager. For more information, have a look at

[the installation instructions](https://docs.astral.sh/uv/getting-started/installation/). With python already installed on your system it should be sufficient to run `pip install uv` from a command prompt.

# Running the code

With UV installed it is sufficient to call the function with

```

uv run statistical_analysis.py <path-to-results-directory>

```

On the first instance of running `uv` it will create a virtual environment and populate it with the requirements as specified in `pyproject.toml`. If you make changes to the code which change the dependencies use `uv sync` to update your virtual environment.

Teo optional arguments can be given `--show_plots` and `--store_plots` to either show the plots in the web browser, or to store the plots as PNG files into the results folder. See the help for more information:

```

> uv run statistical_analysis.py --help

usage: statistical_analysis.py [-h] [--show_plots] [--save_plots] results_folder

positional arguments:

  results_folder  Path to the directory containing the results files

options:

  -h, --help      show this help message and exit

  --show_plots    Shows the plots in the web browser.

  --save_plots    Stores the plots into the results folder, note: requires Chrome / Chromeium to be installed on the local machine.

```

# Interpreting the output

The code optionally generates two graphics, one showing the results per signal type, and one averaged across all signals. For the example data contained in `./example_data/` these are:

![Per Item](./example_data/per-item.png "Data per signal type")

![across items](./example_data/per-renderer.png "Data across signal types")

In the plot averaged across all signals the markers above the plots show pairwise comparisons that are considered statistically significant (p < 0.05).

The code performs two startisical analysis, firstly a 2-way ANOVA, and then a Tukey Honest Squares Difference (HSD) post-hoc test. The output of both is stored into `stats.csv` in the input results folder.

## ANOVA

Three effects are analysed, `Signal`, `Renderer`, and the interaction of the two. The results for the example data are shown below:

|   |sum_sq|df|F|PR(>F)|

|---|---|---|---|---|

|C(Signal)|1234.5416666666592|11.0|0.7687058945620543|0.6706403030269787|

|C(Renderer)|154214.0416666672|3.0|352.08685312024477|4.4697831679421755e-66|

|C(Signal):C(Renderer)|3049.083333333332|33.0|0.6328524975785247|0.9373952536868595|

|Residual|21024.0|144.0|||

In this data `Signal` is considered to be non-significant (p>0.05 and a small F-value), while the effect of `Renderer` is highly significant (p<0.05 and a large F-value). There is no interaction effect of `signal` and `renderer`. This means that all of the variance of the data is assumed to come from the `Renderer` only. To understand the interaction between the different renderers a post-hoc test must be performed.

## Tukey HSD post-hoc

Tukey Honest Squares Difference post-hoc test makes pairwise comparisons between all groups, and determines whether the mean difference is significantly different from zero. If the mean difference between groups is different from zero then you can reject the null hypothesis, meaning that that pair are considered to be different from each other.

In the example data given below the first 3 comparisons are considered to be significant, while the second three are non-significant. This translates to all other renderers being different from `IVAS_Binaural`, and not different from each other. This is shown by the markers in the plot above.

|         group1        |             group2             |meandiff|p-adj | lower | upper |reject|

|---|---|---|---|---|---|---|

|IVAS_Binaural          |IVAS_BinauralRoomIR             |66.8958 |0.0   |60.7568|73.0349|True  |

|IVAS_Binaural          |IVAS_BinauralRoomReverb         |64.625  |0.0   |58.4859|70.7641|True  |

|IVAS_Binaural          |ReferenceRenderer_BinauralRoomIR|64.7292 |0.0   |58.5901|70.8682|True  |

|IVAS_BinauralRoomIR    |IVAS_BinauralRoomReverb         |-2.2708 |0.7729|-8.4099|3.8682 |False |

|IVAS_BinauralRoomIR    |ReferenceRenderer_BinauralRoomIR|-2.1667 |0.797 |-8.3057|3.9724 |False |

|IVAS_BinauralRoomReverb|ReferenceRenderer_BinauralRoomIR|0.1042  |1.0   |-6.0349|6.2432 |False |

 No newline at end of file

import os

import pandas as pd

import plotly.express as px

import statsmodels.api as sm

from statsmodels.formula.api import ols

from statsmodels.stats.multicomp import pairwise_tukeyhsd

import sys

import numpy as np

from operator import itemgetter

import argparse

''' 

******************************************************************************************************

   (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,

   Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,

   Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,

   Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other

   contributors to this repository. All Rights Reserved.

   This software is protected by copyright law and by international treaties.

   The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,

   Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,

   Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,

   Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other

   contributors to this repository retain full ownership rights in their respective contributions in

   the software. This notice grants no license of any kind, including but not limited to patent

   license, nor is any license granted by implication, estoppel or otherwise.

   Contributors are required to enter into the IVAS codec Public Collaboration agreement before making

   contributions.

   This software is provided "AS IS", without any express or implied warranties. The software is in the

   development stage. It is intended exclusively for experts who have experience with such software and

   solely for the purpose of inspection. All implied warranties of non-infringement, merchantability

   and fitness for a particular purpose are hereby disclaimed and excluded.

   Any dispute, controversy or claim arising under or in relation to providing this software shall be

   submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in

   accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and

   the United Nations Convention on Contracts on the International Sales of Goods.

*******************************************************************************************************

'''

df_format = ({

        'Lab':[],

        'Listener':[],

        'session_file':[],

        'Signal':[],

        'Renderer':[],

        'Score':[],

})

def read_result_file(results_file, test_item_labels,wav_file_prefixes):

    ''' Read a .txt test result file in the format used by ARL Step.

    Parameters:

    ----------

    results_file: string or path

        Path to the file to be read

    test_item_labels: list[string]

        List of Test Item Labels as returned by read_asi_file()

    wav_file_prefixes: list[string]

        Liat of wav file prefixes as returned by read_asi_file()

    Returns:

    ----------

    results: Pandas DataFrame

        pd DataFrame with the columns defined in df_format.

    '''

    df = pd.DataFrame(df_format)

    try:

        with open(results_file,'r') as res_file:

            for line in res_file:

                try:

                    if not line.startswith('Lab'):

                        split_line = line.split('\t')

                        cond_idx = [idx for idx,s in enumerate(wav_file_prefixes) if s in split_line[3]]

                        df.loc[len(df)] = [split_line[0], 

                                        f"{split_line[0]}_{split_line[1]}",

                                        split_line[2],

                                        test_item_labels[cond_idx[0]],

                                        split_line[4],

                                        int(split_line[6])

                                        ]

                except Exception as e2:

                    print(f"Cannot parse line {line}")

    except Exception as e:

        raise ValueError(f"Unable to import file {results_file}: {e}")

    return df

def read_asi_file(asi_file):

    ''' Read a .asi file in the format used by ARL Step.

    Parameters:

    ----------

    asi_file: string or path

        Path to the .asi file to be read

    Returns:

    ----------

    test_item_labels: list[strings]

        List of strings with the Name of each test item (audio stream)

    wav_file_prefixes: list[strings]

        List of prefixes used for each wav file associated with the test items of test_item_labels.

    '''

    with open(asi_file,'r') as af:

        lines = af.readlines()

        label_idxs = [idx for idx, s in enumerate(lines) if s.startswith('#')]

        test_item_labels = []

        wav_file_prefxies = []

        for i in range(len(label_idxs)):

            test_item_labels.append(lines[label_idxs[i]].rstrip().replace('# ','').replace('-',' '))

            if i is not len(label_idxs)-1:

                wav_file_prefxies.append(os.path.commonprefix(lines[label_idxs[i]+1:label_idxs[i+1]]))

            else:

                wav_file_prefxies.append(os.path.commonprefix(lines[label_idxs[i]+1:]))

    return test_item_labels, wav_file_prefxies

def add_p_value_annotation(fig, array_columns, _format=dict(interline=0.03, text_height=1.02, color='black')):

    ''' Adds notations for significance between two pairs in a boxplot

    Based on code taken from https://stackoverflow.com/questions/67505252/plotly-box-p-value-significant-annotation

    Parameters:

    ----------

    fig: figure

        plotly boxplot figure

    array_columns: np.array

        array of which columns to mark as significant 

        e.g.: [[0,1], [1,2]] marks column 0 and 1 as significant and 1 and 2 as significant (with an *)

    _format: dict

        format characteristics for the lines

    Returns:

    -------

    fig: figure

        figure with the added notation

    '''

    # Specify in what y_range to plot for each pair of columns

    y_range = np.zeros([len(array_columns), 2])

    for i in range(len(array_columns)):

        y_range[i] = [1.01+i*_format['interline'], 1.02+i*_format['interline']]

    # Get values from figure

    fig_dict = fig.to_dict()

    # Print the p-values

    for index, column_pair in enumerate(array_columns):

        symbol = '*'

        # Vertical line

        fig.add_shape(type="line",

            xref="x", yref="y domain",

            x0=column_pair[0], y0=y_range[index][0], 

            x1=column_pair[0], y1=y_range[index][1],

            line=dict(color=_format['color'], width=2,)

        )

        # Horizontal line

        fig.add_shape(type="line",

            xref="x", yref="y domain",

            x0=column_pair[0], y0=y_range[index][1], 

            x1=column_pair[1], y1=y_range[index][1],

            line=dict(color=_format['color'], width=2,)

        )

        # Vertical line

        fig.add_shape(type="line",

            xref="x", yref="y domain",

            x0=column_pair[1], y0=y_range[index][0], 

            x1=column_pair[1], y1=y_range[index][1],

            line=dict(color=_format['color'], width=2,)

        )

        ## add text at the correct x, y coordinates

        ## for bars, there is a direct mapping from the bar number to 0, 1, 2...

        fig.add_annotation(dict(font=dict(color=_format['color'],size=14),

            x=(column_pair[0] + column_pair[1])/2,

            y=y_range[index][1]*_format['text_height'],

            showarrow=False,

            text=symbol,

            textangle=0,

            xref="x",

            yref="y domain"

        ))

    return fig

def plot_results(results, HSD_results, args):

    output_folder = args.results_folder

    # y-ticks

    yticks = [10, 30, 50, 70, 90]

    ylabels = ['Very unpreferred', 'Unpreferred', 'Neutral', 'Preferred', 'Very preferred']

    fig = px.box(results, x="Signal", y="Score",color='Renderer')

    for n in range(results['Signal'].nunique()-1):

        fig.add_vline(x=n+0.5,line_dash='dash', line_color='gray',)

    fig.update_xaxes(tickangle=-45)

    fig.update_yaxes(tickvals=yticks, ticktext=ylabels)

    fig.update_layout(height=800, width=1200)

    fig2 = px.box(results, x="Renderer", y="Score", color='Renderer')

    categories = fig2._layout_obj.xaxis.categoryarray

    sig_comparisons = []

    for comparison in HSD_results.data:

        if comparison[6] is np.True_:

            sig_comparisons.append([categories.index(comparison[0]),categories.index(comparison[1])])

            sig_comparisons[-1].sort()

    sig_comparisons = sorted(sig_comparisons, key=itemgetter(0,1))

    fig2 = add_p_value_annotation(fig2,sig_comparisons)

    fig2.update_yaxes(tickvals=yticks, ticktext=ylabels)

    fig2.update_layout(showlegend=False, margin=dict(l=20, r=20, t=40*len(sig_comparisons), b=20),height=600,width=1200)

    if args.show_plots:

        fig.show()

        fig2.show()

    if args.save_plots:

        fig.write_image(os.path.join(output_folder, 'per-item.png'),height=800,width=1200)

        fig2.write_image(os.path.join(output_folder, 'per-renderer.png'),height=600,width=1200)

def main(args):

    ''' Performs 2-way ANOVA and Tukey HSD on the results of a listening test.

    Assumes that the test is a MuSCR style test, (i.e. MUSHRA without a reference)

    Parameters:

    ----------

    args: 

        args.results_folder: String containing path to folder containing the results files in .txt format and .asi test control file

        args.show_plots: Boolean, whether to show the plots in the web browser or not

        args.save_plots: Boolean, whether to save the plots into the results_folder as .png files.

    Returns:

    -------

    None

    '''

    results_folder = args.results_folder 

    files = [file for file in os.listdir(results_folder) if file.endswith('.txt')]

    asi_file = os.path.join(results_folder, 'RoomAcoustics.asi')

    test_item_labels, wav_file_prefixes =  read_asi_file(asi_file)

    results = pd.DataFrame(df_format)

    for file in files:

        results = pd.concat([results, read_result_file(results_file = os.path.join(results_folder,file),test_item_labels=test_item_labels, wav_file_prefixes=wav_file_prefixes)])

    output_file = os.path.join(results_folder,'stats.csv')

    model = ols('Score ~ C(Signal) + C(Renderer) + C(Signal):C(Renderer)', data=results).fit()

    anova = sm.stats.anova_lm(model, typ=2)

    anova.to_csv(output_file)

    print("===============================================================================================")

    print("                                          2-way ANOVA                                          ")

    print("===============================================================================================")

    print(anova)

    print("-----------------------------------------------------------------------------------------------")

    tukey = pairwise_tukeyhsd(endog=results['Score'], groups=results['Renderer'], alpha=0.05)

    with open(output_file, 'a') as op:

        op.write('\n')

        op.write(tukey.summary().as_csv()) 

    print('\n')

    print("===============================================================================================")

    print(tukey)

    if args.show_plots or args.save_plots:

        plot_results(results, tukey.summary(), args)

if __name__ == "__main__":

    parser = argparse.ArgumentParser()

    parser.add_argument("results_folder", help="Path to the directory containing the results files", type=str)

    parser.add_argument("--show_plots", action="store_true", help="Shows the plots in the web browser.")

    parser.add_argument("--save_plots", action="store_true", help="Stores the plots into the results folder, note: requires Chrome / Chromeium to be installed on the local machine.")

    args = parser.parse_args()

    if not os.path.isdir(args.results_folder):

        raise IsADirectoryError(f"Given path is not a valid directory.")

    main(args)

 No newline at end of file

Applies to Numpy

Package: numpy

Version: 1.26.4

License: BSD 3-Clause License

Copyright (c) 2005-2022, NumPy Developers.

Copyright (c) 2005–2025, NumPy Developers.

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are

met:

    * Redistributions of source code must retain the above copyright

       notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above

       copyright notice, this list of conditions and the following

       disclaimer in the documentation and/or other materials provided

       with the distribution.

    * Neither the name of the NumPy Developers nor the names of any

       contributors may be used to endorse or promote products derived

       from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR

A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT

OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,

SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT

LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,

DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY

THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT

(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice,

  this list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,

  this list of conditions and the following disclaimer in the documentation

  and/or other materials provided with the distribution.

* Neither the name of the NumPy Developers nor the names of any contributors

  may be used to endorse or promote products derived from this software

  without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE

DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE

FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,

OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Additional licenses from https://github.com/numpy/numpy/blob/main/LICENSES_bundled.txt:

The NumPy repository and source distributions bundle several libraries that are

Package: pandas

Version: 2.3.1

License: BSD License

BSD 3-Clause License

Copyright (c) 2008-2011, AQR Capital Management, LLC, Lambda Foundry, Inc. and PyData Development Team

All rights reserved.

Copyright (c) 2011-2025, Open source contributors.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this

  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,

  this list of conditions and the following disclaimer in the documentation

  and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names of its

  contributors may be used to endorse or promote products derived from

  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE

DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE

FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,

OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 No newline at end of file

Package: plotly

Version: 6.2.0

License: MIT License

Copyright (c) 2019, Plotly, Inc.

Permission is hereby granted, free of charge, to any person obtaining a copy

of this software and associated documentation files (the “Software”), to deal

in the Software without restriction, including without limitation the rights

to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

copies of the Software, and to permit persons to whom the Software is

furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all

copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE

SOFTWARE.

 No newline at end of file