diff --git a/characterization/room_acoustics/binaural_sdm/BinauralSDM_for_IVAS.m b/characterization/room_acoustics/binaural_sdm/BinauralSDM_for_IVAS.m new file mode 100644 index 0000000000000000000000000000000000000000..f871757f251432bf26d9a80a668115571f1d45fe --- /dev/null +++ b/characterization/room_acoustics/binaural_sdm/BinauralSDM_for_IVAS.m @@ -0,0 +1,416 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Description + +% This script is based on an example as provided in the BinauralSDM +% repository (https://github.com/facebookresearch/BinauralSDM). +% Copyright (c) Facebook, Inc. and its affiliates. +% The BRIRs generated by this script include the steps described in +% Amengual et al. 2020 - DOA Postprocessing (smoothing and quantization) +% and RTMod+AP equalization. + +% The overall process is as follows: +% - Spatial data from a multichannel RIR is obtained using the SDM +% utilizing the functions from the SDM Toolbox (Tervo & Patynen). +% - The spatial information is smoothed and quantized as proposed in +% Amengual et al. 2020. +% - Preliminary BRIRs are synthesized by selecting the closest directions +% from an HRIR dataset. In this case, Orange HRIR data set #53 is used. +% +% - The reverberation time of the preliminary BRIRs is corrected using the +% RTMod method introduced in Amengual et al. 2020. +% - After RTMod correction, a cascade of 3 all-pass filters is applied to +% both left and right channels to improve the diffuseness of the late +% reverberation tail. This process is also introduced in Amengual et al. +% 2020. +% - Responses for arbitrarily defined head orientations are generated +% using the previous steps by rotating the DOA vectors. +% - The data is saved in a user defined folder. +% +% +% References: +% - (Tervo et al. 2013) - "Spatial Decomposition Method for Room Impulse +% Responses", JAES, 2013. +% - (SDM Toolbox) - "SDM Toolbox", S. Tervo and J. Patynen, Mathworks +% Exchange +% - (Amengual et al. 2020) - "Optimizations of the Spatial Decomposition +% Method for Binaural Reproduction", JAES 2020. +% Author: Sebastia Amengual (samengual@fb.com) +% +% Adapted: Marek Szczerba (marek.szczerba@philips.com) (MS) +% Last modified: 08/25/2025 +% +% HOW TO +% - Clone the BinauralSDM repository +% (https://github.com/facebookresearch/BinauralSDM). +% - Place this file in the ./BinauralSDM/Src/Examples folder. +% - Update the settings on top of the BinauralSDM_for_IVAS function to +% reflect the actual folder structure. +% - Run this script. +% +% NOTES (MS) +% Eigenmike geometry is supported out-of-the-box. Check create_MicGeometry +% in BinauralSDM/Src folder + + +function BinauralSDM_for_IVAS + + % BinauralSDM_for_IVAS, Binaural SDM processing for IVAS + % + % refactored into functional programming with local variables + + % Parameters and settings + % RIR + RIR_input_path = '..\..\..\RIRs\'; + RIR_room_name = 'Auditorium'; + RIR_source_pos = 'S1'; + + % HRIR + HRIR_input_path = '..\..\..\ivas-codec\scripts\binauralRenderer_interface\HRIRs_sofa\'; + HRIR_Filename = 'HRIR_128_Meth5_IRC_53_Q10_symL_Itrp1_48000.sofa'; + HRIR_Subject = 'HRIR-53'; % Name of the HRIR subject (only used for naming purposes while saving). + HRIR_Type = 'SOFA'; % File format of the HRIR. Only SOFA is supported for now. + + % BRIR + BRIR_output_path='.\'; + BRIR_angles_mode='9.1.6'; % 7.1.4 or 9.1.6 or All + + % settings is a struct that summarizes all releveant settings for processing files + settings.RIR_database_input_path = fullfile(RIR_input_path); + settings.RIR_room = RIR_room_name; + settings.RIR_SourcePos = RIR_source_pos; + + settings.HRIR_input_path = fullfile(HRIR_input_path); + settings.HRIR_filename = HRIR_Filename; + settings.HRIR_subject = HRIR_Subject; + settings.HRIR_type = HRIR_Type; + + settings.BRIR_output_path = fullfile(BRIR_output_path); + settings.BRIR_AnglesMode = BRIR_angles_mode; + + % setup processing + time_start = tic(); + SRIR_data = init_SRIR_data(settings); + SDM_Struct = init_SDM_Struct(SRIR_data); + BRIR_data = init_BRIR_data(settings, SRIR_data); + Plot_data = init_Plot_data(SRIR_data, BRIR_data); + + % do processing + SRIR_data = analyze(SRIR_data, SDM_Struct, Plot_data); + synthesize(SRIR_data, BRIR_data, Plot_data); + + % done + time_exec = toc(time_start); + fprintf('\n... completed in %.0fh %.0fm %.0fs.\n', time_exec/60/60, time_exec/60, mod(time_exec, 60)); +end + +function SRIR_data = init_SRIR_data(settings) + % Analysis parameters + MicArray = 'Eigenmike'; % FRL Array is 7 mics, 10cm diameter, with central sensor. Supported geometries are FRL_10cm, FRL_5cm, + % Tetramic and Eigenmike. Modify the file create_MicGeometry to add other geometries (or contact us, we'd be happy to help). + Room = settings.RIR_room; % Name of the room. RIR file name must follow the convention RoomName_SX_RX.wav + SourcePos = settings.RIR_SourcePos; % Source Position. RIR file name must follow the convention RoomName_SX_RX.wav + ReceiverPos = 'R1'; % Receiver Position. RIR file name must follow the convention RoomName_SX_RX.wav + Database_Path = settings.RIR_database_input_path; % Relative path to folder containing the multichannel RIR + fs = 48e3; % Sampling Rate (in Hz). Only 48 kHz is recommended. Other sampling rates have not been tested. + MixingTime = 0.08; % Mixing time (in seconds) of the room for rendering. Data after the mixing time will be rendered + % as a single direction independent reverb tail and AP rendering will be applied. + DOASmooth = 16; % Window length (in samples) for smoothing of DOA information. 16 samples is a good compromise for noise + % reduction and time resolution. + DOAOnsetLength = 128; % Length (in samples) for assignment of a constant (averaged) DOA for the onset / direct sound. + DenoiseFlag = true; % Flag to perform noise floor compensation on the multichannel RIR. If set, it ensures that the RIR decays + % progressively and removes rendering artifacts due to high noise floor in the RIR. + PlotDenoisedRIR = false; + FilterRawFlag = true; % Flag to perform band pass filtering on the multichannel RIR prior to DOA estimation. If set, only + % information between 200 Hz and 8 kHz (by default) will be used for DOA estimation. This helps increasing + % robustness of the estimation. See create_BRIR_data.m for customization of the filtering. + AlignDOAFlag = true; % If this flag is set, the DOA data will be rotated so the direct sound is aligned to 0,0 (az, el). + + % Initialize SRIR data struct + SRIR_data = create_SRIR_data('MicArray',MicArray,... + 'Room',Room,... + 'SourcePos',SourcePos,... + 'ReceiverPos',ReceiverPos,... + 'Database_Path',Database_Path,... + 'fs',fs,... + 'MixingTime',MixingTime,... + 'DOASmooth',DOASmooth,... + 'DOAOnsetLength',DOAOnsetLength,... + 'Denoise',DenoiseFlag,... + 'PlotDenoisedRIR',PlotDenoisedRIR,... + 'FilterRaw',FilterRawFlag,... + 'AlignDOA',AlignDOAFlag); +end + + +function SDM_Struct = init_SDM_Struct(SRIR_data) + SpeedSound = 345; % Speed of sound in m/s (for SDM Toolbox DOA analysis) + WinLen = 24; %62; % Window Length (in samples) for SDM DOA analysis. For fs = 48kHz, sizes between 36 and 64 seem appropriate. + % The optimal size might be room dependent. See Tervo et al. 2013 and Amengual et al. 2020 for a discussion. + + % Initialize SDM analysis struct (from SDM Toolbox) + SDM_Struct = createSDMStruct('c',SpeedSound,... + 'fs',SRIR_data.fs,... + 'micLocs',SRIR_data.ArrayGeometry,... + 'winLen',WinLen); +end + + +function BRIR_data = init_BRIR_data(settings, SRIR_data) + %% Rendering parameters + DOAAzOffset = 0.0; % Azimuth rotation in degrees aplied after DOA estmation and before DOA quantization / BRIR rendering. + DOAElOffset = 0.0; % Elevation rotation in degrees aplied after DOA estmation and before DOA quantization / BRIR rendering. + QuantizeDOAFlag = true; % Flag to determine if DOA information must me quantized. + DOADirections = 50; % Number of directions to which the spatial information will be quantized using a Lebedev grid. + BandsPerOctave = 1; % Bands per octave used for RT60 equalization + EqTxx = 20; % Txx used for RT60 equalization. For very small rooms or low SNR measurements, T20 is recommended. Otherwise, T30 is recommended. + RTModRegFreq = false; % Regularization frequncy in Hz above which the RTmod modification of the late reverberation should be regularized. + NamingCond = sprintf('Quantized%dDOA', DOADirections); % String used for naming purposes, useful when rendering variations of the sasme RIR. + BRIR_Length = 0.0; % Length of BRIR in seconds, will be chosen by analysis of the room reverberation time if unspecified. + BRIR_Atten = 24; % Attenuation of the rendered BRIRs (in dB). Useful to avoid clipping. Use the same value when rendering various + % positions in the same room to maintain relative level differences. + AnglesMode = settings.BRIR_AnglesMode; % Angle processing mode: '7.1.4', '9.1.6' or 'All' + + % AzOrient = (-180:2:180)'; % Render BRIRs every 2 degrees in azimuth. + % ElOrient = (-90:5:90)'; % Render BRIRs every 5 degrees in elevation. + switch AnglesMode + case 'All' + AzOrient = (-180:2:180)'; % Render BRIRs every 2 degrees in azimuth. + ElOrient = (-90:6:90)'; % Render BRIRs every 6 degrees in elevation. + case '7.1.4' % 7.1+4 loudspeaker setup + AzOrient = [ 30.0, -30.0, 0.0, 135.0, -135.0, 90.0, -90.0, 30.0, -30.0, 135.0, -135.0]; + ElOrient = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 35.0, 35.0, 35.0, 35.0]; + case '9.1.6' % 9.1+6 loudspeaker setup + AzOrient = [ 30.0, -30.0, 0.0, 135.0, -135.0, 110.0, -110.0, 90.0, -90.0, 30.0, -30.0, 110.0, -110.0, 135.0, -135.0]; + ElOrient = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0]; + end + + ExportSofaFlag = true; % Flag to determine if BRIRs should be exported in a SOFA-file. + ExportWavFlag = false; % Flag to determine if BRIRs should be exported in indivudal WAV-files. + ExportDSERcFlag = true; % Flag to determine if BRIRs should be exported with direct sound and early reflections combined (WAV-files only). + ExportDSERsFlag = true; % Flag to determine if BRIRs should be exported with direct sound and early reflections separated (WAV-files only). + DestinationPath = settings.BRIR_output_path; % Folder where the resulting BRIRs will be saved. + + [~,~] = mkdir(settings.HRIR_folder); % ignore warning if directory already exists + HRIR_Path = fullfile(settings.HRIR_folder, settings.HRIR_filename); + + % Append configuration to destination path + DestinationPath = fullfile(DestinationPath, strrep(settings.HRIR_subject, ' ', '_'), ... + sprintf('%s_%s_%s', SRIR_data.Room, SRIR_data.SourcePos, SRIR_data.ReceiverPos), ... + NamingCond); + + % Initialize re-synthesized BRIR struct + BRIR_data = create_BRIR_data('MixingTime',SRIR_data.MixingTime,... + 'HRTF_Subject',settings.HRIR_Subject,... + 'HRTF_Type',settings.HRIR_Type,... + 'HRTF_Path',HRIR_Path,... + 'BandsPerOctave',BandsPerOctave,... + 'EqTxx',EqTxx,... + 'RTModRegFreq',RTModRegFreq,... + 'Length',BRIR_Length,... + 'RenderingCondition',NamingCond,... + 'Attenuation',BRIR_Atten,... + 'AzOrient',AzOrient,... + 'ElOrient',ElOrient,... + 'DOAAzOffset',DOAAzOffset,... + 'DOAElOffset',DOAElOffset,... + 'QuantizeDOAFlag',QuantizeDOAFlag,... + 'DOADirections',DOADirections,... + 'ExportSofaFlag',ExportSofaFlag,... + 'ExportWavFlag',ExportWavFlag,... + 'ExportDSERcFlag',ExportDSERcFlag,... + 'ExportDSERsFlag',ExportDSERsFlag,... + 'DestinationPath',DestinationPath,... + 'fs',SRIR_data.fs); + + if strcmp(AnglesMode, '7.1.4') || strcmp(AnglesMode, '9.1.6') + BRIR_data.Directions = [AzOrient; ElOrient]'; + end + + BRIR_data.AnglesMode = AnglesMode; +end + + +function Plot_data = init_Plot_data(SRIR_data, BRIR_data) + % Initialize visualization struct (from SDM Toolbox, with additions) + Plot_data = createVisualizationStruct(... + 'fs', SRIR_data.fs, 'DefaultRoom', 'Small', ... % or 'VerySmall, 'Medium', 'Large' + 'name', sprintf('%s_%s_%s', SRIR_data.Room, SRIR_data.SourcePos, SRIR_data.ReceiverPos)); + + Plot_data.PlotAnalysisFlag = false; % Flag to determine if analysis results should be plotted. + Plot_data.PlotExportFlag = false; % Flag to determine if analysis result plots should be exported. + Plot_data.DestinationPath = BRIR_data.DestinationPath; +end + + +function SRIR_data = analyze(SRIR_data, SDM_Struct, Plot_data) + %% Normalize RIRs + SRIR_data.P_RIR = SRIR_data.P_RIR ./ max(SRIR_data.P_RIR); + SRIR_data.Raw_RIR = SRIR_data.Raw_RIR ./ max(SRIR_data.Raw_RIR); + + %% Analysis + % Estimate directional information using SDM. This function is a wrapper of + % the SDM Toolbox DOA estimation (using TDOA analysis) to include some + % post-processing. The actual DOA estimation is performed by the SDMPar.m + % function of the SDM Toolbox. + SRIR_data = Analyze_SRIR(SRIR_data, SDM_Struct); + + if Plot_data.PlotAnalysisFlag + % Plot analysis results + plot_name = [Plot_data.name, '_spatio_temporal']; + if SRIR_data.AlignDOA + plot_name = [plot_name, '_aligned']; + end + + Plot_Spec(SRIR_data, Plot_data, [Plot_data.name, '_time_frequency']); + Plot_DOA(SRIR_data, Plot_data, plot_name); + end +end + + +function BRIR_data=synthesize(SRIR_data, BRIR_data, Plot_data) + %% Synthesis + % 1. Pre-processing operations (massage HRTF directions, resolve DOA NaNs). + + % Read HRTF dataset for re-synthesis + HRIR_data = Read_HRTF(BRIR_data); + + [SRIR_data, BRIR_data, HRIR_data, HRIR] = PreProcess_Synthesize_SDM_Binaural(SRIR_data, BRIR_data, HRIR_data); + + % ----------------------------------------------------------------------- + % 2. Rotate DOA information, if required + if BRIR_data.DOAAzOffset || BRIR_data.DOAElOffset + SRIR_data = Rotate_DOA(SRIR_data, BRIR_data.DOAAzOffset, BRIR_data.DOAElOffset); + + if Plot_data.PlotAnalysisFlag + Plot_DOA(SRIR_data, Plot_data, [Plot_data.name, '_spatio_temporal_rotated']); + end + end + + % ----------------------------------------------------------------------- + % 3. Quantize DOA information, if required + if BRIR_data.QuantizeDOAFlag + [SRIR_data, ~] = QuantizeDOA(SRIR_data, ... + BRIR_data.DOADirections, SRIR_data.DOAOnsetLength); + + if Plot_data.PlotAnalysisFlag + Plot_DOA(SRIR_data, Plot_data, [Plot_data.name, '_spatio_temporal_quantized']); + end + end + + % ----------------------------------------------------------------------- + % 4. Compute parameters for RTMod Compensation + + % Synthesize one direction to extract the reverb compensation - solving the + % SDM synthesis spectral whitening + BRIR_Pre = Synthesize_SDM_Binaural(SRIR_data, BRIR_data, HRIR, [0, 0], true); + + % Using the pressure RIR as a reference for the reverberation compensation + BRIR_data.ReferenceBRIR = [SRIR_data.P_RIR, SRIR_data.P_RIR]; + + % Get the desired T30 from the Pressure RIR and the actual T30 from one + % rendered BRIR + [BRIR_data.DesiredT30, BRIR_data.OriginalT30, BRIR_data.RTFreqVector] = GetReverbTime(SRIR_data, BRIR_Pre, BRIR_data.BandsPerOctave, BRIR_data.EqTxx); + + % ----------------------------------------------------------------------- + % 5. Render BRIRs with RTMod compensation for the specified directions + + nDirs = size(BRIR_data.Directions, 1); + + % Render early reflections + hbar = parfor_progressbar(nDirs, 'Please wait, rendering (step 1/2) ...', 'Name', Plot_data.name); + parfor iDir = 1 : nDirs + hbar.iterate(); %#ok + BRIR_early_temp = Synthesize_SDM_Binaural( ... + SRIR_data, BRIR_data, HRIR, BRIR_data.Directions(iDir, :), false); + BRIR_early(:, :, iDir) = Modify_Reverb_Slope(BRIR_data, BRIR_early_temp); + end + close(hbar); + + % Render late reverb + BRIR_late = Synthesize_SDM_Binaural(SRIR_data, BRIR_data, HRIR, [0, 0], true); + BRIR_late = Modify_Reverb_Slope(BRIR_data, BRIR_late, Plot_data); + + % Remove leading zeros + save('dbg_rem_dly.mat','BRIR_early', 'BRIR_late') + [BRIR_early, BRIR_late] = Remove_BRIR_Delay(BRIR_early, BRIR_late, -20); + + % ----------------------------------------------------------------------- + % 6. Apply AP processing for the late reverb + + % AllPass filtering for the late reverb (increasing diffuseness and + % smoothing out the EDC) + BRIR_data.allpass_delays = [37, 113, 215]; % in samples + BRIR_data.allpass_RT = [0.1, 0.1, 0.1]; % in seconds + + BRIR_late = Apply_Allpass(BRIR_late, BRIR_data); + + + % ----------------------------------------------------------------------- + % 7. Split the BRIRs into components by time windowing + [BRIR_DSER, BRIR_LR, BRIR_DS, BRIR_ER] = Split_BRIR(BRIR_early, BRIR_late, BRIR_data.MixingTime, BRIR_data.fs, 256); + if Plot_data.PlotAnalysisFlag + Plot_BRIR(BRIR_data, BRIR_DS, BRIR_ER, BRIR_LR, Plot_data); + end + + % ----------------------------------------------------------------------- + % 8. Save BRIRs + nFiles = BRIR_data.ExportSofaFlag + (BRIR_data.ExportWavFlag * nDirs) + 1; + hbar = parfor_progressbar(nFiles, 'Please wait, saving (step 2/2) ...', 'Name', Plot_data.name); + + if BRIR_data.ExportSofaFlag + hbar.iterate(); + Save_BRIR_sofa(BRIR_data, BRIR_DSER, BRIR_LR, sprintf('BRIR_%s_%s_%s_%s.sofa', SRIR_data.Room, SRIR_data.SourcePos, SRIR_data.ReceiverPos, BRIR_data.AnglesMode)); + end + + if BRIR_data.ExportWavFlag + for iDir = 1 : nDirs + hbar.iterate(); + if iDir == 1 % export identical late reverberation only once + BRIR_LR_export = BRIR_LR; + else + BRIR_LR_export = []; + end + Save_BRIR_wav(BRIR_data, BRIR_DS(:, :, iDir), BRIR_DSER(:, :, iDir), BRIR_ER(:, :, iDir), BRIR_LR_export, BRIR_data.Directions(iDir, :)); + end + end + hbar.iterate(); + SaveRenderingStructs(SRIR_data, BRIR_data); + close(hbar); +end + diff --git a/characterization/room_acoustics/binaural_sdm/README.md b/characterization/room_acoustics/binaural_sdm/README.md new file mode 100644 index 0000000000000000000000000000000000000000..908c12a984dcd385907a02d901d8382075228628 --- /dev/null +++ b/characterization/room_acoustics/binaural_sdm/README.md @@ -0,0 +1,89 @@ + + +# Binaural SDM for IVAS # + +This script is based on an example provided in the [BinauralSDM repository](https://github.com/facebookresearch/BinauralSDM). + +The BRIRs generated by this script include the steps described in *Amengual et al. 2020* — DOA Postprocessing (smoothing and quantization) and RTMod+AP equalization. + +--- + +### Overall Process + +1. **Spatial data extraction** + - Multichannel RIR data is processed using the SDM Toolbox (Tervo & Patynen). + +2. **DOA Postprocessing** + - Spatial information is smoothed and quantized as proposed in *Amengual et al. 2020*. + +3. **Preliminary BRIR synthesis** + - Closest directions are selected from an HRIR dataset (Orange HRIR dataset #53). + +4. **RTMod correction** + - Reverberation time is corrected using the RTMod method (*Amengual et al. 2020*). + +5. **All-pass filtering** + - A cascade of 3 all-pass filters is applied to both channels to improve diffuseness of the late reverberation tail. + +6. **Head orientation support** + - Responses for arbitrary head orientations are generated by rotating DOA vectors. + +7. **Data export** + - Results are saved in a user-defined folder. + +--- + +### How to use the script + +- Clone the [BinauralSDM repository](https://github.com/facebooklace). +- Copy the `BinauralSDM_for_IVAS.m` file into `./BinauralSDM/Src/Examples` folder. +- Update the settings at the top of the `BinauralSDM_for_IVAS` function to reflect the actual folder structure. +- Run the `BinauralSDM_for_IVAS` function to generate the BRIRs based on provided HRIRs and spatially recorded room impulse responses (RIRs). + +--- + +### References + +- Tervo et al. (2013) — *"Spatial Decomposition Method for Room Impulse Responses"*, JAES. +- SDM Toolbox — S. Tervo and J. Patynen, Mathworks Exchange. +- Amengual et al. (2020) — *"Optimizations of the Spatial Decomposition Method for Binaural Reproduction"*, JAES. + +Original Author: Sebastia Amengual, adapted by: Marek Szczerba (MS) + +Last modified: 08/25/2025 + +--- + +### Notes (MS) + +- Eigenmike geometry is supported out-of-the-box. Check `create_MicGeometry` in the `BinauralSDM/Src` folder. diff --git a/characterization/room_acoustics/reverb/calcPredelayRT60DSR.m b/characterization/room_acoustics/reverb/calcPredelayRT60DSR.m new file mode 100644 index 0000000000000000000000000000000000000000..ed133460d8ff0e7ae954f3813f763107532cf0b9 --- /dev/null +++ b/characterization/room_acoustics/reverb/calcPredelayRT60DSR.m @@ -0,0 +1,742 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [predelay, rt60, dsr, f0, refData] = calcPredelayRT60DSR(room, fs, predelayMethod, figureOffset, plotVerbosity) + % ##### Advanced configuration variables #### + + % Speed of sound + speedOfSound = 343; + + % Whether to do DC removal for predelay calculations. For the energy + % density profiles it is important for DC offsets to be removed. + % Specifically at the end of IRs where the values are small. This features + % attempts to do this. Typically only needed for measured RIRs. + dcRemovalForPredelay = 0; + + % Allows to optimize for speed when set to 0, only for debugging code after + % predelay calculation. + doPredelayCalcAlways = 1; + + % ##### Initialization #### + + % Make sure there is something to output + predelay = NaN; + rt60 = NaN; + dsr = NaN; + refData = NaN; + f0 = NaN; + + % Struct for controlling the plotting + if (nargin < 4) + figureOffset = 10000; + else + figureOffset = figureOffset * 10000; + end + + description = ''; + for i = 1:length(room.file_filters) + description = strcat(description, room.file_filters{i}); + if (i < length(room.file_filters)) + description = strcat(description, '|'); + end + end + + plotControl.figureOffset = figureOffset; + plotControl.description = description; + plotControl.plotVerbosity = plotVerbosity; + plotControl.maxSubplots = 8; + + % init echo density profile + edpConfig.T_edp = 60; % Max EDP length (sec) + edpConfig.winlen_edp = round(.05 * fs); % smoothing window + + + % init octave band filtering (center frequencies) + b = 1/3; + f0 = [20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250, 315, 400,... + 500, 630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000, 5000,... + 6300]; + if (fs >= 32e3) + f0 = [f0, 8000, 10000, 12500]; + end + if (fs == 48e3) + f0 = [f0, 16000, 20000]; + end + + % ########################## + % ### FIND AND LOAD RIRS ### + % ########################## + + % Process and filter file list + rir_dir = [room.irRoot, room.rname, '\']; + disp(rir_dir); + filteredFiles = filterFiles(rir_dir, room.file_filters); + + % Load the RIRs + rirs = loadRIRs(filteredFiles, room.sim_delay, fs); + + if (size(rirs, 2) == 0) + warndlg(['No RIRs found for [' description '].']); + return; + end + + if (size(rirs, 1) > 528128) + rirs = rirs(1:528128, :); + warndlg('RIRs longer than 11 seconds. Cropping to first 528128 samples (~11 sec).'); + end + + % ################################################ + % ### ENSURE A REFERENCE DISTANCE IS AVAILABLE ### + % ################################################ + + % If no reference distance is given, calculate it from the RIR. + if (isnan(room.rirRefDistance)) + [room.rirRefDistance, rirs, dirResponse] = calculateRefDistance(rirs, room.srcMicDistance, fs, speedOfSound, plotControl); + dirResponse(8192) = 0; + + plotControl.figureOffset = plotControl.figureOffset + size(rirs, 2); + else + dirResponse = [1; zeros(8192, 1)]; + end + + % ############################ + % ### PERFORM RIR ANALYSIS ### + % ############################ + + % Process RIR data to extract predelay, rt60, and ddr. + [predelay, rt60, dsr, refData] = getPredelayRT60DSR(room, edpConfig, filteredFiles, rirs, ... + b, f0, fs, speedOfSound, ... + room.rirRefDistance, dirResponse, ... + room.t60Offset, doPredelayCalcAlways, dcRemovalForPredelay, predelayMethod, plotControl); +end + +% ============================= Functions ================================= + +function [rirRefDistance, rirs, dirResponse] = calculateRefDistance(rirs, srcMicDistance, fs, speedOfSound, plotControl) + + % Window for measuring direct response energy + dirRespWindowLen = round(3e-3 * fs); % samples + srcMicDist = zeros(size(rirs, 2), 1); + minLength = length(rirs); + + for rirIdx = 1:size(rirs, 2) + % if (plotControl.plotVerbosity >= 9) + % figure(plotControl.figureOffset + rirIdx); clf; + % end + fprintf(1, 'File %d / %d\n', rirIdx, size(rirs, 2)); + + rir = rirs(:, rirIdx); + + % Find largest peak + [~, srcMicLag] = max(abs(rir)); + + % Find -20dB point before that peak. + find20dBPoint = true; + while (find20dBPoint) + t = srcMicLag; + t0 = t; + stop = false; + while (t0 > 1 && ~stop) + t0 = t0 - 1; + if (abs(rir(t)) / abs(rir(t0)) >= 10^(20/20)) + stop = true; + end + end + t1 = t0 + dirRespWindowLen; + + % If another peak before overall largest peak, use that and start over. + if (t0 > 1 && max(abs(rir(1:t0))) > abs(rir(srcMicLag)) * 10^(-3/20)) + [~, srcMicLag] = max(abs(rir(1:t0))); + else + find20dBPoint = false; + end + end + + if (isnan(srcMicDistance)) + % If not provided, calculate source to mic distance from the Time of + % Flight of the direct path response. + srcMicDist(rirIdx) = srcMicLag / fs * speedOfSound; + else + % If provided, use correct values. + %srcMicDist(rirIdx) = srcMicDistance(rirIdx); + srcMicDist(rirIdx) = srcMicDistance; + + % Align RIR to have direct response at expected lag. + srcMicLag_expected = round(srcMicDist(rirIdx) / speedOfSound * fs); + + if (srcMicLag_expected > srcMicLag) + rir = [zeros(srcMicLag_expected - srcMicLag, 1); rir(1:end - (srcMicLag_expected - srcMicLag))]; + rirs(:, rirIdx) = rir; + elseif (srcMicLag_expected < srcMicLag) + rir = [rir(srcMicLag - srcMicLag_expected + 1:end); zeros(srcMicLag - srcMicLag_expected, 1)]; + rirs(:, rirIdx) = rir; + end + + t0 = t0 + (srcMicLag_expected - srcMicLag); + t1 = t1 + (srcMicLag_expected - srcMicLag); + + minLength = min(minLength, length(rirs) + srcMicLag_expected - srcMicLag); + + srcMicLag = srcMicLag_expected; + end + + % Measure the direct response energy in a range around the lag with + % maximum amplitude. A 3ms window starting at the -20 dB point before the + % first peak. + dirRespWindow = [t0:t1]; + dirResponse(:, rirIdx) = rir(dirRespWindow); + E_dirResp(rirIdx) = sum(dirResponse(:, rirIdx).^2); + + fprintf(1, 'Source to mic distance: %.2f m\n', srcMicDist(rirIdx)); + fprintf(1, 'Direct response energy: %.2g\n', E_dirResp(rirIdx)); + end + + % Reference distance is the source-mic distance where the direct path + % response has unit energy. + % I.e. the distance where the PCM signal filtered by the RIR essentially is + % not modified. + % + % Normalize all the rirs equally such that their loudest direct response energy is 1. + + % Truncate zero-padded samples at the end + rirs(minLength + 1:end, :) = []; + + [E_maxDirResp, maxDirRespIdx] = max(E_dirResp); + rirs = rirs .* sqrt(1./ E_maxDirResp); + + dirResponse = dirResponse(:, maxDirRespIdx) .* sqrt(1./ E_maxDirResp); + + + % That makes the corresponding source-mic distance the reference distance + rirRefDistance = srcMicDist(maxDirRespIdx); + + fprintf(1, 'Reference distance: %.2f m\n', rirRefDistance); +end + +% ========================================================================= + +function [predelay_s, rt60, ddr, refData] = getPredelayRT60DSR(room, edpConfig, files, rirs, ... + b, f0, fs, c, ... + rirRefDistance, dirResponse, ... + t60Offset, forcePredelayCalc, dcRemovalForPredelay, predelayMethod, plotControl) + % Main function to get predelay, rt60 and ddr + + maxDatenum = 0; + for fIdx = 1:length(files) + d = dir(files{fIdx}); + if (d.datenum > maxDatenum) + maxDatenum = d.datenum; + end + end + + iniFile = mfilename('fullpath'); + iniFile = [iniFile(1:end-4) '.mat']; + if (exist(iniFile, 'file')) + iniData = load(iniFile); + else + iniData.storedDatenum = 0; + end + + % get pre-delay + if (iniData.storedDatenum == maxDatenum && ~forcePredelayCalc) + warning('Using the previously calculated predelay') + predelay_s = iniData.predelay_s; + else + % Determine predelay from analysis + [predelay_s, doPredelayPlot] = get_predelay(rirs, room.room_dim, edpConfig, c, fs, dcRemovalForPredelay, predelayMethod, plotControl); + + switch (room.predelayScheme) + case 'ShiftFactor4' + predelay_s_meas = predelay_s / 4; % P13 style. Shifts reverb respose to 1/4 of the EIF predelay + case 'DSRScaling' + predelay_s_meas = predelay_s; % P27-style. No change needed. + case 'Fixed_50ms' + predelay_s_meas = 0.05; + case 'Fixed_16ms' + predelay_s_meas = 0.016; + otherwise + error('Selected predelay scheme not supported'); + end + + if (doPredelayPlot && plotControl.plotVerbosity >= 9) + % Plot RIR + fig = figure(plotControl.figureOffset + 1); + clf; + plot((0:length(rirs(:, 1))-1)./fs, rirs(:, 1)); + hold on; + plot(predelay_s * ones(20,1), linspace(-.5, 1, 20), '--', 'Color', fig.Children.XColor); + grid on; + xlim([0 1]); + title(plotControl.description, 'Interpreter', 'none'); + end + end + + fprintf(1, 'Predelay: %.0f ms\n', predelay_s * 1000); + + iniData.predelay_s = predelay_s; + iniData.storedDatenum = maxDatenum; + + save(iniFile, '-struct', 'iniData'); + + % Get T60 and DSR + [rt60, ddr, refData] = getRT60DSR(rirs, predelay_s_meas, t60Offset, rirRefDistance, dirResponse, b, f0, fs, room, plotControl); +end + +% ========================================================================= + +function [predelay_s, doPlot] = get_predelay(rirs, room_dim, edpConfig, c, fs, dcRemovalForPredelay, predelayMethod, plotControl) + + % get predelay from energy decay curve + + EDP = []; + T_edp = edpConfig.T_edp; + winlen_edp = edpConfig.winlen_edp; + minlen = round(T_edp*fs); + L = round(60*fs); + + % extract room data + max_diag = sqrt(sum(room_dim.^2)); + min_delay = max_diag/c; + + minPredelay = Inf; + + doPlot = true; + switch predelayMethod + case 'AbelHuangEDP' + % echo density profile + fig = figure(1001); + clf; + hold on; + for jj = 1:size(rirs, 2) + fprintf(1, 'File %d / %d\n', jj, size(rirs, 2)); + + rir = rirs(:, jj); + + [~, dirRespLag] = max(abs(rir)); + minPredelay = min(dirRespLag, minPredelay); + + % 2. get and collect the echo density profile + edpConfig = get_echo_density_profile(rir, winlen_edp, dcRemovalForPredelay)'; + edpConfig(edpConfig>1) = 1; % suppress spurious spikes + edpConfig = [edpConfig,zeros(1,L-length(edpConfig))]; + EDP = [EDP;edpConfig]; + minlen = min(minlen, length(edpConfig)); + + figure(1001); + plot(edpConfig(1:min(end, 2 * fs)), 'k') + end + + % 3. Get predelay from ensemble average echo density profile + mean_edp = mean(EDP(:,1:minlen),1); + + figure(1001); + plot(mean_edp(1:min(end, 2 * fs)), 'r', 'LineWidth', 2); + + predelayLag = find(mean_edp(minPredelay:end) >= 1, 1, 'first') + (minPredelay - 1); + + fprintf(1, 'First instance mean EDP increases past 1: %d\n', predelayLag) + + if (~isempty(predelayLag)) + predelay_s = predelayLag ./ fs; + + % Plot EDP and RIR + fig = figure(1002); + clf; + plot((0:length(mean_edp)-1)./fs, mean_edp); hold on; % xlim([0 xdata(end)]); + hold on; + plot((0:length(rir)-1)./fs, rir); + plot(predelay_s * ones(20,1), linspace(-.5, 1, 20), '--', 'Color', fig.Children.XColor); + title(plotControl.description, 'Echo density profile'); + grid on; + xlim([0 1]) + else + mean_edp = mean_edp - mean_edp(1); + + xdata = (0:minlen-1)./fs; + ydata = mean_edp; + + [~, idx] = max(ydata); + % thresh = xdata(idx + 0.05*fs); % limit the window over which to line-fit + thresh = xdata(idx + 0.15*fs); % Limit the window over which to line-fit - JK: increasing end point for linear fit + ydata = ydata(xdata < thresh); + xdata = xdata(xdata < thresh); + + + % Plot EDP and RIR + fig = figure(1002); + clf; + plot((0:length(mean_edp)-1)./fs, mean_edp); + % xlim([0 xdata(end)]); + hold on; + plot((0:length(rir)-1)./fs, rir); + + + % Piece-wise Linear fit on the data. First segment must be first order + % with offset 0 and second segment must be zero order. + F = @(B,xdata) min(B(1)*xdata,B(2)); % Format of the equation + IC = [0.6 0.8]; % Initial guess + B = lsqcurvefit(F, IC, xdata, ydata, [1 0.1],[50 5.0]); + plot(xdata, F(B, xdata)); + + + % Calculate the predelay based on piecewise linear model and + % the maximum room diagonal delay. Assume predelay >= min_delay + predelay_s = max(B(2)/B(1), min_delay); + + % convert predelay to samples and add to EDP plot + pd = round(fs*min(predelay_s, xdata(end))); + plot(pd/fs*ones(20,1), linspace(-.5, 1, 20), '--', 'Color', fig.Children.XColor); + title(plotControl.description, 'Echo density profile');grid on; xlim([0 max(xdata)]) + end + doPlot = false; + case 'LindauEquation' + V = prod(room_dim); + predelay_s = (0.58 * sqrt(V) + 21.2) * 1e-3; + + case 'AndreasEquation' + V = prod(room_dim); + predelay_s = (0.89 * sqrt(V)) * 1e-3; + + case 'ToFEquation' + D = max(room_dim); + predelay_s = 4 * D / c; + + otherwise + error('Unsupported method'); + end +end + +% ========================================================================= + + +function [rt60_avg, DSR_avg, refData] = getRT60DSR(rirs, predelay_sec, t60Offset, rirRefDistance, dirResponse, b, f0, fs, room, plotControl) + % requires predelay value to reprocess the RIRs + + if (isfield(room, 'referenceEIF') && exist(room.referenceEIF, 'file')) + refData = getReferenceData(room, fs); + else + refData = []; + end + + if (isfield(room, 'revDistanceGain')) + revDistanceGain = room.revDistanceGain; + else + revDistanceGain = 1; + end + + os = round(fs * t60Offset); % t60 measurement offset in samples + pd = round(fs * predelay_sec); % predelay in samples + nrRIRs = size(rirs, 2); + + % 4. calculate ddr using predelay (averaged across the calculated rir) + E_diffuse = zeros(length(f0), nrRIRs); + E_bnd = zeros(length(f0), nrRIRs); + rt60_os = zeros(length(f0), nrRIRs); + rt60 = zeros(length(f0), nrRIRs); + + [~, bnd_1kHz] = min(abs(f0 - 1000)); + + if (plotControl.plotVerbosity >= 9) + figure(plotControl.figureOffset + 2); + clf; + end + + nSubplots = min(nrRIRs, plotControl.maxSubplots); + nrY = 1; + nrY(nSubplots > 3) = 2; + nrX = ceil(nSubplots / nrY); + + nTruncate = fs * room.tTruncate; + + % Compute RT60 per frequency band + for rirIdx = 1:nrRIRs + fprintf(1, 'File %d / %d\n', rirIdx, nrRIRs); + + rir = rirs(:, rirIdx); + rir_os = rir(os:end); + + T60_os = zeros(length(f0), 1); + T60 = zeros(length(f0), 1); + + if (plotControl.plotVerbosity >= 9 && rirIdx <= nSubplots) + figure(plotControl.figureOffset + 2); + axesHandle = subplot(nrY, nrX, rirIdx); + vertLineObj(1) = plot(axesHandle, [t60Offset t60Offset], [-100 0], 'r--'); + hold on; + vertLineObj(2) = plot(axesHandle, [predelay_sec predelay_sec], [-100 0], '--', 'Color', axesHandle.XColor); + end + validBndList = []; + for bnd=1:length(f0) + [hf, ~] = octaveBandFilterResponse(rir, fs, f0(bnd), b); + hf_os = hf(os:end); + hf_pd = hf(pd:end); + + if (any(isnan(hf_os)) || sum(hf_os.^2) > sum(rir_os.^2)) + warning('Something went wrong with the octavebandfilterresponse()'); + T60_os(bnd) = 0; + else + validBndList = [validBndList; bnd]; + [T60_os(bnd), EDC_os] = processGetRT60(hf_os, fs, [0 -30], nTruncate); + end + + if (any(isnan(hf)) || sum(hf.^2) > sum(rir.^2)) + warning('Something went wrong with the octavebandfilterresponse()'); + T60(bnd) = 0; + else + [T60(bnd), EDC] = processGetRT60(hf, fs, [-10 -40], nTruncate); + end + + if (plotControl.plotVerbosity >= 9 && rirIdx <= nSubplots) + % plot(axesHandle, [0:length(EDC)-1] ./ fs, EDC - EDC(1)) + plot(axesHandle, [0:length(EDC)-1] ./ fs, EDC) + hold on; + measEnd = min(round(os + (T60_os(bnd) / 2) * fs), nTruncate); + % plot(axesHandle, measEnd ./ fs, EDC(min(measEnd, end)) - EDC(1), 'rx') + plot(axesHandle, measEnd ./ fs, EDC(min(measEnd, end)), 'rx'); + + grid on; + if (rirIdx == 1) + % title(plotControl.description, 'Normalized EDC for each band'); + title(plotControl.description, 'EDC for each band'); + end + xlabel('Time [s]'); + ylabel('EDC [dB]'); + ylim([-100 0]); + end + + rt60_os(bnd, rirIdx) = rt60_os(bnd, rirIdx) + T60_os(bnd); + rt60(bnd, rirIdx) = rt60(bnd, rirIdx) + T60(bnd); + + E_bnd(bnd, rirIdx) = sum(hf_os.^2); + E_diffuse(bnd, rirIdx) = sum(abs(hf_pd).^2) ./ revDistanceGain.^2; % Compensate for distance gain applied in measurement + + if (bnd == bnd_1kHz && plotControl.plotVerbosity >= 9) + figure(plotControl.figureOffset + 3); + % clf; + plot((0:length(EDC)-1) ./ fs, EDC, 'LineWidth', 1) + hold on + plot(((0:length(EDC_os)-1) + os) ./ fs, EDC_os, 'LineWidth', 1) + plot(t60Offset * [1 1], [-100 0], 'r--'); + ylim([-100 0]) + legend(['EDC(1kHz) - Whole IR - T60=' num2str(T60(bnd))], ['EDC(1kHz) - after offset - T60=' num2str(T60_os(bnd))]); + grid on; + title(plotControl.description, 'EDC at 1 kHz'); + figure(plotControl.figureOffset + 2); + end + end + fprintf(1, 'Captured %.1f %% of the IR''s energy in the bands.\n', sum(E_bnd(validBndList, rirIdx)) ./ sum(rir_os.^2) * 100); + % if (plotControl.plotVerbosity >= 9) + % legend(vertLineObj, 'T60 measurement offset', 'DSR measurement offset (=predelay)'); + % end + end + + % 5. Get the total source energy for DDR calculation + sourceList{1}.dBPreGain = 0; + sourceList{1}.refDistance = rirRefDistance; % Distance for which direct path response's energy is 1. + sourceList{1}.directivityElem = ''; % = Omni-directional + + normSourceEnergy = calculateNormalizedSourceEnergy(sourceList{1}); + + + % 6. Normalize power spectra using total source energy (incl. octave + % filter response) get DSR according to the agreed-upon definition + if (isfield(room, 'measuredSource')) + E_source = interp1(room.measuredSource.frequencies, room.measuredSource.avgDirecEnergyValues, f0(:)); + E_source = repmat(E_source, 1, nrRIRs); + else + E_source = ones(length(f0),nrRIRs); + end + + DSR = zeros(length(f0),nrRIRs); + for bnd = 1:length(f0) + [hf_src, ~] = octaveBandFilterResponse(dirResponse, fs, f0(bnd), b); + E_source(bnd, 1) = E_source(bnd, 1) * normSourceEnergy * sum(abs(hf_src.^2)); + for rirIdx = 1:nrRIRs + DSR(bnd, rirIdx) = DSR(bnd, rirIdx) + E_diffuse(bnd, rirIdx) / E_source(bnd); + + if (rirIdx == 1 && f0(bnd) == 1000) + fprintf(1, 'At 1 kHz. E_diffuse = %.2g -- E_source = %.2g\n', E_diffuse(bnd), E_source(bnd)); + end + end + end + + % Compute average Schröder frequency for all channels + freqSchroeder = 2000 * sqrt(mean(rt60(bnd_1kHz,:))./prod(room.room_dim)); + freqSchroeder_idx = find(f0 >= 4.0 * freqSchroeder, 1); + + for rirIdx = 1:nrRIRs + rt60_lo = rt60(freqSchroeder_idx, rirIdx); + rt60_os_lo = rt60_os(freqSchroeder_idx, rirIdx); + DSR_lo = DSR(bnd, rirIdx); + for bnd = 1:length(find(f0 < 4.0 * freqSchroeder)) + rt60(bnd, rirIdx) = rt60_lo; + rt60_os(bnd, rirIdx) = rt60_os_lo; + DSR(bnd, rirIdx) = DSR_lo; + end + end + + % Take the average of all files + rt60_os_avg = mean(standardizeMissing(rt60_os, [Inf, -Inf]), 2, "omitnan"); + rt60_avg = mean(standardizeMissing(rt60, [Inf, -Inf]), 2, "omitnan"); + + if (plotControl.plotVerbosity >= 3) + figure(plotControl.figureOffset + 4); clf; + % set(gcf, 'Position', [1336 732 990 609]); + subplot(211); + semilogx(f0, rt60_os_avg) + hold on; + plot(f0, rt60_avg) + % plot(refData.frequencies, refData.t60); + title(plotControl.description, 'Average RT60'); + legend('RT60 - 0-30 dB decay extrapolation after offset', ...'RT60 - 10-40 dB decay extrapolation on whole IR', + 'EIF reference T60', 'Location', 'Best'); + xlabel('Frequency [Hz]') + ylabel('RT60 [s]') + xlim([20 24000]) + grid on; + + subplot(212); + semilogx(f0, rt60_os) + title('', 'RT60 of all RIRs') + xlabel('Frequency [Hz]') + ylabel('RT60 [s]') + xlim([20 24000]) + grid on; + + if (nrRIRs > 1 && plotControl.plotVerbosity >= 9) + figure(plotControl.figureOffset + 5); clf; + % semilogx(mean(rt60_os, 1)) + semilogx(rt60_os_avg); + title(plotControl.description, 'Broadband RT60 per RIR'); + xlabel('RIR index'); + ylabel('RT60 [s]'); + grid on + end + end + + if (plotControl.plotVerbosity >= 2) + figure(plotControl.figureOffset + 6); clf; + % set(gcf, 'Position', [1335 48 990 609]); + subplot(221); + semilogx(f0, 10*log10(E_diffuse(:, 1))) + hold on + semilogx(f0, 10*log10(E_source(:, 1))) + semilogx(f0, 10*log10(DSR(:, 1))) + legend('Diffuse energy', 'Source energy', 'DSR'); + xlabel('Frequency [Hz]') + xlim([20 24000]) + grid on; + title(plotControl.description); + + subplot(222); + semilogx(f0, 10*log10(DSR)) + title('', 'DSR of all RIRs'); + xlabel('Frequency [Hz]') + xlim([20 24000]) + ylabel('DSR [dB]'); + grid on; + + subplot(223); + semilogx(f0, 10*log10(mean(DSR, 2))); + hold on; + plot(f0, 10*log10(mean(DSR, 2))); + % plot(refData.frequencies, 10*log10(refData.ddr)) + title('', 'Mean DSR'); + legend('Mean DSR', 'EIF reference DSR', 'Location', 'Best'); + xlabel('Frequency [Hz]') + xlim([20 24000]) + ylabel('DSR [dB]'); + grid on; + + subplot(224); + refDSRInterp = interp1(f0, DSR, f0); + % refDSRInterp = interp1(refData.frequencies, refData.ddr, f0); + % semilogx(f0, 10*log10(mean(DSR, 2)) - 10*log10(refDSRInterp(:))) + semilogx(f0, 10*log10(mean(DSR, 2)) - 10*log10(mean(refDSRInterp, 2))) + title('', 'DSR difference'); + legend('RIR DSR - refDSR', 'Location', 'Best'); + xlabel('Frequency [Hz]') + xlim([20 24000]) + ylabel('DSR [dB]'); + grid on; + end + + if (plotControl.plotVerbosity >= 1) + figure(plotControl.figureOffset + 7); clf; + % set(gcf, 'Position', [1336 732 990 609]); + subplot(211); + semilogx(f0, rt60_os_avg) + hold on; + plot(f0, rt60_avg) + % plot(refData.frequencies, refData.t60); + title(plotControl.description, 'Average RT60'); + legend('RT60 - 0-30 dB decay extrapolation after offset', ...'RT60 - 10-40 dB decay extrapolation on whole IR', + 'EIF reference T60', 'Location', 'Best'); + xlabel('Frequency [Hz]') + ylabel('RT60 [s]') + xlim([20 24000]) + grid on; + + subplot(212); + semilogx(f0, 10*log10(mean(DSR, 2))); + hold on; + % plot(f0, 10*log10(mean(DSR))); + plot(f0, 10*log10(mean(DSR, 2))); + % plot(refData.frequencies, 10*log10(refData.ddr)) + title('', 'Mean DSR'); + legend('Mean DSR', 'EIF reference DSR', 'Location', 'Best'); + xlabel('Frequency [Hz]') + xlim([20 24000]) + ylabel('DSR [dB]'); + grid on; + end + + % Average over all RIRs + E_diffuse = mean(E_diffuse, 2); + DSR_avg = mean(DSR, 2); +end + +% ========================================================================= + +function refData = getReferenceData(room, fs) + eifData = parseEIF(room.referenceEIF, fs); + refData = eifData.acousticEnvironment{room.referenceRoomIdx}.acousticParameters{1}; +end + +% ========================================================================= + + diff --git a/characterization/room_acoustics/reverb/calculateNormalizedSourceEnergy.m b/characterization/room_acoustics/reverb/calculateNormalizedSourceEnergy.m new file mode 100644 index 0000000000000000000000000000000000000000..d09684b40950c6decac0101f15956202a98b418a --- /dev/null +++ b/characterization/room_acoustics/reverb/calculateNormalizedSourceEnergy.m @@ -0,0 +1,51 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function srcNormEnergyWideband = calculateNormalizedSourceEnergy(source) + +% Surface integral of sphere around omnidirectional source at refDistance +% results in the surface area of a sphere with r = refDistance. Normalize +% to surface area associated to PCM signal (human ear), to get scaling +% factor. +dirEnergyScale = 4 * pi * source.refDistance^2 ./ (0.05 * 0.04 / 2); + + +% Include effect of pregain. +srcNormEnergyWideband = dirEnergyScale .* 10^(source.dBPreGain / 10); + + + diff --git a/characterization/room_acoustics/reverb/filterFiles.m b/characterization/room_acoustics/reverb/filterFiles.m new file mode 100644 index 0000000000000000000000000000000000000000..616a1f5d04f67ac27e3f04bf476bcfa3d59e7564 --- /dev/null +++ b/characterization/room_acoustics/reverb/filterFiles.m @@ -0,0 +1,64 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +function filtered = filterFiles(file_dir, file_filters) +% filter out files from git repo based on specific keywords + +files = dir(file_dir); +% filter out only desired files.. +filtered = {}; +nfiles = 0; +for jj=1:length(files) + + fname = files(jj).name; + res = 1; + for kk=1:length(file_filters) + if isempty(strfind(fname, file_filters{kk})) + res = 0; + end + end + + if res == 1 + filtered{nfiles+1} = strcat(file_dir, fname); + nfiles = nfiles+1; + end +end + +% ========================================================================= + + diff --git a/characterization/room_acoustics/reverb/loadRIRs.m b/characterization/room_acoustics/reverb/loadRIRs.m new file mode 100644 index 0000000000000000000000000000000000000000..674dfc05816bf1296978bd9cbf40057ea2f4a690 --- /dev/null +++ b/characterization/room_acoustics/reverb/loadRIRs.m @@ -0,0 +1,73 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function rirSet = loadRIRs(files, sim_delay, fs) + +nCh = 0; +maxLength = 0; +for fileIdx = 1:length(files) + [~, fname, ext] = fileparts(files{fileIdx}); + if (strcmpi(ext, '.wav')) + info = audioinfo(files{fileIdx}); + maxLength = max(maxLength, info.TotalSamples); + nCh = nCh + info.NumChannels; + else + error('Unsupported file format: %s', ext); + end +end + +rirSet = zeros(maxLength, nCh); +chIdx = 1; % Channel index initialization +for fileIdx = 1:length(files) + [~, fname, ext] = fileparts(files{fileIdx}); + + fprintf(1, 'RIR index: %2d -- File %2d / %2d -- %s\n', fileIdx, fileIdx, length(files), [fname ext]); + [rir, fs_file] = audioread(files{fileIdx}); + + if (fs ~= fs_file) + error('Samplerate in wav file differs from configured sample rate.'); + end + + % Remove known delay + rir = rir(sim_delay + 1:end, :); + + % Extract the channels + for n = 1:size(rir, 2) + rirSet(1:length(rir), chIdx) = rir(:, n); + chIdx = chIdx + 1; + end +end \ No newline at end of file diff --git a/characterization/room_acoustics/reverb/octaveBandFilter.m b/characterization/room_acoustics/reverb/octaveBandFilter.m new file mode 100644 index 0000000000000000000000000000000000000000..0969302fe1971ab3e9bd8066e8d299c2272eab3d --- /dev/null +++ b/characterization/room_acoustics/reverb/octaveBandFilter.m @@ -0,0 +1,77 @@ + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Contact: Jeroen Koppens (jeroen.koppens@philips.com) +% Marek Szczerba (marek.szczerba@philips.com) +% +% function [B, A] = octaveBandFilter(fc, width, fs, [N]) +% +% Returns the coefficients of a bandfilter, based on a 2Nth Butterworth +% filter [1]. +% +% fc : Center frequency (Hz) +% width : Filter width (octaves), for example 1/3 = one-third-octave +% fs : Sampling frequency (Hz) +% N : Filter order (default: 3) +% +% REFERENCES: +% +% [1] Couvreur, C. (1997?), "Implementation of a one-third-octave filter +% bank in MATLAB", Applied Acoustics + +function [B, A] = octaveBandFilter(fc, width, fs, N) + +% Check arguments + +if nargin < 3 + error('Too few arguments') +end + +% Default order? +% +if nargin < 4 + N = 3; + end + +half_width = width/2; + +f1 = fc / (2^half_width); +f2 = fc * (2^half_width); + +Wn = 2*[f1 f2]/fs; +[B,A] = butter(N,Wn); diff --git a/characterization/room_acoustics/reverb/octaveBandFilterDelay.m b/characterization/room_acoustics/reverb/octaveBandFilterDelay.m new file mode 100644 index 0000000000000000000000000000000000000000..d1605a2153b2bcbeadac53f5e87545a661efc3f3 --- /dev/null +++ b/characterization/room_acoustics/reverb/octaveBandFilterDelay.m @@ -0,0 +1,57 @@ + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Contact: Jeroen Koppens (jeroen.koppens@philips.com) +% Marek Szczerba (marek.szczerba@phhilips.com) +% +% function N_delay = octaveBandFilterDelay(B,A) +% +% Date: Aug 2010 +% Author: Jasper van Dorp Schuitman +% +% Estimates the delay of filter with coefficients B and A + +function N_delay = octaveBandFilterDelay(B,A) + +N = 1024; +x = zeros(N,1); +x(1) = 1; +y = filter(B,A,x); + +[mx,M] = max(abs(y)); +N_delay = M-1; diff --git a/characterization/room_acoustics/reverb/octaveBandFilterResponse.m b/characterization/room_acoustics/reverb/octaveBandFilterResponse.m new file mode 100644 index 0000000000000000000000000000000000000000..7baa4456a8cef140bbf0bee09ca1d63bfd331343 --- /dev/null +++ b/characterization/room_acoustics/reverb/octaveBandFilterResponse.m @@ -0,0 +1,85 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Contact: Jeroen Koppens (jeroen.koppens@philips.com) +% Marek Szczerba (marek.szczerba@philips.com) +% +% function [hf, N_delay] = octaveBandFilterResponse(h, fs, fc, [b]) +% +% Date: Aug 2010 +% Author: Jasper van Dorp Schuitman +% +% Filters impulse response 'h' +% +% fs: Sampling frequency (Hz) +% fc: Center frequency (Hz) +% b: Bandwidth (octaves, default: 1) +% N_delay: Delay of filter in samples + +function [hf, N_delay] = octaveBandFilterResponse(h, fs, fc, b) + +%% Check input arguments + +if nargin < 3 + error('Too few input arguments, see help for details'); +end + +% Default bandwidth +if nargin < 4 + b = 1; +else + if isempty(b) + b = 1; + end +end + +%% Construct filter +if (fc < 100) + [B,A] = octaveBandFilter(fc, b, fs, 2); +else + [B,A] = octaveBandFilter(fc, b, fs); +end + +%% Pad h with extra zeros to account for the filter delay + +Nc = size(h,2); +N_delay = octaveBandFilterDelay(B,A); +h = [h; zeros(N_delay, Nc)]; + +%% Filter response + +hf = filter(B, A, h); diff --git a/characterization/room_acoustics/reverb/processGetRT60.m b/characterization/room_acoustics/reverb/processGetRT60.m new file mode 100644 index 0000000000000000000000000000000000000000..bd02cc5591951644e1ceaf3c7e8fe11c887168f0 --- /dev/null +++ b/characterization/room_acoustics/reverb/processGetRT60.m @@ -0,0 +1,64 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Contact: Jeroen Koppens (jeroen.koppens@philips.com) +% Marek Szczerba (marek.szczerba@philips.com) +% + +function [rt60, EDC] = processGetRT60(ir, fs, assessAtLevels, nTruncate) + +lvl_hi = assessAtLevels(1); +lvl_lo = assessAtLevels(2); + +EDC = 10*log10(flipud(cumsum(flipud(ir.^2)))); + +if nargin == 4 && ~isempty(nTruncate) + EDC = EDC(1:nTruncate); +end + +EDC_norm = EDC - EDC(1); + +[~, idx_hiLvl] = min(abs(EDC_norm - lvl_hi)); +[~, idx_loLvl] = min(abs(EDC_norm - lvl_lo)); + +% In case it never reaches high or low +if EDC_norm(end) > min(lvl_lo, lvl_hi) + p = [-0, EDC_norm(end)]; +else + p = polyfit((idx_hiLvl:idx_loLvl).', EDC_norm(idx_hiLvl:idx_loLvl), 1); +end +rt60 = -60 ./ p(1) / fs; diff --git a/characterization/room_acoustics/reverb/roomAcousticsCalc.m b/characterization/room_acoustics/reverb/roomAcousticsCalc.m new file mode 100644 index 0000000000000000000000000000000000000000..bf6656b2d7e28cf2d9164409da4f5a24b79ccfd1 --- /dev/null +++ b/characterization/room_acoustics/reverb/roomAcousticsCalc.m @@ -0,0 +1,287 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% (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. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% +% Contact: Jeroen Koppens (jeroen.koppens@philips.com) +% Marek Szczerba (marek.szczerba@philips.com) +% +% USAGE +% ===== +% +% Provide input by editing the variables in this m file. Provide the +% following data: +% +% Reference to one or more RIR files. +% - irRoot +% - room.rname Subfolder in irRoot with RIR files. +% - file_filters Cell array of strings, used to filter the files in +% the given room{end + 1}.rname subfolder. E.g. +% {'_48000'} +% {'RIR', '.wav'} +% {'524288rays', '48kHz', '.csv'} +% Files should be provided as CSV files or wav files. +% +% +% Required metadata in addition to RIR +% - room.room_dim Room dimensions. Needed for some predelay methods. +% - room.aparams String to be used in EIF element printout for +% 'position' attribute of AcousticParameters. +% - sim_delay Simulator or recording delay in samples. Used to +% compensate for the delay imposed by the simulator. +% If srcMicDistance is also given, actually compensated +% delay may differ. +% - rirRefDistance Indicates how RIR levels are normalized. I.e. for what +% source-mic distance the direct response is 1. +% If unknown, set to NaN and provide srcMicDistance +% when available. +% - srcMicDistance Indicates the distance in meters between the source +% and microphone during the measurement. Provide vector +% when multiple measurements available. +% If unknown, set to NaN and ensure that leading lag in +% RIR is representative of the time-of-flight delay for +% the direct part. E.g. by synchronously playing test +% sequence and recording microphone output. +% (Contact author for a Matlab mex function that plays +% test sequence and records response simultaneously) +% Parameter is not needed (i.e. set to NaN) when +% rirRefDistance is provided. +% +% Configuration +% - fs Sample rate of the RIRs, and for which to +% generate banded metadata. +% - predelayMethod Method to determine predelay. Possible values: +% {'LindauEquation', 'AndreasEquation', 'AbelHuangEDP', 'ToFEquation'} +% LindauEquation : (0.58 * sqrt(V) + 21.2) / 1000 +% AndreasEquation : (0.89 * sqrt(V)) / 1000 +% ToFEquation : 4 * d(longest dim) / 343 +% AbelHuangEDP : Based on echo density +% profile. First lag EDP >= 1. If no value +% above 1: knee of piecewise linear fit. +% - b, f0 Bandwidth in terms of octave bands and center +% frequencies for which to derive RT60 and DDR. +% +% +% Configuration suggested to keep at default values. +% - dcRemovalForPredelay Whether to do DC removal or not. Default = 0 +% - doPredelayCalcAlways Whether to reuse previous predelay calculation. Default = 1. +% - speedOfSound The speed of sound in m/s. Default = 343. +% + + +function roomAcousticsCalc() + + % ##### Configuration variables #### + + % Controls how DSR is measured, depending on how the plugin handled changing predelay. + % Options: {'ShiftFactor4', 'DSRScaling', 'Fixed_50ms'} + % ShiftFactor4 - For plugins that divide predelay by four (and don't scale DSR) + % DSRScaling - For plugins that scale DSR so that the reverb level + % after predelay remains stable, regardless of predelay change. + % Fixed_50ms - Starts measurement for reverb energy at 50 ms. + predelayScheme = 'DSRScaling'; + + % Delay of the simulation that should be removed (e.g. in the input file + % and/or HRTF). + sim_delay = 0; % In samples + + % Reference distance in meters. Use NaN to estimate from IR. + rirRefDistance = NaN; + + % Source to microphone distance in meters + srcMicDistance = 1.5; + + % Filters for selecting certain versions of the IRs + + % Root folder for the IRs + irRoot = '.\IR_recordings\Qualcomm'; + fileFilters = {'ir_48kHz.wav'}; + predelayScheme = 'DSRScaling'; + rendCfgPath = '.\conf'; + + % Sampling rate + fs = 48e3; + + % Truncation [s] + tTruncate = 1; + + % Which method to use for determining predelay + % [Default: time-of-flight] + predelayMethod = 'LindauEquation'; % {'LindauEquation', 'AndreasEquation', 'ToFEquation', 'AbelHuangEDP'} + + % How much of the plots should be made. + % 0 - No plots + % 1 - Only summary plot of T60 and DSR + % 2 - Summary plot, DSR overview plot + % 3 - Summary plot, DSR overview plot, T60 overview plot + % 9 - All plots + plotVerbosity = 9; + + % Init room(s) (select subsets by uncommenting/commenting) + room = cell(0,1); + + % =============== Qualcomm: Auditorium =============== + room{end + 1}.rname = ''; + room{end}.IDs = {'auditorium_L'}; + room{end}.room_dim = [19.7 18.625 5.175]; + room{end}.aparams = '0 1.5 0'; + room{end}.referenceRoomIdx = 1; + room{end}.t60Offset = 0.04; % In seconds + + % =============== Qualcomm: Kitchen with curtains =============== + room{end + 1}.rname = ''; + room{end}.IDs = {'kitchen_L_curtains'}; + room{end}.room_dim = [6.5 4 3]; + room{end}.aparams = '0 1.5 0'; + room{end}.referenceRoomIdx = 1; + room{end}.t60Offset = 0.04; % In seconds + + % =============== Qualcomm: Kitchen without curtains =============== + room{end + 1}.rname = ''; + room{end}.IDs = {'kitchen_L_nocurtains'}; + room{end}.room_dim = [6.5 4 3]; + room{end}.aparams = '0 1.5 0'; + room{end}.referenceRoomIdx = 1; + room{end}.t60Offset = 0.04; % In seconds + + % =============== Qualcomm: Listening lab =============== + room{end + 1}.rname = ''; + room{end}.IDs = {'listeningLab_L'}; + room{end}.room_dim = [7 4 5]; + room{end}.aparams = '0 1.5 0'; + room{end}.referenceRoomIdx = 1; + room{end}.t60Offset = 0.04; % In seconds + + % =============== Qualcomm: Lounge with doors open =============== + room{end + 1}.rname = ''; + room{end}.IDs = {'lounge_L_doorsclosed'}; + room{end}.room_dim = [6 7 3]; + room{end}.aparams = '0 1.5 0'; + room{end}.referenceRoomIdx = 1; + room{end}.t60Offset = 0.04; % In seconds + + % =============== Qualcomm: Lounge with doors closed =============== + room{end + 1}.rname = ''; + room{end}.IDs = {'lounge_L_doorsopen'}; + room{end}.room_dim = [6 7 3]; + room{end}.aparams = '0 1.5 0'; + room{end}.referenceRoomIdx = 1; + room{end}.t60Offset = 0.04; % In seconds + + % ------------------------------------------------------------------------- + + for i = 1:length(room) + % Add config data to the room definition, unless it is overridden by the + % room. + if (~isfield(room{i}, 'irRoot')) + room{i}.irRoot = irRoot; + end + if (~isfield(room{i}, 'predelayScheme')) + room{i}.predelayScheme = predelayScheme; + end + if (~isfield(room{i}, 'file_filters')) + room{i}.file_filters = [room{i}.IDs, fileFilters]; + end + if (~isfield(room{i}, 'sim_delay')) + room{i}.sim_delay = sim_delay; + end + if (~isfield(room{i}, 'rirRefDistance')) + room{i}.rirRefDistance = rirRefDistance; + end + if (~isfield(room{i}, 'srcMicDistance')) + room{i}.srcMicDistance = srcMicDistance; + end + if (~isfield(room{i}, 'tTruncate')) + room{i}.tTruncate = tTruncate; + end + + % Set figure offset to allow generating a set of figures for each room + figureOffset = i; + + % Calculate predelay, RT60 and DSR. + [predelay, rt60, dsr, fc, refData] = calcPredelayRT60DSR(room{i}, fs, predelayMethod, figureOffset, plotVerbosity); + + room{i}.predelay = predelay; + room{i}.RT60 = rt60; + room{i}.DSR = dsr; + room{i}.fc = fc; + room{i}.refData = refData; + + saveRenderConfigData(rendCfgPath, strjoin(room{i}.IDs, '-'), room{i}) + + end + + description = ''; + for i = 1:length(fileFilters) + description = strcat(description, fileFilters{i}); + if (i < length(fileFilters)) + description = strcat(description, '_'); + end + end + + logFilename = [irRoot '\ReverbAnalysisLog_' description '.mat']; + save(logFilename, 'room'); + +end + +% ========================================================================= +% Saves the acoustic environment data to IVAS compatible text format file +% Note: only single acoustic environment and single frequency grid supported + +function saveRenderConfigData(path, roomName, roomData) + indent = " "; + fid = fopen(fullfile(path, strcat(roomName, ".cfg")), "wt"); + + % Header + fprintf(fid, "[roomAcoustics]\n"); + fprintf(fid, "frequencyGridCount = 1;\n"); + fprintf(fid, "acousticEnvironmentCount = 1;\n\n"); + + % Frequency grid + fprintf(fid, "[frequencyGrid:0]\n"); + fprintf(fid, "method = individualFrequencies;\n"); + fprintf(fid, "nrBands = %d;\n", length(roomData.fc)); + fprintf(fid, "frequencies = [ %s ];\n", regexprep(strjoin(compose('%f', roomData.fc), ', '), '(.*?, ){4}', '$0\n' + indent)); + + % Acoustics environment + fprintf(fid, "[acousticEnvironment:0]\n"); + fprintf(fid, "frequencyGridIndex = 0;\n"); + fprintf(fid, "predelay = %f;\n", roomData.predelay); + fprintf(fid, "rt60 = [ %s ];\n", regexprep(strjoin(compose("%f", roomData.RT60), ', '), '(.*?, ){4}', '$0\n' + indent)); + fprintf(fid, "dsr = [ %s ];\n", regexprep(strjoin(compose("%e", roomData.DSR), ', '), '(.*?, ){4}', '$0\n' + indent)); + fclose(fid); +end +