classdef itaRavenProject < handle %RavenProject - The class for working with RAVEN. % This class allows you to create configurations of settings, to % run simulations and to retrieve the results. % % Using: % rpf = itaRavenProject(raven_project_file) % % Public Properties: % Enter this command to get the properties: % >> properties itaRavenProject % % Public Methods: % Enter this command to get the methods: % >> methods itaRavenProject % % Enter this command to get more info of method: % >> help itaRavenProject/methodname % % % % Examples: % rpf = itaRavenProject(); % rpf.SetModel('cave.ac'); % rpf.getEquationBasedReverbTime(); returns reverb time based on % sabine equation % rpf.run(); run simulation % rpf.getT30(); get T30 based on simulation % rpf.plotMaterialsAbsorption(); plots absorption coefficients % % Author: Soenke Pelzer (spe@akustik.rwth-aachen.de) % Lukas Aspöck (las@akustik.rwth-aachen.de) % % Version: 0.2 % First release: 01.11.10 % Last revision: 11.07.18 % Copyright: Institute of Technical Acoustics, RWTH Aachen University % % % This file is part of the application Raven for the ITA-Toolbox. All rights reserved. % You can find the license for this m-file in the application folder. % properties(Constant) % CONSTANTS freqLabel3rd = { ' 20 Hz', ' 25 Hz', ' 31 Hz', ' 40 Hz', ' 50 Hz', ' 63 Hz', ' 80 Hz', ' 100 Hz', ' 125 Hz', ' 160 Hz', ' 200 Hz', ... ' 250 Hz', ' 315 Hz', ' 400 Hz', ' 500 Hz', ' 630 Hz', ' 800 Hz', ' 1 kHz', '1.25 kHz', ' 1.6 kHz', ' 2 kHz', ... ' 2.5 kHz', '3.15 kHz', ' 4 kHz', ' 5 kHz', ' 6.3 kHz', ' 8 kHz', ' 10 kHz', '12.5 kHz', ' 16 kHz', ' 20 kHz'}; freqLabelOct = { ' 31 Hz', ' 63 Hz', ' 125 Hz', ' 250 Hz', ' 500 Hz', ' 1 kHz', ' 2 kHz', ' 4 kHz', ' 8 kHz', ' 16 kHz'}; freqVector3rd = [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 8000 10000 12500 16000 20000]; freqVectorOct = [31.5 63 125 250 500 1000 2000 4000 8000 16000]; MODE_BSP = 0; MODE_HASH = 1; MODE_BRUTE = 2; COORD_TRAFO_SKETCHUP2RAVEN = [1 3 -2]; COORD_TRAFO_RAVEN2SKETCHUP = [1 -3 2]; end properties (GetAccess = 'public', SetAccess = 'private') % raven ravenExe ravenLogFile = 'RavenLog.txt' ravenProjectFile ravenIniFile projectName % general sampleRate = 44100 % paths pathResults pathDirectivities pathMaterials fileHRTF fileSpeakers % model modelFileList = [] model = [] % [Global] simulationTypeIS simulationTypeRT generateRIR generateBRIR generateISHOA generateRTHOA generateISVBAP generateRTVBAP exportFilter exportHistogram exportWallHitLog exportPlaneWaveList accelerationType logPerformance % [PrimarySources] % sourceNames sourceDirectivity sourcePositions sourceViewVectors sourceUpVectors sourceSoundStates sourceSoundLevels % [Receiver] % receiverNames receiverPositions receiverViewVectors receiverUpVectors receiverStates % [ImageSources] % ISOrder_PS ISOrder_SS ISSkipDirectSound % [RayTracing] % numParticles_Sphere numParticles_Portal energyLoss_Sphere energyLoss_Portal detectionSphereAziResolution detectionSphereEleResolution filterLength timeSlotLength radiusSphere fixReflectionPattern % [Filter] % fixPoissonSequence poissonSequenceNumber filterResolution maxReflectionDensity ambisonicsOrder numberSpreadedSources spreadingStdDeviation fftDegreeForWallFilterInterpolation % [PlaneWaveLists] % planeWaveList_IS = [] planeWaveList_RT = [] % [Performance] performance = struct('ISFilterMonaural',[],'ISFilterBinaural',[],'RTFilterMonaural',[],'RTFilterBinaural',[],'ISGenerateImageSources',[],'ISTransformationMatrix',[],'ISAudibilityTest',[],'RTTotal',[],'RTBands',[]); end properties (GetAccess = 'public', SetAccess = 'public') % [WallHitLogs] % wallHitLog_IS = [] wallHitLog_RT = [] initialParticleEnergy = [] end properties (GetAccess = 'private', SetAccess = 'private') ravenExe64 = '..\bin64\RavenConsole64.exe' ravenExe32 = '..\bin32\RavenConsole.exe' rpf_ini raven_ini projectID projectTag projectLoaded = false simulationDone = false keepOutputFiles plotModelHandle = []; % [PrimarySources] % sourceNameString sourceDirectivityString % [Receiver] % receiverNameString uniformReceiverGridX = [] uniformReceiverGridY = [] uniformReceiverGridZ = [] % RESULTS % monoIR = [] monoIR_IS = [] monoIR_RT = [] binauralIR = [] binauralIR_IS = [] binauralIR_RT = [] ambisonicsIR = [] ambisonicsIR_IS = [] ambisonicsIR_RT = [] vbapIR = [] vbapIR_IS = [] vbapIR_RT = [] histogram = [] histogramRT = [] end %---------------------- PUBLIC METHODS -------------------------------% methods %------------------------------------------------------------------ function obj = itaRavenProject(raven_project_file) %RavenProject - constructor % To Create a new project with empty default configuration. % % Using: % rpf = RavenProject(raven_project_file) % % Input: % [optional] existing raven project file % % Output: % obj - an instance of class RavenProject % if strcmp(computer('arch'), 'win32') obj.ravenExe = obj.ravenExe32; elseif strcmp(computer('arch'), 'win64') obj.ravenExe = obj.ravenExe64; else error('Only Windows OS are supported.'); end itaRavenProjectPath = which('itaRavenProject.m'); obj.ravenIniFile = [ itaRavenProjectPath(1:end-9) '.ini']; ravenIniExists = exist(obj.ravenIniFile,'file'); if (ravenIniExists) % load path from itaRaven.ini obj.raven_ini = IniConfig(); obj.raven_ini.ReadFile(obj.ravenIniFile); obj.ravenExe = obj.raven_ini.GetValues('Global', 'PathRavenExe', obj.ravenExe); if (ischar(obj.ravenExe) && exist(obj.ravenExe,'file')) obj.raven_ini.WriteFile(obj.ravenIniFile); end end if (~exist(obj.ravenExe,'file')) % neither the default relative raven console path nor the path in % itaRaven.ini was found, try to locate RavenConsole locatedRavenExe = which(obj.ravenExe(10:end)); % try default raven exe path after instalation if isempty(locatedRavenExe) if strcmp(computer('arch'), 'win32') defaultInstallationPathRavenExe = 'C:\ITASoftware\Raven\bin32\RavenConsole.exe'; elseif strcmp(computer('arch'), 'win64') defaultInstallationPathRavenExe = 'C:\ITASoftware\Raven\bin64\RavenConsole64.exe'; else error('Only Windows OS are supported.'); end if exist(defaultInstallationPathRavenExe,'file') locatedRavenExe = defaultInstallationPathRavenExe; obj.ravenExe = defaultInstallationPathRavenExe; end end if isempty(locatedRavenExe) disp('[itaRaven]:'); disp('No RAVEN binary was found! Please select path to the RAVEN console application (RavenConsole.exe/RavenConsole64.exe)!'); disp('To run RAVEN simulations, an installation of the RAVEN software is required. '); disp('Individual licenses for academic purposes are available on request.'); disp('Please contact: las@akustik.rwth-aachen.de'); [ selectedRavenExe, selectedRavenPath] = uigetfile('*.exe',' No RAVEN binary was found! Please select path to RavenConsole.exe'); if (ischar(selectedRavenExe)&& ischar(selectedRavenPath)) obj.ravenExe = [ selectedRavenPath selectedRavenExe]; else warning('WARNING: RAVEN binary path was set to default, but binary was not found. Please specify the correct path in itaRaven.ini'); obj.ravenExe = defaultInstallationPathRavenExe; end end if (~ravenIniExists) obj.raven_ini = IniConfig(); obj.raven_ini.AddSections({'Global'}); obj.raven_ini.AddKeys('Global', {'PathRavenExe'}, {obj.ravenExe}); obj.raven_ini.WriteFile(obj.ravenIniFile); else obj.raven_ini.SetValues('Global', {'PathRavenExe'}, {obj.ravenExe}); obj.raven_ini.WriteFile(obj.ravenIniFile); end end % check if raven project file exists if (nargin > 0) && exist(raven_project_file, 'file') obj.loadRavenConfig(raven_project_file); else error('No raven project file given or file not found.'); end end function delete(obj) obj.deleteResultsInRavenFolder(); end %------------------------------------------------------------------ function setRavenIniPath(obj, newPath) obj.ravenIniFile = newPath; end %------------------------------------------------------------------ function copyProjectToNewRPFFile(obj, newPath) if (exist(fileparts(newPath),'dir')) obj.ravenProjectFile = newPath; obj.rpf_ini.WriteFile(obj.ravenProjectFile); else error('itaRavenProject/copyProjectToNewRPFFile: invalid path'); end end %------------------------------------------------------------------ function setRavenExe(obj, newRavenExe) if (exist(newRavenExe,'file')) obj.ravenExe = newRavenExe; if (exist(obj.ravenIniFile,'file')) obj.raven_ini.SetValues('Global', {'PathRavenExe'}, {obj.ravenExe}); obj.raven_ini.WriteFile(obj.ravenIniFile); else obj.raven_ini = IniConfig(); obj.raven_ini.AddSections({'Global'}); obj.raven_ini.AddKeys('Global', {'PathRavenExe'}, {obj.ravenExe}); end else error('[itaRaven]: Error: Path to new Raven binary not found!'); end end %------------------------------------------------------------------ function loadRavenConfig(obj, filename) %loadRavenConfig - Reads an existing raven project file % % change relative to absolute path if (~strcmp(filename(2),':')) obj.ravenProjectFile = [pwd '\' filename]; else obj.ravenProjectFile = filename; end obj.rpf_ini = IniConfig(); obj.rpf_ini.ReadFile(filename); % [Global] % obj.projectName = obj.rpf_ini.GetValues('Global', 'ProjectName', 'Matlab'); obj.projectTag = obj.projectName; obj.pathResults = obj.rpf_ini.GetValues('Global', 'ProjectPath_Output', '..\RavenOutput'); obj.pathDirectivities = obj.rpf_ini.GetValues('Global', 'ProjectPath_DirectivityDB', '..\RavenDatabase\DirectivityDatabase'); obj.pathMaterials = obj.rpf_ini.GetValues('Global', 'ProjectPath_MaterialDB', '..\RavenDatabase\MaterialDatabase'); obj.fileHRTF = obj.rpf_ini.GetValues('Global', 'ProjectPath_HRTFDB', '..\RavenDatabase\HRTF\ITA-Kunstkopf_HRIR_AP11_Pressure_Equalized_3x3_256.daff'); obj.fileSpeakers = obj.rpf_ini.GetValues('Global', 'SpeakerConfigFile', 'Speakers.ini'); obj.simulationTypeIS = obj.rpf_ini.GetValues('Global', 'simulationTypeIS', 1); obj.simulationTypeRT = obj.rpf_ini.GetValues('Global', 'simulationTypeRT', 1); obj.generateRIR = obj.rpf_ini.GetValues('Global', 'generateRIR', 1); obj.generateBRIR = obj.rpf_ini.GetValues('Global', 'generateBRIR', 1); obj.generateISHOA = obj.rpf_ini.GetValues('Global', 'generateISHOA', 0); obj.generateRTHOA = obj.rpf_ini.GetValues('Global', 'generateRTHOA', 0); obj.generateISVBAP = obj.rpf_ini.GetValues('Global', 'generateISVBAP', 0); obj.generateRTVBAP = obj.rpf_ini.GetValues('Global', 'generateRTVBAP', 0); obj.exportFilter = obj.rpf_ini.GetValues('Global', 'exportFilter', 1); obj.exportHistogram = obj.rpf_ini.GetValues('Global', 'exportHistogram', 1); obj.exportWallHitLog = obj.rpf_ini.GetValues('Global', 'exportWallHitList', 0); obj.exportPlaneWaveList = obj.rpf_ini.GetValues('Global', 'exportPlaneWaveList', 0); obj.accelerationType = obj.rpf_ini.GetValues('Global', 'accelerationType', 0); % default 0 = MODE_BSP obj.logPerformance = obj.rpf_ini.GetValues('Global', 'logPerformance', 0); obj.keepOutputFiles = obj.rpf_ini.GetValues('Global', 'keepOutputFiles', 0); % change relative to absolute paths if obj.ravenExe(2) == ':' % absolute path ravenBasePath = fileparts(fileparts(obj.ravenExe)); % base path of raven if (strcmp(obj.pathResults(1:2),'..')), obj.pathResults = [ ravenBasePath obj.pathResults(3:end) ]; end if (strcmp(obj.pathDirectivities(1:2),'..')), obj.pathDirectivities = [ ravenBasePath obj.pathDirectivities(3:end) ]; end if (strcmp(obj.pathMaterials(1:2),'..')), obj.pathMaterials = [ ravenBasePath obj.pathMaterials(3:end) ]; end if (strcmp(obj.fileHRTF(1:2),'..')), obj.fileHRTF = [ ravenBasePath obj.fileHRTF(3:end) ]; end end % [Rooms] % model_string = obj.rpf_ini.GetValues('Rooms', 'Model'); obj.modelFileList = textscan(model_string, '%s', 'Delimiter' , ','); obj.modelFileList = obj.modelFileList{1}; % textscan implementation issue if numel(obj.modelFileList) == 1 obj.modelFileList = obj.modelFileList{1}; % de-cell if only 1 room given if obj.ravenExe(2) == ':' % convert to absolute path if (strcmp(obj.modelFileList(1:2),'..')), obj.modelFileList = [ ravenBasePath obj.modelFileList(3:end) ]; end end end % [PrimarySources] % obj.sourceNameString = obj.rpf_ini.GetValues('PrimarySources', 'sourceNames', 'Sender'); obj.sourceNames = textscan(obj.sourceNameString, '%s', 'Delimiter', ','); obj.sourceNames = obj.sourceNames{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht obj.sourceDirectivityString = obj.rpf_ini.GetValues('PrimarySources', 'sourceDirectivity', ''); if ~isempty(obj.sourceDirectivityString) obj.sourceDirectivity = textscan(obj.sourceDirectivityString, '%s', 'Delimiter', ','); obj.sourceDirectivity = obj.sourceDirectivity{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht else obj.sourceDirectivity = {}; end obj.sourcePositions = obj.rpf_ini.GetValues('PrimarySources', 'sourcePositions'); obj.sourcePositions = reshape(obj.sourcePositions, 3, numel(obj.sourcePositions)/3)'; obj.sourceViewVectors = obj.rpf_ini.GetValues('PrimarySources', 'sourceViewVectors', '0, 0, 1'); obj.sourceViewVectors = reshape(obj.sourceViewVectors, 3, numel(obj.sourceViewVectors)/3)'; obj.sourceUpVectors = obj.rpf_ini.GetValues('PrimarySources', 'sourceUpVectors', '0, 1, 0'); obj.sourceUpVectors = reshape(obj.sourceUpVectors, 3, numel(obj.sourceUpVectors)/3)'; obj.sourceSoundStates = obj.rpf_ini.GetValues('PrimarySources', 'sourceSoundStates', '1'); obj.sourceSoundLevels = obj.rpf_ini.GetValues('PrimarySources', 'sourceSoundLevels', '100'); % [Receiver] % obj.receiverNameString = obj.rpf_ini.GetValues('Receiver', 'receiverNames', 'Receiver'); obj.receiverNames = textscan(obj.receiverNameString, '%s', 'Delimiter', ','); obj.receiverNames = obj.receiverNames{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht obj.receiverPositions = obj.rpf_ini.GetValues('Receiver', 'receiverPositions'); obj.receiverPositions = reshape(obj.receiverPositions, 3, numel(obj.receiverPositions)/3)'; obj.receiverViewVectors = obj.rpf_ini.GetValues('Receiver', 'receiverViewVectors', '1, 0 ,0'); obj.receiverViewVectors = reshape(obj.receiverViewVectors, 3, numel(obj.receiverViewVectors)/3)'; obj.receiverUpVectors = obj.rpf_ini.GetValues('Receiver', 'receiverUpVectors', '0, 1, 0'); obj.receiverUpVectors = reshape(obj.receiverUpVectors, 3, numel(obj.receiverUpVectors)/3)'; % [ImageSources] % obj.ISOrder_PS = obj.rpf_ini.GetValues('ImageSources', 'ISOrder_PS', 2); obj.ISOrder_SS = obj.rpf_ini.GetValues('ImageSources', 'ISOrder_SS', 2); obj.ISSkipDirectSound = obj.rpf_ini.GetValues('ImageSources', 'ISSkipDirectSound', 0); % [RayTracing] % obj.numParticles_Sphere = obj.rpf_ini.GetValues('RayTracing', 'numberOfParticles_DetectionSphere', 10000); obj.numParticles_Portal = obj.rpf_ini.GetValues('RayTracing', 'numberOfParticles_Portal', 10000); obj.energyLoss_Sphere = obj.rpf_ini.GetValues('RayTracing', 'energyLoss_DetectionSphere', 63); obj.energyLoss_Portal = obj.rpf_ini.GetValues('RayTracing', 'energyLoss_Portal', 63); obj.filterLength = obj.rpf_ini.GetValues('RayTracing', 'filterLength_DetectionSphere', 2000); obj.timeSlotLength = obj.rpf_ini.GetValues('RayTracing', 'timeResolution_DetectionSphere', 6); obj.radiusSphere = obj.rpf_ini.GetValues('RayTracing', 'radius_DetectionSphere', 0.5); obj.fixReflectionPattern= obj.rpf_ini.GetValues('RayTracing', 'fixReflectionPattern', 0); % [Filter] % obj.fixPoissonSequence = obj.rpf_ini.GetValues('Filter', 'setFixPoissonSequence', 0); obj.poissonSequenceNumber = obj.rpf_ini.GetValues('Filter', 'poissonSequenceNumber', 9876); obj.filterResolution = obj.rpf_ini.GetValues('Filter', 'filterResolution', 1); % 1 = Octave, 0 = 3rd Octave obj.ambisonicsOrder = obj.rpf_ini.GetValues('Filter', 'ambisonicsOrder', -1); obj.numberSpreadedSources = obj.rpf_ini.GetValues('Filter', 'numberSpreadedSources', 0); obj.spreadingStdDeviation = obj.rpf_ini.GetValues('Filter', 'spreadingStdDeviation', 0.2); obj.fftDegreeForWallFilterInterpolation = obj.rpf_ini.GetValues('Filter', 'fftDegreeForWallFilterInterpolation', 8); %% check input data files (HRTF, source directivity, room model) % HRTF (also shows if DAFF file is v15 or v17) if (exist(obj.fileHRTF,'file') == 0) warning('HRTF file not found at the specified location.'); else if (contains(obj.fileHRTF,'.v17.daff')) disp('itaRavenProject: Loaded DAFFv17 HRTF file'); end if (contains(obj.fileHRTF,'.v15.daff')) disp('itaRavenProject: Loaded DAFFv15 HRTF file'); end end % Source directivity (only checks first directivity file) if (~isempty(obj.sourceDirectivity)) % first check if any directivity was defined if (exist([ obj.pathDirectivities '\' obj.sourceDirectivity{1}],'file') == 0) warning('Source directivity file not found at the specified location.'); else if (contains(obj.sourceDirectivity{1},'.v17.daff')) disp('itaRavenProject: Loaded DAFFv17 source directivity file'); end if (contains(obj.sourceDirectivity{1},'.v15.daff')) disp('itaRavenProject: Loaded DAFFv15 source directivity file'); end end end % room model if (numel(obj.modelFileList) == 0 ) warning('No room model file (*.ac) specified in project file.'); end if (numel(obj.modelFileList) == 1 && exist(obj.modelFileList,'file') == 0) warning('Room model file (*.ac) not found at the specified location.'); end obj.projectLoaded = true; end %------------------------------------------------------------------ function reloadRavenConfig(obj, filename) if nargin < 2 filename = obj.ravenProjectFile; end obj.loadRavenConfig(filename); end %------------------------------------------------------------------ function saveRavenConfig(obj, filename) %saveRavenConfig - Saves the current object in a (different) %raven project file. Can be used as a log if various %simulations are run with multiple configurations % obj.rpf_ini.WriteFile(filename); end %------------------------------------------------------------------ function run(obj) if numel(obj.projectName) > 70 warning('Long project name. This might cause an error while writing output files. Consider resetting your project name!') end obj.simulationDone = false; if obj.projectLoaded savedProjectName = obj.projectName; obj.projectID = datestr(now, 30); obj.projectTag = [obj.projectName obj.projectID]; % give the project name a date and time string to help to identify the results obj.setProjectName(obj.projectTag); % run the simulation disp(['Running simulation... (' obj.ravenExe ')']); if exist(obj.ravenLogFile, 'file') delete(obj.ravenLogFile); end % system([obj.ravenExe ' "' obj.ravenProjectFile '" >> ' obj.ravenLogFile]); if (~exist(obj.ravenExe,'file')) error('[itaRaven]: Error: Cannot find Raven binary file! Please check itaRaven.ini file!'); end prevPath = pwd; cd(fileparts(obj.ravenExe)); dos(['"' obj.ravenExe '"' ' "' obj.ravenProjectFile '"'],'-echo'); cd(prevPath); % restore the initial project name obj.setProjectName(savedProjectName); % gather results disp('[R] Simulation seems to be finished. Getting results...'); obj.gatherResults(); disp('[R] Done.'); obj.simulationDone = true; % delete results in raven folder structure -> they are copied now into this class if (obj.keepOutputFiles == 0) obj.deleteResultsInRavenFolder(); end else disp('No projected defined yet.'); end end %------------------------------------------------------------------ function runNoGathering(obj) % basicly the same as the run method, but without the call of % obj.gatherResults(); obj.simulationDone = false; if obj.projectLoaded savedProjectName = obj.projectName; obj.projectID = datestr(now, 30); obj.projectTag = [obj.projectName obj.projectID]; % give the project name a date and time string to help to identify the results obj.setProjectName(obj.projectTag); % set filter length to the length of the reverberation % obj.setFilterLengthToReverbTime(); % run the simulation disp(['Running simulation... (' obj.ravenExe ')']); if exist(obj.ravenLogFile, 'file') delete(obj.ravenLogFile); end % system([obj.ravenExe ' "' obj.ravenProjectFile '" >> ' obj.ravenLogFile]); prevPath = pwd; cd(fileparts(obj.ravenExe)); dos(['"' obj.ravenExe '"' ' "' obj.ravenProjectFile '"'],'-echo'); disp('Done.'); cd(prevPath); % restore the initial project name obj.setProjectName(savedProjectName); % gather results disp('This function does _not_ gather the results. Please provide arguments to getWallHitLogBand(band)'); % obj.gatherResults(); % disp('Done.'); obj.simulationDone = true; % delete results in raven folder structure -> they are copied now into this class % if (obj.keepOutputFiles == 0) % obj.deleteResultsInRavenFolder(); % end else disp('No projected defined yet.'); end end %------------------------------------------------------------------ function keepImpulseResponseFiles(obj, keepFiles) % keepImpulseResponseFiles % % By default, RAVEN Impulse Responses are deleted from the hard disk % after simulation and are only available in your rpf project % % By setting keepOutFiles to 1 / true, results are kept in % the Output-Folder (obj.pathResults) obj.keepOutputFiles = keepFiles; obj.rpf_ini.SetValues('Global', 'keepOutputFiles', keepFiles); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function openOutputFolder(obj) % opens the output folder in windows explorer if (exist(obj.pathResults,'dir')) dos(['C:\Windows\Explorer.exe ' obj.pathResults]); else disp('Output Folder does not exist!'); end end %------------------------------------------------------------------ function numReceivers = createReceiverArray(obj, xpositions, zpositions, yheight) numReceivers = numel(xpositions) * numel(zpositions); [obj.uniformReceiverGridX, obj.uniformReceiverGridZ] = meshgrid(xpositions, zpositions); obj.uniformReceiverGridY = ones(size(obj.uniformReceiverGridX,1), size(obj.uniformReceiverGridX,2)) * yheight; rec_names = obj.makeNumberedNames('Array', numReceivers); rec_states = str2num(obj.writeXtimes_num(1, numReceivers)); rec_view = str2num(obj.writeXtimes('1,0,0', numReceivers)); rec_up = str2num(obj.writeXtimes('0,1,0', numReceivers)); positions = zeros(numReceivers * 3, 1); positions(1:3:end) = obj.uniformReceiverGridX(:); positions(2:3:end) = obj.uniformReceiverGridY(:); positions(3:3:end) = obj.uniformReceiverGridZ(:); obj.setReceiverNames(rec_names); obj.setReceiverPositions(positions); obj.setReceiverViewVectors(rec_view); obj.setReceiverUpVectors(rec_up); obj.setReceiverStates(rec_states); disp(['Receiver grid created successfully. ' num2str(numReceivers) ' receivers placed.']); end %------------------------------------------------------------------ function createReceiverArrayAuto(obj, yheight, distance, roomID) if nargin < 4 roomID = 0; % default = first room (ID = 0) end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end if nargin < 3 distance = 1; % default = 1 meter end if nargin < 2 yheight = min(roommodel.nodes(:,2)) + 1.2; % default listener 1.2 meters above ground if yheight > max(roommodel.nodes(:,2)) yheight = mean(roommodel.nodes(:,2)); % if 1.2 meters is out of ceiling, try putting receivers in the middle between floor and ceiling end end xpositions = min(roommodel.nodes(:,1)) + distance/2 : distance : max(roommodel.nodes(:,1)) - distance/2; zpositions = min(roommodel.nodes(:,3)) + distance/2 : distance : max(roommodel.nodes(:,3)) - distance/2; obj.createReceiverArray(xpositions, zpositions, yheight); end %------------------------------------------------------------------ function rebuildReceiverGrid(obj) numReceivers = numel(obj.receiverPositions) / 3; %xyz if numReceivers < 1 error('No receiver position data present.'); end x = obj.receiverPositions(:,1); y = obj.receiverPositions(:,2); z = obj.receiverPositions(:,3); averageDistance = sqrt((max(x)-min(x))*(max(z)-min(z)) / numReceivers); obj.createReceiverArray(min(x) : averageDistance : max(x), min(z) : averageDistance : max(z), mean(y)); end %------------------------------------------------------------------ function setModel(obj, filename) if (nargin < 2) error('Not enough input arguments.') end if ~ischar(filename) error('Requires string input.') else if iscell(filename) obj.modelFileList = obj.cat_cell_of_strings(filename); else obj.modelFileList = filename; end obj.rpf_ini.SetValues('Rooms', 'Model', obj.modelFileList); end obj.model = []; obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function [outputFilePath] = setModelToShoebox(obj, length,width,height) % setModelToShoebox % creates a shoebox room model and automatically sets the model % of the current project to the defined shoebox room. % The shoebox room has 6 materials assigned (matShoebox1 .. % matShoebox6), which can be set using % rpf.setMaterial(matShoebox1,abs,scat) % % surface order: % (1) floor, (2) ceiling, % (3) larger wall (length x height; left, view from origin) % (4) smaller wall (width x height; front) % (5) larger wall (length x height; right) % (6) smaller wall (width x height; back) % % Using: % outputFileName = rpf.setModelToShoebox(myLength,myWidth,myHeight); % % see ita_raven_demo_shoebox for more details on usage % % Input: % length, width, height (all in meters) % % Output: % outputFilePath: full path of save ac3d model % outputFileName = ['Shoebox_' num2str(length) 'x' num2str(width) 'x' num2str(height) '.ac' ]; if strfind(obj.modelFileList,'RavenModels') outputFilePath = [ obj.modelFileList(1:strfind(obj.modelFileList,'RavenModels')+10) '\' outputFileName ]; else % if model wasn't stored in RavenModels folder, save it in % the folder of the current rpf file. outputFilePath = [ fileparts(obj.ravenProjectFile) '\' outputFileName ]; end fid = fopen(outputFilePath,'w'); fprintf(fid,'AC3Db\n'); for iMat=1:6 fprintf(fid,['MATERIAL "matShoebox' num2str(iMat) '" rgb ' num2str(0.0) ' ' num2str(0.0) ' ' num2str(0.1*iMat) ' amb 0.2 0.2 0.2 emis 0 0 0 spec 0.2 0.2 0.2 shi 128 trans 0 \n']); end fprintf(fid,'OBJECT world\n'); fprintf(fid,'kids 6\n'); lS=sprintf('%1.4f',length); wS=sprintf('%1.4f',-width); hS=sprintf('%1.4f',height); for iWall=1:6 fprintf(fid,'OBJECT poly\n'); fprintf(fid,'name "polygon_object"\n'); fprintf(fid,'numvert 4\n'); if iWall==1 fprintf(fid,'0 0 0\n'); fprintf(fid,[lS ' 0 0\n']); fprintf(fid,[lS ' 0 ' wS '\n']); fprintf(fid,['0 0 ' wS '\n']); end if iWall==2 fprintf(fid,[lS ' ' hS ' 0\n']); fprintf(fid,['0 ' hS ' 0\n']); fprintf(fid,['0 ' hS ' ' wS '\n']); fprintf(fid,[lS ' ' hS ' ' wS '\n']); end if iWall==3 fprintf(fid,[lS ' 0 0\n']); fprintf(fid,'0 0 0\n'); fprintf(fid,['0 ' hS ' 0\n']); fprintf(fid,[lS ' ' hS ' 0\n']); end if iWall==4 fprintf(fid,'0 0 0\n'); fprintf(fid,['0 0 ' wS '\n']); fprintf(fid,['0 ' hS ' ' wS '\n']); fprintf(fid,['0 ' hS ' 0\n']); end if iWall==5 fprintf(fid,['0 0 ' wS '\n']); fprintf(fid,[lS ' 0 ' wS '\n']); fprintf(fid,[lS ' ' hS ' ' wS '\n']); fprintf(fid,['0 ' hS ' ' wS '\n']); end if iWall==6 fprintf(fid,[lS ' 0 ' wS '\n']); fprintf(fid,[lS ' 0 0\n']); fprintf(fid,[lS ' ' hS ' 0\n']); fprintf(fid,[lS ' ' hS ' ' wS '\n']); end fprintf(fid,'numsurf 1\n'); fprintf(fid,'SURF 0x10\n'); fprintf(fid,['mat ' num2str(iWall-1) '\n']); fprintf(fid,'refs 4\n'); fprintf(fid,'3 0 0\n'); fprintf(fid,'2 0 0\n'); fprintf(fid,'1 0 0\n'); fprintf(fid,'0 0 0\n'); fprintf(fid,'kids 0\n'); end fclose(fid); obj.setModel(outputFilePath); end %------------------------------------------------------------------ function [outputFilePath] = setModelToFaces(obj,points,faces,materials) % setModelToFaces % creates a room model using points and faces and automatically % sets the model of the current project to the so defined room. % The room has one material assigned to each face (mat1 .. % mat6), which can be set using % rpf.setMaterial(mat1,abs,scat) % % Using: % outputFileName = rpf.setModelToFaces(points,faces,materials); % % see ita_raven_demo_faces_room for more details on usage % % Input: % points (matrix Nx3 in meters, columns represnt x, y, and z) % faces (cell array, faces are defined pointing to the rows % of matrix points) the first element points the material % materials (cell array) This array does not necesarlly match % the number of surfaces % Output: % FilePath: full path of save ac3d model % nW=length(faces); outputFileName = ['room_' num2str(nW) '_faces']; if strfind(obj.modelFileList,'RavenModels') outputFilePath = [ obj.modelFileList(1:strfind(obj.modelFileList,'RavenModels')+10) '\' outputFileName ]; else % if model wasn't stored in RavenModels folder, save it in % the folder of the current rpf file. outputFilePath = [ fileparts(obj.ravenProjectFile) '\' outputFileName ]; end fn = floor(size(points,2)/4); % first element of a three elements point if size(points,2)>3 idx = (points(:,1)); else idx=1:size(points,1); end fid = fopen(outputFilePath,'w'); fprintf(fid,'AC3Db\n'); for iW=1:nW fprintf(fid,['MATERIAL "' materials{faces{iW}(1)} '" rgb ' num2str(0.0) ' ' num2str(0.0) ' ' num2str(0.9/nW*iW) ' amb 0.2 0.2 0.2 emis 0 0 0 spec 0.2 0.2 0.2 shi 128 trans 0 \n']); end fprintf(fid,'OBJECT world\n'); fprintf(fid,['kids ' num2str(nW) '\n']); for iW=1:nW nP=length(faces{iW})-1; fprintf(fid,'OBJECT poly\n'); fprintf(fid,'name "polygon_object"\n'); fprintf(fid,['numvert ' num2str(nP) '\n']); for iP=1:nP fprintf(fid,[sprintf('%1.4f %1.4f %1.4f',points(faces{iW}(iP+1)==idx,fn+(1:3))) '\n']); end fprintf(fid,'numsurf 1\n'); fprintf(fid,'SURF 0x10\n'); fprintf(fid,['mat ' num2str(iW-1) '\n']); fprintf(fid,'refs 4\n'); fprintf(fid,'3 0 0\n'); fprintf(fid,'2 0 0\n'); fprintf(fid,'1 0 0\n'); fprintf(fid,'0 0 0\n'); fprintf(fid,'kids 0\n'); end fclose(fid); obj.setModel(outputFilePath); end %------------------------------------------------------------------ function figureHandle=plotModel(obj, tgtAxes, comp2axesMapping, wireframe) if isempty(obj.modelFileList) return; end if nargin < 4 wireframe = 0; else if ischar(wireframe) wireframe = isequal(wireframe, 'wireframe'); end end if nargin < 3 comp2axesMapping = [1 -3 2]; end if nargin < 2 if (ishandle(obj.plotModelHandle)) tgtAxes = obj.plotModelHandle; else figure; tgtAxes = gca; end end obj.plotModelHandle = tgtAxes; if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); obj.model{iRoom}.plotModel(tgtAxes, comp2axesMapping, wireframe); hold on; end else obj.model = itaAc3dModel(obj.modelFileList); obj.model.plotModel(tgtAxes, comp2axesMapping, wireframe); end else if iscell(obj.model) for iRoom = 1 : numel(obj.model) obj.model{iRoom}.plotModel(tgtAxes, comp2axesMapping, wireframe); hold on; end else obj.model.plotModel(tgtAxes, comp2axesMapping, wireframe); end end % get and transform source and receiver data invertAxes = sign(comp2axesMapping); comp2axesMapping = abs(comp2axesMapping); spos = obj.getSourcePosition; spos = spos(:, comp2axesMapping).*repmat(invertAxes,size(spos,1),1); sview = obj.getSourceViewVectors; sview = sview(:, comp2axesMapping).*repmat(invertAxes,size(sview,1),1); sup = obj.getSourceUpVectors; sup = sup(:, comp2axesMapping).*repmat(invertAxes,size(sup,1),1); snames = obj.getSourceNames; rpos = obj.getReceiverPosition; rpos = rpos(:, comp2axesMapping).*repmat(invertAxes,size(rpos,1),1); rview = obj.getReceiverViewVectors; rview = rview(:, comp2axesMapping).*repmat(invertAxes,size(rview,1),1); rup = obj.getReceiverUpVectors; rup = rup(:, comp2axesMapping).*repmat(invertAxes,size(rup,1),1); rnames = obj.getReceiverNames; % plot source and receivers plot3(spos(:,1),spos(:,2),spos(:,3),'marker','o','markersize',9,'linestyle','none','linewidth',1.5) % hold on; view(0,90); plot3(rpos(:,1),rpos(:,2),rpos(:,3),'marker','x','markersize',9,'linestyle','none','linewidth',1.5) % plot view vectors (red) of sources quiver3(spos(:,1),spos(:,2),spos(:,3),sview(:,1),sview(:,2),sview(:,3),0,'color','r','maxheadsize',1.5,'linewidth',1.5); % plot view/up vectors (red) of receivers quiver3(rpos(:,1),rpos(:,2),rpos(:,3),rview(:,1),rview(:,2),rview(:,3),0,'color','r','maxheadsize',1.5,'linewidth',1.5); % plot up vectors (green) of sources and receivers ) quiver3(spos(:,1),spos(:,2),spos(:,3),0.5*sup(:,1),0.5*sup(:,2),0.5*sup(:,3),0,'color','g','maxheadsize',1.5,'linewidth',1.5); quiver3(rpos(:,1),rpos(:,2),rpos(:,3),0.5*rup(:,1),0.5*rup(:,2),0.5*rup(:,3),0,'color','g','maxheadsize',1.5,'linewidth',1.5); % plot names text(spos(:,1)-0.2,spos(:,2),spos(:,3),snames,'FontName','Times','FontSize',20,'HorizontalAlignment','right') text(rpos(:,1)+0.2,rpos(:,2),rpos(:,3),rnames,'FontName','Times','FontSize',20) figureHandle=gca; end %------------------------------------------------------------------ function tgtAxes = plotModelRoom(obj, tgtAxes, comp2axesMapping, wireframe) % identical to plotModel, without sound sources and receivers if isempty(obj.modelFileList) return; end if nargin < 4 wireframe = 0; else if ischar(wireframe) wireframe = isequal(wireframe, 'wireframe'); end end if nargin < 3 comp2axesMapping = [1 -3 2]; end if nargin < 2 if (ishandle(obj.plotModelHandle)) tgtAxes = obj.plotModelHandle; else figure; tgtAxes = gca; end end obj.plotModelHandle = tgtAxes; if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); obj.model{iRoom}.plotModel(tgtAxes, comp2axesMapping, wireframe); hold on; end else obj.model = itaAc3dModel(obj.modelFileList); obj.model.plotModel(tgtAxes, comp2axesMapping, wireframe); end else if iscell(obj.model) for iRoom = 1 : numel(obj.model) obj.model{iRoom}.plotModel(tgtAxes, comp2axesMapping, wireframe); hold on; end else obj.model.plotModel(tgtAxes, comp2axesMapping, wireframe); end end end %------------------------------------------------------------------ function plotMaterialsAbsorption(obj, exportPlot, fileType, outputPath) if nargin < 2 exportPlot = false; end if nargin < 3 fileType='png'; end if nargin < 4 outputPath= obj.pathResults; end % strip final slash or backslash if ( strcmp(outputPath(end),'\') || strcmp(outputPath(end),'/') ) outputPath = outputPath(1:end-1); end freqVector = [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 8000 10000 12500 16000 20000]; freqLabel3rdVisual = { '', '', '31.5 Hz', '', '', '' '', ' ', ' 125 Hz', ' ', ' ', ... '', ' ', ' ', ' 500 Hz', ' ', ' ', '', '', ' ', ' 2 kHz', ... ' ', '', '', '', '', ' 8 kHz', ' ', '', '', '20 kHz'}; yticks = { '0.0','', '0.20','','0.40','','0.60','','0.80','','1.0'}; allMaterials = obj.getRoomMaterialNames; numberMaterials = length(allMaterials); currentMaterial = itaResult; currentMaterial.freqVector = freqVector; currentMaterial.freqData = []; for iMat=1:numberMaterials [absorp scatter ] = obj.getMaterial(allMaterials{iMat}); currentMaterial.freqData = [ currentMaterial.freqData absorp' ]; currentSurfaceArea = obj.getSurfaceAreaOfMaterial(allMaterials{iMat}); allMaterials{iMat} = strrep(allMaterials{iMat},'_',' '); allMaterials{iMat} = [ allMaterials{iMat} ' (S = ' num2str(currentSurfaceArea,'%5.2f') ' m² ;']; allMaterials{iMat} = [ allMaterials{iMat} ' A (Eyring, f=1000 Hz) = ' num2str(-currentSurfaceArea*log(1-currentMaterial.freqData(18,iMat)),'%5.2f') ' m² )']; end currentMaterial.channelNames = allMaterials; currentMaterial.allowDBPlot = false; currentPlot = ita_plot_freq(currentMaterial,'LineWidth',2); myAxes = gca; for iMat=1:numberMaterials myAxes.Children(iMat+2).Marker = 's'; end % change format of plot title(''); ylabel('Absorption coefficient'); xlabel('Frequency in Hz'); set(myAxes,'XLim',[20 20000]); set(myAxes,'YLim',[0 1.05]); leg = findobj(gcf,'Tag','legend'); set(leg,'Location','NorthWest'); set(leg,'FontSize',12); % remove [1] in legend entry for iMat=1:numberMaterials leg.String{iMat} = leg.String{iMat}(1:end-4); end % export plot to raven output if (exportPlot) [pathstr,name,ext] = fileparts(obj.ravenProjectFile); dateTimeStr = datestr(now,30); dateTimeStr = strrep(dateTimeStr,'T','_'); fileNamePNG = [ outputPath '\Absorption_' name '_' dateTimeStr '.png']; fileNamePDF = [ outputPath '\Absorption_' name '_' dateTimeStr '.pdf']; set(gcf,'units','normalized','outerposition',[0 0 1 1]); set(gcf,'PaperUnits','inches','PaperPosition',1.34*[0 0 8 5]); if strcmpi(fileType,'pdf') set(gcf,'PaperUnits','centimeters','PaperSize',[25.7, 16.5]); print('-dpdf','-r200', fileNamePDF); close(currentPlot); else if ~strcmpi(fileType,'png') warning('plotMaterialsAbsorption: Specified fileType not supported. Using PNG.'); end print('-dpng','-r200', fileNamePNG); end end end %------------------------------------------------------------------ function averageRoomScattering = plotMaterialsScattering(obj, exportPlot, fileType, outputPath) if nargin < 2 exportPlot = false; end if nargin < 3 fileType='png'; end if nargin < 4 outputPath= obj.pathResults; end % strip final slash or backslash if ( strcmp(outputPath(end),'\') || strcmp(outputPath(end),'/') ) outputPath = outputPath(1:end-1); end freqVector = [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 8000 10000 12500 16000 20000]; freqLabel3rdVisual = { '', '', '31.5 Hz', '', '', '' '', ' ', ' 125 Hz', ' ', ' ', ... '', ' ', ' ', ' 500 Hz', ' ', ' ', '', '', ' ', ' 2 kHz', ... ' ', '', '', '', '', ' 8 kHz', ' ', '', '', '20 kHz'}; yticks = { '0.0','', '0.20','','0.40','','0.60','','0.80','','1.0'}; allMaterials = obj.getRoomMaterialNames; numberMaterials = length(allMaterials); currentMaterial = itaResult; currentMaterial.freqVector = freqVector; currentMaterial.freqData = []; averageRoomScattering = zeros(1,31); for i=1:numberMaterials [absorp scatter ] = obj.getMaterial(allMaterials{i}); currentMaterial.freqData = [ currentMaterial.freqData scatter' ]; currentSurfaceArea = obj.getSurfaceAreaOfMaterial(allMaterials{i}); allMaterials{i} = strrep(allMaterials{i},'_',' '); allMaterials{i} = [ allMaterials{i} ' (S = ' num2str(currentSurfaceArea,'%5.2f') ' m² )']; averageRoomScattering = averageRoomScattering + (scatter)*currentSurfaceArea; end averageRoomScattering = averageRoomScattering / obj.getRoomSurfaceArea; currentMaterial.channelNames = allMaterials; currentMaterial.allowDBPlot = false; currentPlot=ita_plot_freq(currentMaterial,'LineWidth',2); myAxes = gca; for iMat=1:numberMaterials myAxes.Children(iMat+2).Marker = 's'; end % change format of plot title(''); ylabel('Scattering coefficient'); xlabel('Frequency in Hz'); set(myAxes,'XLim',[20 20000]); set(myAxes,'YLim',[0 1.05]); leg = findobj(gcf,'Tag','legend'); set(leg,'Location','NorthWest'); set(leg,'FontSize',12); % remove [1] in legend entry for iMat=1:numberMaterials leg.String{iMat} = leg.String{iMat}(1:end-4); end % export plot to raven output if (exportPlot) [pathstr,name,ext] = fileparts(obj.ravenProjectFile); dateTimeStr = datestr(now,30); dateTimeStr = strrep(dateTimeStr,'T','_'); fileNamePNG = [ outputPath '\Scattering_' name '_' dateTimeStr '.png']; fileNamePDF = [ outputPath '\Scattering_' name '_' dateTimeStr '.pdf']; set(gcf,'units','normalized','outerposition',[0 0 1 1]); set(gcf,'PaperUnits','inches','PaperPosition',1.34*[0 0 8 5]); if strcmpi(fileType,'pdf') set(gcf,'PaperUnits','centimeters','PaperSize',[25.7, 16.5]); print('-dpdf','-r200', fileNamePDF); close(currentPlot); else if ~strcmpi(fileType,'png') warning('plotMaterialsScattering: Specified fileType not supported. Using PNG.'); end print('-dpng','-r200', fileNamePNG); end end end % [Global] % %------------------------------------------------------------------ function setProjectName(obj, projectName) obj.projectName = projectName; obj.rpf_ini.SetValues('Global', 'ProjectName', projectName); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end function setPathMaterials(obj, pathMaterials) obj.pathMaterials = pathMaterials; obj.rpf_ini.SetValues('Global', 'ProjectPath_MaterialDB', pathMaterials); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSimulationTypeIS(obj, typeIS) obj.simulationTypeIS = typeIS; obj.rpf_ini.SetValues('Global', 'simulationTypeIS', typeIS); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSimulationTypeRT(obj, typeRT) obj.simulationTypeRT = typeRT; obj.rpf_ini.SetValues('Global', 'simulationTypeRT', typeRT); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setGenerateRIR(obj, genRIR) obj.generateRIR = genRIR; obj.rpf_ini.SetValues('Global', 'generateRIR', genRIR); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setGenerateBRIR(obj, genBRIR) obj.generateBRIR = genBRIR; obj.rpf_ini.SetValues('Global', 'generateBRIR', genBRIR); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setGenerateISHOA(obj, genISHOA) obj.generateISHOA = genISHOA; obj.rpf_ini.SetValues('Global', 'generateISHOA', genISHOA); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setGenerateRTHOA(obj, genRTHOA) obj.generateRTHOA = genRTHOA; obj.rpf_ini.SetValues('Global', 'generateRTHOA', genRTHOA); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setGenerateISVBAP(obj, genISVBAP) obj.generateISVBAP = genISVBAP; obj.rpf_ini.SetValues('Global', 'generateISVBAP', genISVBAP); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setGenerateRTVBAP(obj, genRTVBAP) obj.generateRTVBAP = genRTVBAP; obj.rpf_ini.SetValues('Global', 'generateRTVBAP', genRTVBAP); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setExportFilter(obj, exportFilter) obj.exportFilter = exportFilter; obj.rpf_ini.SetValues('Global', 'exportFilter', exportFilter); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setExportHistogram(obj, exportHisto) obj.exportHistogram = exportHisto; obj.rpf_ini.SetValues('Global', 'exportHistogram', exportHisto); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setExportWallHitLog(obj, exportWallHitLog) obj.exportWallHitLog = exportWallHitLog; obj.rpf_ini.SetValues('Global', 'exportWallHitList', exportWallHitLog); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setExportPlaneWaveList(obj, exportPlaneWaveList) obj.exportPlaneWaveList = exportPlaneWaveList; obj.rpf_ini.SetValues('Global', 'exportPlaneWaveList', exportPlaneWaveList); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setAccelerationType(obj, accelerationType) obj.accelerationType = accelerationType; obj.rpf_ini.SetValues('Global', 'AccelerationType', accelerationType); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setLogPerformance(obj, logPerformance) obj.logPerformance = logPerformance; obj.rpf_ini.SetValues('Global', 'logPerformance', logPerformance); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setPortals(obj, portalTypes, portalStatus) if nargin < 3 % if only 1 parameter is given (portalTypes), interprete % them as portal states! (more likely to be used like this) obj.rpf_ini.SetValues('Portals', 'portalStatus', obj.make_proper_string(portalTypes)); % 1 = open, 0 = closed else obj.rpf_ini.SetValues('Portals', 'portalType', portalTypes); obj.rpf_ini.SetValues('Portals', 'portalStatus', obj.make_proper_string(portalStatus)); % 1 = open, 0 = closed end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setTemperature(obj, temp) obj.rpf_ini.SetValues('Rooms', 'Temperature', num2str(temp)); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function temp = getTemperature(obj) temp = obj.rpf_ini.GetValues('Rooms', 'Temperature', -1); end %------------------------------------------------------------------ function setHumidity(obj, humid) obj.rpf_ini.SetValues('Rooms', 'Humidity', num2str(humid)); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function humid = getHumidity(obj) humid = obj.rpf_ini.GetValues('Rooms', 'Humidity', -1); end %------------------------------------------------------------------ function setPressure(obj, press) obj.rpf_ini.SetValues('Rooms', 'Pressure', num2str(press)); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function press = getPressure(obj) press = obj.rpf_ini.GetValues('Rooms', 'Pressure', -1); end %------------------------------------------------------------------ function enableAirAbsorption(obj) obj.rpf_ini.SetValues('Rooms', 'noAirAbsorption', 0); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function disableAirAbsorption(obj) obj.rpf_ini.SetValues('Rooms', 'noAirAbsorption', 1); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function airAbsEnabled = getAirAbsorptionEnabled(obj) airAbsEnabled = ~obj.rpf_ini.GetValues('Rooms', 'noAirAbsorption', 0); end % [PrimarySources] % %------------------------------------------------------------------ function setSourceNames(obj, sourceNames) if iscell(sourceNames) obj.sourceNameString = obj.cat_cell_of_strings(sourceNames); obj.sourceNames = sourceNames; else obj.sourceNameString = sourceNames; obj.sourceNames = textscan(obj.sourceNameString, '%s', 'Delimiter', ','); obj.sourceNames = obj.sourceNames{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht end obj.rpf_ini.SetValues('PrimarySources', 'sourceNames', obj.sourceNameString); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function directivityName = getSourceDirectivity(obj, srcID) readstring = obj.rpf_ini.GetValues('PrimarySources', 'sourceDirectivity', ''); if ~isempty(readstring) readstring = textscan(readstring, '%s', 'Delimiter', ','); readstring = readstring{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht if nargin<2 directivityName=readstring; else directivityName=readstring(srcID); end else obj.sourceDirectivity = {}; end % % Original % if iscell(obj.sourceDirectivity) % directivityName = obj.sourceDirectivity{sourceID + 1}; % else % obj.sourceDirectivity = textscan(obj.sourceDirectivityString, '%s', 'Delimiter', ','); % obj.sourceDirectivity = obj.sourceDirectivity{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht % directivityName = obj.sourceDirectivity{sourceID + 1}; % end end %------------------------------------------------------------------ function setSourceDirectivity(obj, directivity) if iscell(directivity) obj.sourceDirectivityString = obj.cat_cell_of_strings(directivity); obj.sourceDirectivity = directivity; else obj.sourceDirectivityString = directivity; if ~isempty(directivity) obj.sourceDirectivity = textscan(obj.sourceDirectivityString, '%s', 'Delimiter', ','); obj.sourceDirectivity = obj.sourceDirectivity{1}; % textscan liefert cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht end end obj.rpf_ini.SetValues('PrimarySources', 'sourceDirectivity', obj.sourceDirectivityString); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSourceDirectivityPath(obj, myPath) % % setSourceDirectivityPath(path) % % Set the RAVEN source directivity database path. Here all % directivities (.daff files) should be located % % To generate new files, check out opendaff in the % \applications\VirtualAcoustics\openDAFF\ % % RAVEN currently only supports OpenDAFFv1.5 files % obj.pathDirectivities = myPath; obj.rpf_ini.SetValues('Global', 'ProjectPath_DirectivityDB', myPath); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function pos = getSourcePosition(obj, srcID) all_pos = obj.rpf_ini.GetValues('PrimarySources', 'sourcePositions'); if (nargin < 2) pos = reshape(all_pos',3,length(all_pos)/3)'; else pos = all_pos((srcID*3 + 1) : (srcID*3 + 3)); end end %------------------------------------------------------------------ function setSourcePositions(obj, pos, coord_trafo) %setSourcePositions(pos, [coord_trafo]) % coord_trafo: 3-element vector % 1. element: which of the given axis to use as x-coordinates % 2. element: which of the given axis to use as y-coordinates % 3. element: which of the given axis to use as z-coordinates % if input is string convert to matrix if ischar(pos) pos = reshape(sscanf(pos, '%f,'), 3, numel(pos)/3)'; end % if input is itaCoordinates convert to matrix if isa(pos, 'itaCoordinates') pos = pos.cart; end % check dimensions if size(pos, 2) ~= 3 error('Positions must be given in vector (1x3) or matrix (Nx3) format.'); end % store matrix in RavenProject object obj.sourcePositions = pos; % store positions in .rpf file on disk if nargin > 2 obj.rpf_ini.SetValues('PrimarySources', 'sourcePositions', obj.make_proper_string( ... [sign(coord_trafo(1)) .* pos(:, abs(coord_trafo(1))), ... sign(coord_trafo(2)) .* pos(:, abs(coord_trafo(2))), ... sign(coord_trafo(3)) .* pos(:, abs(coord_trafo(3)))] )); else obj.rpf_ini.SetValues('PrimarySources', 'sourcePositions', obj.make_proper_string(pos)); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSourceViewVectors(obj, view) obj.sourceViewVectors = view; if isnumeric(view) || iscell(view) obj.rpf_ini.SetValues('PrimarySources', 'sourceViewVectors', obj.make_proper_string(view)); else obj.rpf_ini.SetValues('PrimarySources', 'sourceViewVectors', view); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSourceUpVectors(obj, up) obj.sourceUpVectors = up; if isnumeric(up) || iscell(up) obj.rpf_ini.SetValues('PrimarySources', 'sourceUpVectors', obj.make_proper_string(up)); else obj.rpf_ini.SetValues('PrimarySources', 'sourceUpVectors', up); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSourceSoundStates(obj, state) obj.sourceSoundStates = state; if isnumeric(state) || iscell(state) obj.rpf_ini.SetValues('PrimarySources', 'sourceSoundStates', obj.make_proper_string(state)); else obj.rpf_ini.SetValues('PrimarySources', 'sourceSoundStates', state); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSourceLevels(obj, levels) obj.sourceSoundLevels = levels; if isnumeric(levels) || iscell(levels) obj.rpf_ini.SetValues('PrimarySources', 'sourceSoundLevels', obj.make_proper_string(levels)); else obj.rpf_ini.SetValues('PrimarySources', 'sourceSoundLevels', levels); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end % [Receiver] % %------------------------------------------------------------------ function setReceiverNames(obj, rec_names) if iscell(rec_names) obj.receiverNameString = obj.cat_cell_of_strings(rec_names); obj.receiverNames = rec_names; else obj.receiverNameString = rec_names; obj.receiverNames = textscan(obj.receiverNameString, '%s', 'Delimiter', ','); obj.receiverNames = obj.receiverNames{1}; % textscan liefer cell array in nochmal einer zelle, diese doppelkapselung wird hier rückgängig gemacht end obj.rpf_ini.SetValues('Receiver', 'receiverNames', obj.receiverNameString); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function pos = getReceiverPosition(obj, recID) all_pos = obj.rpf_ini.GetValues('Receiver', 'receiverPositions'); if (nargin < 2) pos = reshape(all_pos',3,length(all_pos)/3)'; else pos = all_pos((recID*3 + 1) : (recID*3 + 3)); end end %------------------------------------------------------------------ function rvv = getReceiverViewVectors(obj, recID) all_pos = obj.rpf_ini.GetValues('Receiver', 'receiverViewVectors'); if (nargin < 2) rvv=reshape(all_pos',3,length(all_pos)/3)'; else rvv = all_pos((recID*3 + 1) : (recID*3 + 3)); end end %------------------------------------------------------------------ function ruv = getReceiverUpVectors(obj, recID) all_pos = obj.rpf_ini.GetValues('Receiver', 'receiverUpVectors'); if (nargin < 2) ruv = reshape(all_pos',3,length(all_pos)/3)'; else ruv = all_pos((recID*3 + 1) : (recID*3 + 3)); end end %------------------------------------------------------------------ function rss = getReceiverStates(obj,recID) rss = obj.rpf_ini.GetValues('Receiver', 'receiverStates', -1); if nargin==2 rss=rss(recID); end end %------------------------------------------------------------------ function rn = getReceiverNames(obj,recID) rn = obj.rpf_ini.GetValues('Receiver', 'receiverNames', -1); rn = textscan(rn, '%s', 'Delimiter', ','); rn = rn{1}; if nargin==2 rn=rn{1}{recID}; end end %------------------------------------------------------------------ function setReceiverPositions(obj, pos, coord_trafo) %setReceiverPositions(pos, coord_trafo) % coord_trafo: 3-element vector % 1. element: which of the given axis to use as x-coordinates % 2. element: which of the given axis to use as y-coordinates % 3. element: which of the given axis to use as z-coordinates % if input is string convert to matrix if ischar(pos) pos = reshape(sscanf(pos, '%f,'), 3, numel(pos)/3)'; end % if input is itaCoordinates convert to matrix if isa(pos, 'itaCoordinates') pos = pos.cart; end % check dimensions if size(pos, 2) ~= 3 error('Positions must be given in vector (1x3) or matrix (Nx3) format.'); end % store matrix in RavenProject object obj.receiverPositions = pos; % store positions in .rpf file on disk if nargin > 2 obj.rpf_ini.SetValues('Receiver', 'receiverPositions', obj.make_proper_string( ... [sign(coord_trafo(1)) .* pos(:, abs(coord_trafo(1))), ... sign(coord_trafo(2)) .* pos(:, abs(coord_trafo(2))), ... sign(coord_trafo(3)) .* pos(:, abs(coord_trafo(3)))] )); else obj.rpf_ini.SetValues('Receiver', 'receiverPositions', obj.make_proper_string(pos)); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setReceiverViewVectors(obj, view) obj.receiverViewVectors = view; if isnumeric(view) || iscell(view) obj.rpf_ini.SetValues('Receiver', 'receiverViewVectors', obj.make_proper_string(view)); else obj.rpf_ini.SetValues('Receiver', 'receiverViewVectors', view); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setReceiverUpVectors(obj, up) obj.receiverUpVectors = up; if isnumeric(up) || iscell(up) obj.rpf_ini.SetValues('Receiver', 'receiverUpVectors', obj.make_proper_string(up)); else obj.rpf_ini.SetValues('Receiver', 'receiverUpVectors', up); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setReceiverStates(obj, states) obj.receiverStates = states; if isnumeric(states) || iscell(states) obj.rpf_ini.SetValues('Receiver', 'receiverStates', obj.make_proper_string(states)); else obj.rpf_ini.SetValues('Receiver', 'receiverStates', states); end obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function directivityName = getReceiverHRTF(obj) directivityName = obj.rpf_ini.GetValues('Global', 'ProjectPath_HRTFDB'); end %------------------------------------------------------------------ function setReceiverHRTF(obj, HRTF) [pathHRTF,~,~] = fileparts(obj.fileHRTF); newFile = fullfile(pathHRTF, HRTF); % check if file path is absolute or relative path % (if ':" is in path, path is absolute) if ~isempty(strfind(HRTF, ':')) newFile = HRTF; else newFile = fullfile(pathHRTF, HRTF); end if exist(newFile, 'file') obj.fileHRTF = newFile; obj.rpf_ini.SetValues('Global', 'ProjectPath_HRTFDB', obj.fileHRTF); obj.rpf_ini.WriteFile(obj.ravenProjectFile); else error(['File not found: ' newFile]); end end % [ImageSources] % %------------------------------------------------------------------ function setISOrder_PS(obj, is_order) % setISOrder_PS % sets image source order % -1 = disabled % 0 = direct sound only % N = for N>1 : image source calculation up to order N % Do not set above 4, especially for complex models obj.ISOrder_PS = is_order; obj.rpf_ini.SetValues('ImageSources', 'ISOrder_PS', is_order); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setISOrder_SS(obj, is_order) obj.ISOrder_SS = is_order; obj.rpf_ini.SetValues('ImageSources', 'ISOrder_SS', is_order); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSkipDirectSound(obj, skipit) obj.ISSkipDirectSound = skipit; obj.rpf_ini.SetValues('ImageSources', 'ISSkipDirectSound', skipit); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end % [RayTracing] % %------------------------------------------------------------------ function setNumParticles(obj, num_particles) obj.numParticles_Sphere = num_particles; obj.rpf_ini.SetValues('RayTracing', 'numberOfParticles_DetectionSphere', num_particles); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setNumParticles_Portal(obj, num_particles) obj.numParticles_Portal = num_particles; obj.rpf_ini.SetValues('RayTracing', 'numberOfParticles_Portal', num_particles); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setFixReflectionPattern(obj, fixit) obj.fixReflectionPattern = fixit; obj.rpf_ini.SetValues('RayTracing', 'fixReflectionPattern', fixit); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setScatterModel(obj, diffuseOrHit) %setFilterResolution(obj, diffuseOrHit) % 0 = Diffuse Rain (or 'diffuse'), 1 = Hit Oriented (or 'hit') if ~isnumeric(diffuseOrHit) diffuseOrHit = double(strcmp(diffuseOrHit, 'hit')); end %obj.scatterModel = diffuseOrHit; obj.rpf_ini.SetValues('RayTracing', 'scatterModel_DetectionSphere', diffuseOrHit); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setScatterModel_Portal(obj, diffuseOrHit) %setFilterResolution(obj, diffuseOrHit) % 0 = Diffuse Rain (or 'diffuse'), 1 = Hit Oriented (or 'hit') if ~isnumeric(diffuseOrHit) diffuseOrHit = double(strcmp(diffuseOrHit, 'hit')); end %obj.scatterModel_Portal = diffuseOrHit; obj.rpf_ini.SetValues('RayTracing', 'scatterModel_Portal', diffuseOrHit); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setEnergyLoss(obj, energyloss) % setEnergyLoss % Define maximum energy loss for ray tracing (in dB) % typical value 60 dB % (for room acoustic parameter evaluation) % Increase to ~70-80 dB if high quality auralizations are desired % significantly increases simulation time (ray tracing) obj.energyLoss_Sphere = energyloss; obj.rpf_ini.SetValues('RayTracing', 'energyLoss_DetectionSphere', energyloss); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setEnergyLoss_Portal(obj, energyloss) obj.energyLoss_Portal = energyloss; obj.rpf_ini.SetValues('RayTracing', 'energyLoss_Portal', energyloss); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setTimeSlotLength(obj, slotlength) % setTimeSlotLength % set timeslot lengt in ms % (histogram resolution, for raytracing algorithm) obj.timeSlotLength = slotlength; obj.rpf_ini.SetValues('RayTracing', 'timeResolution_DetectionSphere', slotlength); obj.rpf_ini.SetValues('RayTracing', 'timeResolution_Portal', slotlength); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setRadiusDetectionSphere(obj, radius) obj.radiusSphere = radius; obj.rpf_ini.SetValues('RayTracing', 'radius_DetectionSphere', radius); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setDetectionSphereResolutionAzimuth(obj, resAzi) % setDetectionSphereResolutionAzimuth % set detection sphere azimuth resolution in degree % default value = 10 % samples sphere into grid, creating "directivity groups" % see Phd thesis by Dirk Schröder for details obj.detectionSphereAziResolution = resAzi; obj.rpf_ini.SetValues('RayTracing', 'resolutionAzimuth_DetectionSphere', resAzi); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setDetectionSphereResolutionElevation(obj, resEle) % setDetectionSphereResolutionElevation % set detection sphere elevation resolution in degree % default value = 10 % samples sphere into grid, creating "directivity groups" % see Phd thesis by Dirk Schröder for details obj.detectionSphereEleResolution = resEle; obj.rpf_ini.SetValues('RayTracing', 'resolutionElevation_DetectionSphere', resEle); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setFilterLength(obj, filter_length) % set filter length in ms if filter_length < 3 filter_length = filter_length * 1000; filter_length = fix(filter_length); disp(['Filter length should be set in milliseconds. Filter Length now set to ' num2str(filter_length) ' ms.']); else filter_length = fix(filter_length); end obj.filterLength = filter_length; obj.rpf_ini.SetValues('RayTracing', 'filterLength_DetectionSphere', filter_length); obj.rpf_ini.SetValues('RayTracing', 'filterLength_Portal', filter_length); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setFilterLengthToReverbTime(obj, roomID) % get approximate reverberation time and set filter length accordingly if iscell(obj.modelFileList) if (nargin < 2) disp('Please define which room should be used. Now using the first room (roomID = 0).'); roomID = 0; end if isempty(obj.model) roommodel = itaAc3dModel(obj.modelFileList{roomID + 1}); else roommodel = obj.model{roomID + 1}; end else if isempty(obj.model) roommodel = itaAc3dModel(obj.modelFileList); else roommodel = obj.model; end end RT = roommodel.getReverbTime(obj.pathMaterials); filter_length = max(RT) * 1.2 * 1000; % filter length parameter needed in [ms] and 20% additional headroom obj.setFilterLength(filter_length); end % [Filter] % %------------------------------------------------------------------ function setSamplingFrequency(obj, samplingRate) obj.sampleRate = samplingRate; obj.rpf_ini.SetValues('Filter', 'samplingFrequency', samplingRate); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setFixPoissonSequence(obj, fixit) obj.fixPoissonSequence = fixit; obj.rpf_ini.SetValues('Filter', 'setFixPoissonSequence', fixit); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setPoissonSequenceNumber(obj, number) obj.poissonSequenceNumber = number; obj.rpf_ini.SetValues('Filter', 'poissonSequenceNumber', number); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setFilterResolution(obj, octOr3rd) %setFilterResolution(obj, octOr3rd) % 1 = Octave (or 'oct'), 0 = 3rd Octave (or '3rd') if ~isnumeric(octOr3rd) octOr3rd = double(strcmp(octOr3rd, 'oct')); if ~octOr3rd octOr3rd = double(strcmp(octOr3rd, '3rd')); end end obj.filterResolution = octOr3rd; obj.rpf_ini.SetValues('Filter', 'filterResolution', octOr3rd); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setMaximumReflectionDensity(obj, maxReflectionDensity) % setMaximumReflectionDensity % sets the maxium reflection density for filter synthesis. Default % value is 20000. Typical values between 10000 and 20000 % for details see Diss Dirk Schröder or Aspöck @ DAGA 2013 / DAGA 2017 obj.maxReflectionDensity = maxReflectionDensity; obj.rpf_ini.SetValues('Filter', 'maximumReflectionDensity', maxReflectionDensity); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setAmbisonicsOrder(obj, order) obj.ambisonicsOrder = order; obj.rpf_ini.SetValues('Filter', 'ambisonicsOrder', order); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setNumberSpreadedSources(obj, numberSpreadedSrc) obj.numberSpreadedSources = numberSpreadedSrc; obj.rpf_ini.SetValues('Filter', 'numberSpreadedSources', numberSpreadedSrc); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setSpreadingStdDeviation(obj, spreadingStdDev) obj.spreadingStdDeviation = spreadingStdDev; obj.rpf_ini.SetValues('Filter', 'spreadingStdDeviation', spreadingStdDev); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setFftDegreeForWallFilterInterpolation(obj, fft_degree) obj.fftDegreeForWallFilterInterpolation = fft_degree; obj.rpf_ini.SetValues('Filter', 'fftDegreeForWallFilterInterpolation', fft_degree); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setOutputPath(obj, path) % % setOutputPath(path) % % Set the RAVEN output path. If you keep the output files % (cf. keepOutputFiles variable / method keepOutputFiles), % Impulse Responses as wave files are stored here : % ///) % % Absorption and scattering plots are also saved here. % % Use method openOutputFolder to open result folder. % obj.pathResults = path; obj.rpf_ini.SetValues('Global', 'ProjectPath_Output', path); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setMaterialDatabasePath(obj, path) % % setMaterialDatabasePath(path) % % Set the RAVEN material database path. Here all materials % (.mat files) should be located, which are contained in % your model (see rpf.getRoomMaterialNames or directly in % your model file .ac % obj.pathMaterials = path; obj.rpf_ini.SetValues('Global', 'ProjectPath_MaterialDB', path); obj.rpf_ini.WriteFile(obj.ravenProjectFile); end %------------------------------------------------------------------ function setMaterial(obj, materialname, alpha, scatter, roomID) % % setMaterial(materialname, alpha, scatter, roomID) % % materialname can be the filename of the .mat file (without % extension or a number to address the material ID if the % current room (can be further specified by roomID). % % alpha and scatter have to be size of % - 31 elements (for one-third octave resolution) % - 10 elements (for octave resolution) % - or a single number (flat spectrum) if nargin < 4 scatter = 0.3; end % check if a material ID was given instead of a material name if isnumeric(materialname) % addressing by material ID (which is dependent on the room that is currently loaded!) matID = materialname; if nargin < 5 roomID = 0; end materialsInCurrentRoom = obj.getRoomMaterialNames(roomID); matname = materialsInCurrentRoom{matID + 1}; else matname = materialname; end mat_filename = fullfile(obj.pathMaterials, [matname '.mat']); % open file for writing fid = fopen(mat_filename, 'wt'); % print header fprintf(fid, '[Material]\r\n'); fprintf(fid, ['name=' matname '\r\n']); fprintf(fid, ['notes=Generated/Modified by MATLAB script (' datestr(now) ')\r\n']); % parse input data if length(alpha) == 1 alpha = ones(1,31) * alpha; elseif length(alpha) ~= 31 && length(alpha) ~= 10 disp('ERROR! alpha has to be a vector with the length 1,10, or 31.'); return; end if length(scatter) == 1 scatter = ones(1,31) * scatter; elseif length(scatter) ~= 31 && length(scatter) ~= 10 disp('ERROR! scatter has to be a vector with the length 1,10, or 31.'); return; end % print absorption values if length(alpha)==31 fprintf(fid, ['absorp=' num2str(alpha(1), '%2.4f')]); for i = 1 : 30 fprintf(fid, [', ' num2str(alpha(i+1), '%2.4f')]); end else fprintf(fid, ['absorp=' num2str(alpha(1), '%2.4f')]); for i = 1 : 30 fprintf(fid, [', ' num2str(alpha(int32(i+1)/3), '%2.4f')]); end end fprintf(fid, '\r\n'); % print scattering values if length(scatter)==31 fprintf(fid, ['scatter=' num2str(scatter(1), '%2.4f')]); for i = 1 : 30 fprintf(fid, [', ' num2str(scatter(i+1), '%2.4f')]); end else fprintf(fid, ['scatter=' num2str(scatter(1), '%2.4f')]); for i = 1 : 30 fprintf(fid, [', ' num2str(scatter(int32(i+1)/3), '%2.4f')]); end end fprintf(fid, '\r\n'); % print interpolation values if length(alpha)==31 fprintf(fid, 'interpol=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\r\n'); else fprintf(fid, 'interpol=1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1\r\n'); end % close file fclose(fid); end %------------------------------------------------------------------ function order = getAmbisonicsOrder(obj) order = obj.rpf_ini.GetValues('Filter', 'ambisonicsOrder', -1); end %------------------------------------------------------------------ function numberSpreadedSources = getNumberSpreadedSources(obj) numberSpreadedSources = obj.rpf_ini.GetValues('Filter', 'numberSpreadedSources', 0); end %------------------------------------------------------------------ function spreadingStdDeviation = getSpreadingStdDeviation(obj) spreadingStdDeviation = obj.rpf_ini.GetValues('Filter', 'spreadingStdDeviation', 0.2); end %------------------------------------------------------------------ function sourceviewvectors = getSourceViewVectors(obj,srcID) all_pos = obj.rpf_ini.GetValues('PrimarySources', 'sourceViewVectors'); if (nargin < 2) sourceviewvectors=reshape(all_pos',3,length(all_pos)/3)'; else sourceviewvectors = all_pos((srcID*3 + 1) : (srcID*3 + 3)); end end %------------------------------------------------------------------ function suv = getSourceUpVectors(obj,srcID) all_pos = obj.rpf_ini.GetValues('PrimarySources', 'sourceUpVectors'); if (nargin < 2) suv = reshape(all_pos',3,length(all_pos)/3)'; else suv = all_pos((srcID*3 + 1) : (srcID*3 + 3)); end end %------------------------------------------------------------------ function sss = getSourceSoundStates(obj,srcID) sss = obj.rpf_ini.GetValues('PrimarySources', 'sourceSoundStates', -1); if nargin==2 sss=sss(srcID); end end %------------------------------------------------------------------ function ssl = getSourceSoundLevels(obj,srcID) ssl = obj.rpf_ini.GetValues('PrimarySources', 'sourceSoundLevels', -1); if nargin==2 ssl=ssl(srcID); end end %------------------------------------------------------------------ function sn = getSourceNames(obj,srcID) sn = obj.rpf_ini.GetValues('PrimarySources', 'sourceNames', -1); sn = textscan(sn, '%s', 'Delimiter', ','); sn = sn{1}; if nargin==2 sn=sn{1}{srcID}; end end %------------------------------------------------------------------ function [absorp, scatter] = getMaterial(obj, materialName, roomID) % check if a material ID was given instead of a material name if isnumeric(materialName) % addressing by material ID (which is dependent on the room that is currently loaded!) matID = materialName; if nargin < 3 roomID = 0; end materialsInCurrentRoom = obj.getRoomMaterialNames(roomID); matname = materialsInCurrentRoom{matID + 1}; else if (numel(materialName) > 4) && isequal(materialName(end-3:end), '.mat') materialName = materialName(1:end-4); end matname = materialName; end ini = IniConfig(); ini.ReadFile(fullfile(obj.pathMaterials, [matname '.mat'])); absorp = ini.GetValues('Material', 'absorp'); if nargout > 1 scatter = ini.GetValues('Material', 'scatter'); end end %------------------------------------------------------------------ function monoIR = getImpulseResponse(obj) if isempty(obj.monoIR) error('No mono room impulse response present.'); end if numel(obj.monoIR) > 1 monoIR = obj.monoIR; else monoIR = obj.monoIR{1,1}; end end %------------------------------------------------------------------ function monoIRitaAudio = getImpulseResponseItaAudio(obj) if isempty(obj.monoIR) error('No mono room impulse response present.'); end monoIRitaAudio = itaAudio(size(obj.monoIR)); for iSrc = 1 : size(obj.monoIR, 1) for iRec = 1 : size(obj.monoIR, 2) monoIRitaAudio(iSrc,iRec).timeData = obj.monoIR{iSrc,iRec}; monoIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; monoIRitaAudio(iSrc,iRec).signalType = 'energy'; monoIRitaAudio(iSrc,iRec).channelNames{1} = obj.receiverNames{iRec}; %Make sure number of samples is even if mod(monoIRitaAudio(iSrc,iRec).nSamples,2)==1 monoIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function monoIR_IS = getImpulseResponseImageSources(obj) if isempty(obj.monoIR_IS) error('No mono image source impulse response present.'); end if numel(obj.monoIR_IS) > 1 monoIR_IS = obj.monoIR_IS; else monoIR_IS = obj.monoIR_IS{1,1}; end end %------------------------------------------------------------------ function monoIR_RT = getImpulseResponseRayTracing(obj) if isempty(obj.monoIR_RT) error('No mono ray tracing impulse response present.'); end if numel(obj.monoIR_RT) > 1 monoIR_RT = obj.monoIR_RT; else monoIR_RT = obj.monoIR_RT{1,1}; end end %------------------------------------------------------------------ function monoIRitaAudio = getImpulseResponseImageSourcesItaAudio(obj) if isempty(obj.monoIR_IS) error('No mono image source impulse response present.'); end monoIRitaAudio = itaAudio(size(obj.monoIR_IS)); for iSrc = 1 : size(obj.monoIR_IS, 1) for iRec = 1 : size(obj.monoIR_IS, 2) monoIRitaAudio(iSrc,iRec).timeData = obj.monoIR_IS{iSrc,iRec}; monoIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; monoIRitaAudio(iSrc,iRec).signalType = 'energy'; monoIRitaAudio(iSrc,iRec).channelNames{1} = obj.receiverNames{iRec}; %Make sure number of samples is even if mod(monoIRitaAudio(iSrc,iRec).nSamples,2)==1 monoIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function monoIRitaAudio = getImpulseResponseRayTracingItaAudio(obj) if isempty(obj.monoIR_RT) error('No mono ray tracing impulse response present.'); end monoIRitaAudio = itaAudio(size(obj.monoIR_RT)); for iSrc = 1 : size(obj.monoIR_RT, 1) for iRec = 1 : size(obj.monoIR_RT, 2) monoIRitaAudio(iSrc,iRec).timeData = obj.monoIR_RT{iSrc,iRec}; monoIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; monoIRitaAudio(iSrc,iRec).signalType = 'energy'; monoIRitaAudio(iSrc,iRec).channelNames{1} = obj.receiverNames{iRec}; %Make sure number of samples is even if mod(monoIRitaAudio(iSrc,iRec).nSamples,2)==1 monoIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function binauralIR = getBinauralImpulseResponse(obj) if isempty(obj.binauralIR) error('No binaural impulse response present.'); end if numel(obj.binauralIR) > 1 binauralIR = obj.binauralIR; else binauralIR = obj.binauralIR{1,1}; end end %------------------------------------------------------------------ function binauralIRitaAudio = getBinauralImpulseResponseItaAudio(obj) if isempty(obj.binauralIR) error('No binaural impulse response present.'); end binauralIRitaAudio = itaAudio(size(obj.binauralIR)); for iSrc = 1 : size(obj.binauralIR, 1) for iRec = 1 : size(obj.binauralIR, 2) binauralIRitaAudio(iSrc,iRec).timeData = obj.binauralIR{iSrc,iRec}; binauralIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; binauralIRitaAudio(iSrc,iRec).signalType = 'energy'; if isempty(binauralIRitaAudio(iSrc,iRec)); continue; end binauralIRitaAudio(iSrc,iRec).channelNames{1} = [obj.receiverNames{iRec} '_Left']; binauralIRitaAudio(iSrc,iRec).channelNames{2} = [obj.receiverNames{iRec} '_Right']; %Make sure number of samples is even if mod(binauralIRitaAudio(iSrc,iRec).nSamples,2)==1 binauralIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function binauralIR_IS = getBinauralImpulseResponseImageSources(obj) if isempty(obj.binauralIR_IS) error('No binaural image source impulse response present.'); end if numel(obj.binauralIR_IS) > 1 binauralIR_IS = obj.binauralIR_IS; else binauralIR_IS = obj.binauralIR_IS{1,1}; end end %------------------------------------------------------------------ function binauralIRitaAudio = getBinauralImpulseResponseImageSourcesItaAudio(obj) if isempty(obj.binauralIR_IS) error('No binaural image source impulse response present.'); end binauralIRitaAudio = itaAudio(size(obj.binauralIR_IS)); for iSrc = 1 : size(obj.binauralIR_IS, 1) for iRec = 1 : size(obj.binauralIR_IS, 2) binauralIRitaAudio(iSrc,iRec).timeData = obj.binauralIR_IS{iSrc,iRec}; binauralIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; binauralIRitaAudio(iSrc,iRec).signalType = 'energy'; if isempty(binauralIRitaAudio(iSrc,iRec)); continue; end binauralIRitaAudio(iSrc,iRec).channelNames{1} = [obj.receiverNames{iRec} '_Left']; binauralIRitaAudio(iSrc,iRec).channelNames{2} = [obj.receiverNames{iRec} '_Right']; %Make sure number of samples is even if mod(binauralIRitaAudio(iSrc,iRec).nSamples,2)==1 binauralIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function binauralIR_RT = getBinauralImpulseResponseRayTracing(obj) if isempty(obj.binauralIR_RT) error('No binaural ray tracing impulse response present.'); end if numel(obj.binauralIR_RT) > 1 binauralIR_RT = obj.binauralIR_RT; else binauralIR_RT = obj.binauralIR_RT{1,1}; end end %------------------------------------------------------------------ function binauralIRitaAudio = getBinauralImpulseResponseRayTracingItaAudio(obj) if isempty(obj.binauralIR_RT) error('No binaural image source impulse response present.'); end binauralIRitaAudio = itaAudio(size(obj.binauralIR_RT)); for iSrc = 1 : size(obj.binauralIR_RT, 1) for iRec = 1 : size(obj.binauralIR_RT, 2) binauralIRitaAudio(iSrc,iRec).timeData = obj.binauralIR_RT{iSrc,iRec}; binauralIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; binauralIRitaAudio(iSrc,iRec).signalType = 'energy'; if isempty(binauralIRitaAudio(iSrc,iRec)); continue; end binauralIRitaAudio(iSrc,iRec).channelNames{1} = [obj.receiverNames{iRec} '_Left']; binauralIRitaAudio(iSrc,iRec).channelNames{2} = [obj.receiverNames{iRec} '_Right']; %Make sure number of samples is even if mod(binauralIRitaAudio(iSrc,iRec).nSamples,2)==1 binauralIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function binauralPoissonSequenceitaAudio = getBinauralPoissonSequenceItaAudio(obj) data = load([ obj.pathResults '\binauralPoissonNoiseProcess.txt']); binauralPoissonSequenceitaAudio = itaAudio; binauralPoissonSequenceitaAudio.samplingRate = obj.sampleRate; binauralPoissonSequenceitaAudio.signalType = 'energy'; binauralPoissonSequenceitaAudio.timeData = data; if isempty(binauralPoissonSequenceitaAudio); return; end binauralPoissonSequenceitaAudio.channelNames{1} = [obj.receiverNames{end} '_Left']; binauralPoissonSequenceitaAudio.channelNames{2} = [obj.receiverNames{end} '_Right']; %Make sure number of samples is even if mod(binauralPoissonSequenceitaAudio.nSamples,2)==1 binauralPoissonSequenceitaAudio.timeData(end, :) = []; end end %------------------------------------------------------------------ function binauralImpulseResponse = generateBinauralImpulseResponseItaAudio(obj) % get noise process (containing poisson distributed HRIRs without temporal or spectral weighting) poisson = obj.getBinauralPoissonSequenceItaAudio(); % add white noise to poisson sequence (for more realistic simulaitons % noiseLevel = 10^(-20/10) * max(max(poisson.timeData)); % addNoise = ita_generate('whitenoise', noiseLevel, poisson.samplingRate, poisson.fftDegree); % addNoise.timeData(:,2) = addNoise.timeData(:,1); % add second channel % poisson = poisson + addNoise; % get ray tracing results % histo = obj.getHistogramRayTracing(); histo = obj.histogramRT(); if iscell(histo) if numel(histo) > 1 disp('WARNING! Generating binaural impulse response only for the last receiver!') end histo = histo{end}; end % low pass late part of the histogram to avoid bubble sound histo.data = obj.applyLowPassToLatePartOfHistogram(histo.data); % convert added energy histogram to sound pressure envelope histogram and % time-interpolate time slot quantized histogram to full number of samples length histoInterpolatedToSamples = interp1(histo.timevector, sqrt(histo.data), poisson.timeVector, 'cubic'); histoInterpolatedToSamples(histoInterpolatedToSamples < 0) = 0; % energyScalingFactor = sum(histo.data) ./ sum(histoInterpolatedToSamples); % is vector of frequencies, sums up the time slots % energyScalingFactor(isnan(energyScalingFactor)) = 0; % histoInterpolatedToSamples = histoInterpolatedToSamples .* repmat(energyScalingFactor, poisson.nSamples, 1); % band-pass poisson sequence if size(histo.data, 2) == 10 poissonFilteredLeft = ita_mpb_filter(poisson.ch(1), 'oct', 1, 'order', 6); poissonFilteredRight = ita_mpb_filter(poisson.ch(2), 'oct', 1, 'order', 6); else poissonFilteredLeft = ita_mpb_filter(poisson.ch(1), 'oct', 3, 'order', 6); poissonFilteredRight = ita_mpb_filter(poisson.ch(2), 'oct', 3, 'order', 6); % toolbox 30 thirdocts vs. ravens 31 thirdocts, delete 20 Hz data histoInterpolatedToSamples(:, 1) = []; end % temporal shaping of the poisson sequence using the sound presssure envelope from the histogram poissonFilteredLeft.timeData(:, :) = poissonFilteredLeft.timeData(:,:) .* histoInterpolatedToSamples; poissonFilteredRight.timeData(:, :) = poissonFilteredRight.timeData(:, :) .* histoInterpolatedToSamples; % create new itaAudio to store resulting filter binauralImpulseResponse = itaAudio; binauralImpulseResponse.samplingRate = obj.sampleRate; binauralImpulseResponse.signalType = 'energy'; binauralImpulseResponse.timeData(:,:) = [sum(poissonFilteredLeft.timeData, 2) sum(poissonFilteredRight.timeData, 2)]; % total energy scaling (total energy histogram == total energy impulse response of one channel) % convert averaged energy histogram to added energy histogram histo.data = obj.powerSpectrum2energySpectrum(histo.data); totalEnergyHistogram = sum(sum(histo.data,2)); numTimeSlots = size(histo.data, 1); maxEnergyHistogram = (1 * numTimeSlots); % 1 == EnergyAdded for one timeslot if all frequencyies are = 1 maxEnergyPoisson = sum(poisson.timeData.^2); % vector(2) for left and right channel totalEnergyImpulseResponse = sum(binauralImpulseResponse.timeData.^2); % vector(2) for left and right channel binauralImpulseResponse.timeData = binauralImpulseResponse.timeData .* sqrt((maxEnergyPoisson / totalEnergyImpulseResponse) * (totalEnergyHistogram / maxEnergyHistogram)); % add image sources BRIR_IS = obj.getBinauralImpulseResponseImageSourcesItaAudio(); [BRIR_IS, binauralImpulseResponse] = ita_extend_dat(BRIR_IS, binauralImpulseResponse); binauralImpulseResponse = BRIR_IS + binauralImpulseResponse; end %------------------------------------------------------------------ function ambisonicsIR = getAmbisonicsImpulseResponse(obj) if isempty(obj.ambisonicsIR) error('No ambisonics b-format impulse response present.'); end if numel(obj.ambisonicsIR) > 1 ambisonicsIR = obj.ambisonicsIR; else ambisonicsIR = obj.ambisonicsIR{1,1}; end end %------------------------------------------------------------------ function ambisonicsIRita = getAmbisonicsImpulseResponseItaAudio(obj) if isempty(obj.ambisonicsIR) error('No ambisonics b-format impulse response present.'); end ambisonicsIRita = itaAudio(size(obj.ambisonicsIR)); for iSrc = 1 : size(obj.ambisonicsIR, 1) for iRec = 1 : size(obj.ambisonicsIR, 2) ambisonicsIRita(iSrc,iRec).timeData = obj.ambisonicsIR{iSrc,iRec}; ambisonicsIRita(iSrc,iRec).samplingRate = obj.sampleRate; ambisonicsIRita(iSrc,iRec).signalType = 'energy'; for linInd = 1 : ambisonicsIRita(iSrc,iRec).nChannels [deg, ord] = ita_sph_linear2degreeorder(linInd); ambisonicsIRita(iSrc,iRec).channelNames{linInd} = [obj.receiverNames{iRec} '_Degree(' num2str(deg) ')_Order(' num2str(ord) ')']; end %Make sure number of samples is even if mod(ambisonicsIRita(iSrc,iRec).nSamples,2)==1 ambisonicsIRita(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function ambisonicsIR_IS = getAmbisonicsImpulseResponseImageSources(obj) if isempty(obj.ambisonicsIR_IS) error('No image sources ambisonics b-format impulse response present.'); end if numel(obj.ambisonicsIR_IS) > 1 ambisonicsIR_IS = obj.ambisonicsIR_IS; else ambisonicsIR_IS = obj.ambisonicsIR_IS{1,1}; end end %------------------------------------------------------------------ function ambisonicsIRita = getAmbisonicsImpulseResponseImageSourcesItaAudio(obj) if isempty(obj.ambisonicsIR_IS) error('No image sources ambisonics b-format impulse response present.'); end ambisonicsIRita = itaAudio(size(obj.ambisonicsIR_IS)); for iSrc = 1 : size(obj.ambisonicsIR_IS, 1) for iRec = 1 : size(obj.ambisonicsIR_IS, 2) ambisonicsIRita(iSrc,iRec).timeData = obj.ambisonicsIR_IS{iSrc,iRec}; ambisonicsIRita(iSrc,iRec).samplingRate = obj.sampleRate; ambisonicsIRita(iSrc,iRec).signalType = 'energy'; for linInd = 1 : ambisonicsIRita(iSrc,iRec).nChannels [deg, ord] = ita_sph_linear2degreeorder(linInd); ambisonicsIRita(iSrc,iRec).channelNames{linInd} = [obj.receiverNames{iRec} '_Degree(' num2str(deg) ')_Order(' num2str(ord) ')']; end %Make sure number of samples is even if mod(ambisonicsIRita(iSrc,iRec).nSamples,2)==1 ambisonicsIRita(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function ambisonicsIR_RT = getAmbisonicsImpulseResponseRayTracing(obj) if isempty(obj.ambisonicsIR_RT) error('No ray tracing ambisonics b-format impulse response present.'); end if numel(obj.ambisonicsIR_RT) > 1 ambisonicsIR_RT = obj.ambisonicsIR_RT; else ambisonicsIR_RT = obj.ambisonicsIR_RT{1,1}; end end %------------------------------------------------------------------ function ambisonicsIRita = getAmbisonicsImpulseResponseRayTracingItaAudio(obj) if isempty(obj.ambisonicsIR_RT) error('No ray tracing ambisonics b-format impulse response present.'); end ambisonicsIRita = itaAudio(size(obj.ambisonicsIR_RT)); for iSrc = 1 : size(obj.ambisonicsIR_RT, 1) for iRec = 1 : size(obj.ambisonicsIR_RT, 2) ambisonicsIRita(iSrc,iRec).timeData = obj.ambisonicsIR_RT{iSrc,iRec}; ambisonicsIRita(iSrc,iRec).samplingRate = obj.sampleRate; ambisonicsIRita(iSrc,iRec).signalType = 'energy'; for linInd = 1 : ambisonicsIRita(iSrc,iRec).nChannels [deg, ord] = ita_sph_linear2degreeorder(linInd); ambisonicsIRita(iSrc,iRec).channelNames{linInd} = [obj.receiverNames{iRec} '_Degree(' num2str(deg) ')_Order(' num2str(ord) ')']; end %Make sure number of samples is even if mod(ambisonicsIRita(iSrc,iRec).nSamples,2)==1 ambisonicsIRita(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function vbapIR = getVBAPImpulseResponse(obj) if isempty(obj.vbapIR) error('No vbap impulse response present.'); end if numel(obj.vbapIR) > 1 vbapIR = obj.vbapIR; else vbapIR = obj.vbapIR{1,1}; end end %------------------------------------------------------------------ function vbapIRita = getVBAPImpulseResponseItaAudio(obj) if isempty(obj.vbapIR) error('No vbap impulse response present.'); end vbapIRita = itaAudio(size(obj.vbapIR)); for iSrc = 1 : size(obj.vbapIR, 1) for iRec = 1 : size(obj.vbapIR, 2) vbapIRita(iSrc,iRec).timeData = obj.vbapIR{iSrc,iRec}; vbapIRita(iSrc,iRec).samplingRate = obj.sampleRate; vbapIRita(iSrc,iRec).signalType = 'energy'; for linInd = 1 : vbapIRita(iSrc,iRec).nChannels vbapIRita(iSrc,iRec).channelNames{linInd} = [obj.receiverNames{iRec} '_Speaker_' num2str(linInd)]; end %Make sure number of samples is even if mod(vbapIRita(iSrc,iRec).nSamples,2)==1 vbapIRita(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function vbapIR_IS = getVBAPImpulseResponseImageSources(obj) if isempty(obj.vbapIR_IS) error('No image sources vbap impulse response present.'); end if numel(obj.vbapIR_IS) > 1 vbapIR_IS = obj.vbapIR_IS; else vbapIR_IS = obj.vbapIR_IS{1,1}; end end %------------------------------------------------------------------ function vbapIRita = getVBAPImpulseResponseImageSourcesItaAudio(obj) if isempty(obj.vbapIR_IS) error('No image sources vbap impulse response present.'); end vbapIRita = itaAudio(size(obj.vbapIR_IS)); for iSrc = 1 : size(obj.vbapIR_IS, 1) for iRec = 1 : size(obj.vbapIR_IS, 2) vbapIRita(iSrc,iRec).timeData = obj.vbapIR_IS{iSrc,iRec}; vbapIRita(iSrc,iRec).samplingRate = obj.sampleRate; vbapIRita(iSrc,iRec).signalType = 'energy'; for linInd = 1 : vbapIRita(iSrc,iRec).nChannels vbapIRita(iSrc,iRec).channelNames{linInd} = [obj.receiverNames{iRec} '_Speaker_' num2str(linInd)]; end %Make sure number of samples is even if mod(vbapIRita(iSrc,iRec).nSamples,2)==1 vbapIRita(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function vbapIR_RT = getVBAPImpulseResponseRayTracing(obj) if isempty(obj.vbapIR_RT) error('No ray tracing vbap impulse response present.'); end if numel(obj.vbapIR_RT) > 1 vbapIR_RT = obj.vbapIR_RT; else vbapIR_RT = obj.vbapIR_RT{1,1}; end end %------------------------------------------------------------------ function vbapIRita = getVBAPImpulseResponseRayTracingItaAudio(obj) if isempty(obj.vbapIR_RT) error('No ray tracing vbap impulse response present.'); end vbapIRita = itaAudio(size(obj.vbapIR_RT)); for iSrc = 1 : size(obj.vbapIR_RT, 1) for iRec = 1 : size(obj.vbapIR_RT, 2) vbapIRita(iSrc,iRec).timeData = obj.vbapIR_RT{iSrc,iRec}; vbapIRita(iSrc,iRec).samplingRate = obj.sampleRate; vbapIRita(iSrc,iRec).signalType = 'energy'; for linInd = 1 : vbapIRita(iSrc,iRec).nChannels vbapIRita(iSrc,iRec).channelNames{linInd} = [obj.receiverNames{iRec} '_Speaker_' num2str(linInd)]; end %Make sure number of samples is even if mod(vbapIRita(iSrc,iRec).nSamples,2)==1 vbapIRita(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function histo = getHistogram(obj) if isempty(obj.histogram) error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end if numel(obj.histogramRT) > 1 histo = obj.histogram; else histo = obj.histogram{1,1}; end end %------------------------------------------------------------------ function histo = getHistogram_itaResult(obj) histo_data = obj.getHistogram(); if ~iscell(histo_data) histo_data = {histo_data}; end % pre-alloc multi-instance for all sources and receivers histo(size(obj.histogram,1), size(obj.histogram,2)) = itaResult(); for iHisto = 1 : numel(obj.histogram) if ~isempty(obj.histogram{iHisto}) histo(iHisto).timeData = histo_data{iHisto}.data; histo(iHisto).timeVector = histo_data{iHisto}.timevector; if obj.filterResolution == 1 histo(iHisto).channelNames = strcat('Histogram ', obj.freqLabelOct); else histo(iHisto).channelNames = strcat('Histogram ', obj.freqLabel3rd); end end end end %------------------------------------------------------------------ function histo = getHistogramRayTracing(obj) if isempty(obj.histogramRT) raytracing_histo_files = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), '_RT', '.hst'); % read the histogram files back from disk obj.histogramRT = obj.readHistogram(raytracing_histo_files); end if numel(obj.histogramRT) > 1 histo = obj.histogramRT; else histo = obj.histogramRT{1,1}; end end %------------------------------------------------------------------ function histo = getHistogramEnergyAdded(obj) if isempty(obj.histogram) error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end numberFrequencyBands = size(obj.histogram{find(~cellfun(@isempty, obj.histogram), 1)}.data, 2); numberTimeSteps = size(obj.histogram{find(~cellfun(@isempty, obj.histogram), 1)}.data, 1); if numberFrequencyBands == 10 centerFreqs = [31.25, 62.50, 125.00, 250.00, 500.00, 1000.00, 2000.00, 4000.00, 8000.00,16000.00]; lowerBandLimits = centerFreqs / sqrt(2); upperBandLimits = centerFreqs * sqrt(2); bandWidths = upperBandLimits - lowerBandLimits; %bandWidths = [0.0010 0.0020 0.0041 0.0082 0.0163 0.0326 0.0652 0.1302 0.2591 0.4870]; %bandWidths = [44.2 44.2 88.4 176.8 353.5 707.1 1414.2 2828.5 5656.9 11313.8]; elseif numberFrequencyBands == 31 centerFreqs = [19.686, 24.803, 31.250, 39.373, 49.606, 62.500, 78.745, 99.213, 125.000, 157.490, 198.425, 250.000, 314.980, 396.850, 500.000, 629.961, 793.701, 1000.000, 1259.921, 1587.401, 2000.000, 2519.842, 3174.802, 4000.000, 5039.684, 6349.604, 8000.000,10079.368,12699.208,16000.000,20158.737]; lowerBandLimits = centerFreqs / sqrt(2); upperBandLimits = centerFreqs * sqrt(2); bandWidths = upperBandLimits - lowerBandLimits; %bandWidths = [0.0002095384,0.0002644467,0.0003336922,0.0004203726,0.0005296132,0.0006672770,0.0008407238,0.0010592415,0.0013345612,0.0016814420,0.0021184840,0.0026691217,0.0033628807,0.0042369610,0.0053382298,0.0067257351,0.0084738691,0.0106763538,0.0134512584,0.0169473147,0.0213518612,0.0269008258,0.0338912528,0.0426969859,0.0537882314,0.0677558315,0.0853411434,0.1074724174,0.1353085098,0.1702908522,0.1861151317]; %bandWidths = [22.4 5.5 7.2 9.1 11.4 14.5 18.3 23.0 28.9 36.5 46.0 57.9 72.9 91.8 115.8 145.9 183.8 231.6 291.7 367.6 463.1 583.5 735.2 926.3 1167.0 1470.4 1852.5 2334.0 2941.2 3705.0 4667.6]; else error('Unkown frequency resolution.'); end bandWidthTotal = sum(bandWidths); histo = obj.histogram; for iHisto = 1 : numel(histo) histo{iHisto}.data = histo{iHisto}.data(:,:) .* repmat(bandWidths./bandWidthTotal, numberTimeSteps, 1); % disp('ACHTUNG! Histogramm ist nicht von energy averaged auf energy added konvertiert worden!'); end if numel(histo) < 2 histo = histo{1,1}; end end %------------------------------------------------------------------ function histo = getHistogramRayTracingEnergyAdded(obj) if isempty(obj.histogramRT) error('No ray tracing histogram present.'); end numberFrequencyBands = size(obj.histogramRT{find(~cellfun(@isempty, obj.histogramRT), 1)}.data, 2); numberTimeSteps = size(obj.histogramRT{find(~cellfun(@isempty, obj.histogramRT), 1)}.data, 1); if numberFrequencyBands == 10 bandWidths = [44.2 44.2 88.4 176.8 353.5 707.1 1414.2 2828.5 5656.9 11313.8]; elseif numberFrequencyBands == 31 bandWidths = [22.4 5.5 7.2 9.1 11.4 14.5 18.3 23.0 28.9 36.5 46.0 57.9 72.9 91.8 115.8 145.9 183.8 231.6 291.7 367.6 463.1 583.5 735.2 926.3 1167.0 1470.4 1852.5 2334.0 2941.2 3705.0 4667.6]; else error('Unkown frequency resolution.'); end bandWidthTotal = sum(bandWidths); histo = obj.histogramRT; for iHisto = 1 : numel(histo) histo{iHisto}.data = histo{iHisto}.data .* repmat(bandWidths./bandWidthTotal, numberTimeSteps, 1); % disp('ACHTUNG! Histogramm ist nicht von energy averaged auf energy added konvertiert worden!'); end if numel(histo) < 2 histo = histo{1,1}; end end %------------------------------------------------------------------ function wall_log = getWallHitLog_IS(obj, logFile) if (nargin < 2) && isempty(obj.wallHitLog_IS) error('No wall hit log present.'); end if nargin > 1 obj.wallHitLog_IS = obj.readWallHitLog_IS(logFile); end if numel(obj.wallHitLog_IS) > 1 wall_log = obj.wallHitLog_IS; else wall_log = obj.wallHitLog_IS{1,1}; end end %------------------------------------------------------------------ function planewaves = getPlaneWaves_IS(obj, sourceID, receiverID) % % planewaves = getPlaneWaves_IS(sourceID, receiverID) % % Returns the plane waves of the image source simulation. % The returned freqData values are linear absolute sound pressure. % sourceID : optional ID of the sound source % receiverID : optional ID of the receiver % (counting starts with zero) if isempty(obj.planeWaveList_IS) error('No plane wave list present. Activate using RavenProject.setExportPlaneWaveList(1).'); end if numel(obj.planeWaveList_IS) > 1 if nargin > 1 if nargin > 2 planewaves = obj.planeWaveList_IS{sourceID+1, receiverID+1}; % 2 arguments else planewaves = obj.planeWaveList_IS{sourceID+1, :}; % 1 argument end else planewaves = obj.planeWaveList_IS; % 0 argument - return whole cell array end else planewaves = obj.planeWaveList_IS{1,1}; % 0 argument - return single entry end end %------------------------------------------------------------------ function planewaves = getPlaneWaves_RT(obj, sourceID, receiverID) % % planewaves = getPlaneWaves_RT(sourceID, receiverID) % % Returns the plane waves of the ray tracing simulation. % The returned freqData values are linear absolute sound pressure. % sourceID : optional ID of the sound source % receiverID : optional ID of the receiver % (counting starts with zero) if isempty(obj.planeWaveList_RT) error('No plane wave list present. Activate using RavenProject.setExportPlaneWaveList(1).'); end if numel(obj.planeWaveList_RT) > 1 if nargin > 1 if nargin > 2 planewaves = obj.planeWaveList_RT{sourceID+1, receiverID+1}; % 2 arguments else planewaves = obj.planeWaveList_RT{sourceID+1, :}; % 1 argument end else planewaves = obj.planeWaveList_RT; % 0 argument - return whole cell array end else planewaves = obj.planeWaveList_RT{1,1}; % 0 argument - return single entry end end %------------------------------------------------------------------ function histo = buildHybridHistogramFromWallHitLogs(obj) % evaluate ray tracer wall hit log files and build histogram histo = obj.buildRTHistogramFromWallHitLog(); % evaluate image sources wall hit log files and add data to ray tracing histogram % histo = obj.addISEnergyToRTHistogram(histo); disp('IS disabled so far...'); end %------------------------------------------------------------------ function histo = buildRTHistogramFromWallHitLog(obj) % check if data is present if numel(obj.wallHitLog_RT) < 1 error('No wall hit log present.'); end % number of time steps numTimeSteps = ceil((obj.filterLength + 100)/ obj.timeSlotLength + 0.5); % calculate sound speed soundspeed = calculateSoundSpeed(obj.getTemperature, obj.getHumidity, obj.getPressure); % get mean particle running distance for each timeslot distance = (1:numTimeSteps)' * (0.5*obj.timeSlotLength/1000) * soundspeed; % calculate air absorption airAbsorption = determineAirAbsorptionParameter(obj.getTemperature, obj.getPressure, obj.getHumidity); % get number of sources and receivers and frequecy bands numSources = size(obj.wallHitLog_RT, 1); numReceivers = size(obj.wallHitLog_RT, 2); numFreqBands = size(obj.wallHitLog_RT, 3); % build material database roomMaterialNames = obj.getRoomMaterialNames(); numMaterials = numel(roomMaterialNames); materialMatrixAbsorption = zeros(numMaterials, 31); materialMatrixScattering = zeros(numMaterials, 31); for iMat = 1 : numMaterials [abso, scat] = readRavenMaterial(roomMaterialNames{iMat}); if ~isempty(abso) materialMatrixAbsorption(iMat,:) = abso; materialMatrixScattering(iMat,:) = scat; else error('Material not found.'); end end % prepare histograms histo = cell(numSources, numReceivers); for iSrc = 1 : numSources for iRec = 1 : numReceivers histo{iSrc,iRec}.sourcename = obj.sourceNames{iSrc}; histo{iSrc,iRec}.receivername = obj.receiverNames{iSrc}; histo{iSrc,iRec}.timevector = (0:(numTimeSteps-1))' * obj.timeSlotLength; histo{iSrc,iRec}.data = zeros(numTimeSteps, numFreqBands); for iBand = 1 : numFreqBands % translate octave/third-octave indices if numFreqBands == 10 thirdOctIndex = (iBand-1)*3+3; % transfer functions, air absorption etc. is always in third octave resolution. this translates octaves to third octaves if necessary else thirdOctIndex = iBand; end % distance: Vector (numTimeslots x 1) vorher definieren mit der jeweiligen % mittleren distanz eines zeitslots % mean energy loss due to air absorption for each time slot airAbs = exp(- (airAbsorption(thirdOctIndex)) * distance); % use logarithm so that number of wall hits can be multiplied %M = exp(obj.wallHitLog_RT{iSrc,iRec,iBand} * log(1 - materialMatrixAbsorption(:,iBand))); M = bsxfun(@power, (1 - materialMatrixAbsorption(:,iBand)'), obj.wallHitLog_RT{iSrc,iRec,iBand}); M(sum(obj.wallHitLog_RT{iSrc,iRec,iBand},2) == 0) = 0; % elimante walls that haven't been hit M = prod(M, 2); % multiply all walls histo{iSrc,iRec}.data(:, iBand) = airAbs' .* sum(bsxfun(@times, reshape(M, obj.numParticles_Sphere, numTimeSteps), obj.initialParticleEnergy{iSrc,iBand})); % todo: mal ohne exp/log probieren... könnte schneller sein % % loop all detections % for iRefl = 1 : numel(obj.wallHitLog_RT{iSrc,iRec,iBand}) % % get particle running distance % distance = obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).time * soundspeed; % % % build wall filter % wallTransferFunction = 1; % for iWall = 1 : numel(obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).materials) % % wallTransferFunction = (1 - alpha) % wallTransferFunction = wallTransferFunction * ... % (1 - materialMatrixAbsorption(obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).materials(iWall)+1, thirdOctIndex)); % end % % % add energy % histo{iSrc,iRec}.data(obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).timeslot + 1, iBand) = ... % histo{iSrc,iRec}.data(obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).timeslot + 1, iBand) + ... % obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).initialEnergy * ... % wallTransferFunction * obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).lambertIntegral * ... % exp(- (airAbsorption(thirdOctIndex)) * distance); % % % % debug % % histo{iSrc,iRec}.data(obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).timeslot + 1, iBand) = ... % % histo{iSrc,iRec}.data(obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).timeslot + 1, iBand) + ... % % obj.wallHitLog_RT{iSrc,iRec,iBand}(iRefl).energy; % end end end end end %------------------------------------------------------------------ function histoRT = addISEnergyToRTHistogram(obj, histoRT) % this functions adds the energy of audible image sources to the spectral energy histogram % of the ray tracer. note: audible image sources and RT histogram must be determined before % calculate sound speed soundspeed = calculateSoundSpeed(obj.getTemperature, obj.getHumidity, obj.getPressure); % calculate air absorption airAbsorption = determineAirAbsorptionParameter(obj.getTemperature, obj.getPressure, obj.getHumidity); % get image sources wall hit log refl_IS = obj.getWallHitLog_IS(); if ~iscell(refl_IS) refl_IS = {refl_IS}; end % get number of sources and receivers and frequecy bands numSources = size(refl_IS, 1); numReceivers = size(refl_IS, 2); numFreqBands = numel(refl_IS{1,1}(1).spectrum); % build material database roomMaterialNames = obj.getRoomMaterialNames(); numMaterials = numel(roomMaterialNames); materialMatrixAbsorption = zeros(numMaterials, 31); materialMatrixScattering = zeros(numMaterials, 31); for iMat = 1 : numMaterials [materialMatrixAbsorption(iMat,:), materialMatrixScattering(iMat,:)] = readRavenMaterial(roomMaterialNames{iMat}); end % prepare histograms for iSrc = 1 : numSources % open directivity if ~isempty(obj.sourceDirectivity) && (numel(obj.sourceDirectivity) >= iSrc) sourceDirectivityFile = obj.sourceDirectivity{iSrc}; daffDirectivity = DAFF('open', sourceDirectivityFile); else sourceDirectivityFile = []; daffDirectivity = []; end for iRec = 1 : numReceivers % process direct sound and image sources for i = 1 : numel(refl_IS{iSrc,iRec}) % get time and distance of image source time = refl_IS{iSrc,iRec}(i).time; distance = time * soundspeed; histogramIndex = refl_IS{iSrc,iRec}(i).timeslot + 1; for band = 1 : numFreqBands if numFreqBands == 10 thirdOctIndex = (band-1)*3+3; % transfer functions, air absorption etc. is always in third octave resolution. this translates octaves to third octaves if necessary else thirdOctIndex = band; end % get source energy if isempty(sourceDirectivityFile) energy = 1.0; else % read DAFF directivity filter data = DAFF('getNearestNeighbourRecord', daffDirectivity, 'obj', refl_IS{iSrc,iRec}(i).azimuth, refl_IS{iSrc,iRec}(i).elevation); energy = data(thirdOctIndex) * data(thirdOctIndex); end % geometrical spreading energy = energy / (distance * distance); % build wall filter wallTransferFunction = 1; for iWall = 1 : numel(refl_IS{iSrc,iRec}(i).materials) % wallTransferFunction = (1 - alpha) * (1 - S) wallTransferFunction = wallTransferFunction * ... (1 - materialMatrixAbsorption(refl_IS{iSrc,iRec}(i).materials(iWall)+1, thirdOctIndex)) * ... (1 - materialMatrixScattering(refl_IS{iSrc,iRec}(i).materials(iWall)+1, thirdOctIndex)); end histoRT{iSrc,iRec}.data(histogramIndex, band) = histoRT{iSrc,iRec}.data(histogramIndex, band) + ... wallTransferFunction * energy * ... exp(- (airAbsorption(thirdOctIndex)) * distance); end % for number of bands end % for image sources % % build debug histogram % histodatadebug = zeros(size(histodata)); % for i = 1 : numel(refl_IS{iSrc,iRec}) % histodatadebug(refl_IS{iSrc,iRec}(i).timeslot, :) = histodatadebug(refl_IS{iSrc,iRec}(i).timeslot, :) + refl_IS{iSrc,iRec}(i).spectrum; % end end % for receivers end % for sources % close directivity file if ~isempty(sourceDirectivityFile) DAFF('close', daffDirectivity) end end %------------------------------------------------------------------ function V = getRoomVolume(obj, roomID) %V = getRoomVolume(roomID) % % returns the volume of the current room. % if nargin < 2 roomID = 0; end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % get all materials in room V = roommodel.getVolume(); end %------------------------------------------------------------------ function S = getRoomSurfaceArea(obj, roomID) %S = getRoomSurfaceArea(roomID) % % returns the total surface area of the current room. % if nargin < 2 roomID = 0; end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % get all materials in room S = roommodel.getSurface(); end %------------------------------------------------------------------ function S = getMaterialSurfaceArea(obj, materialName) %S = getRoomSurfaceArea(roomID) % % returns the surface area of the current material. % % WARNING: NOT IMPLEMENTED! % if nargin < 2 error('Please provide a material name.'); end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % return surface area of given material S = roommodel.getMaterialSurface(materialName); end %------------------------------------------------------------------ function A = getRoomEquivalentAbsorptionArea_Sabine(obj, roomID) %A = getRoomEquivalentAbsorptionArea_Sabine(roomID) % % returns the total equivalent absorption area of the current % room for all frequency bands. % if nargin < 2 roomID = 0; end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % get reverberation time after sabine A = roommodel.getEquivalentAbsorptionArea_sabine(obj.pathMaterials); end %------------------------------------------------------------------ function A = getRoomEquivalentAbsorptionArea_Eyring(obj, roomID) %A = getRoomEquivalentAbsorptionArea_Sabine(roomID) % % returns the total equivalent absorption area of the current % room for all frequency bands. % if nargin < 2 roomID = 0; end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % get reverberation time after eyring A = roommodel.getEquivalentAbsorptionArea_eyring(obj.pathMaterials); end %------------------------------------------------------------------ function RT = getEquationBasedReverbTime(obj, eyring, roomID) %RT = getEquationBasedReverbTime(eyring, roomID) % % returns the estimated reverberation time of the current room. % if nargin < 3 roomID = 0; end if nargin < 2 eyring = 1; % default to eyring end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % determine factor (default: 0.161) and air absorption % factor = 24*log(10) / calculateSoundSpeed(obj.getTemperature, obj.getHumidity,obj.getPressure); factor = 0.161; airAbsorption = determineAirAbsorptionParameter(obj.getTemperature, obj.getPressure, obj.getHumidity); % get reverberation time if (eyring == 1) RT = roommodel.getReverbTime(obj.pathMaterials, 'eyring', factor, airAbsorption); else RT = roommodel.getReverbTime(obj.pathMaterials, 'sabine', factor, airAbsorption); end end %------------------------------------------------------------------ function RT = getReverbTime_Sabine(obj) % getReverbTime_Sabine % returns the _calculated_ reverberation time based on the % volume, surface area and the absorption coefficients of the room % according to Sabine equation RT = obj.getEquationBasedReverbTime(0); end %------------------------------------------------------------------ function RT = getReverbTime_Eyring(obj) % getReverbTime_Eyring % returns the _calculated_ reverberation time based on the % volume, surface area and the absorption coefficients of the room % according to Eyring equation RT = obj.getEquationBasedReverbTime(1); end %------------------------------------------------------------------ function S = getSurfaceAreaOfMaterial(obj,material) if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end S = roommodel.getSurfaceArea(material); end %------------------------------------------------------------------ function materialNames = getRoomMaterialNames(obj, roomID) %materialNames = getRoomMaterialNames(roomID) % % returns the names of materials that are used in the % current room. % if nargin < 2 roomID = 0; end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % get all materials in room materialNames = roommodel.getMaterialNames(); end %------------------------------------------------------------------ function materialNames = setRoomMaterialNames(obj, materialNames, roomID) %setRoomMaterialNames(materialNames, roomID) % % returns the names of materials that are used in the % current room. % if nargin < 3 roomID = 0; end if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end % get all materials in room roommodel.setMaterialNames(materialNames); end %------------------------------------------------------------------ function IACC = getIACC(obj, brir) if nargin < 2 brir = obj.getBinauralImpulseResponseItaAudio(); end % apply filter bank third_octs = ita_mpb_filter(brir, 'oct', 3); timeData = third_octs.timeData; numFreqs = size(timeData, 2)/2; % preallocation IACC vector IACC = zeros(numFreqs, 1); samples_in_1ms = round(brir.samplingRate / 1000); wb=itaWaitbar(samples_in_1ms*2*31); for i3rd = 1 : numFreqs % calculate interaural correlation function IACF = sum(timeData(:, i3rd*2-1).^2) * sum(timeData(:, i3rd*2).^2); % search for maximum of cross correlation function in interval -1ms to +1ms % samples_in_1ms = round(brir.samplingRate / 1000); % shift left channel for lambda = 1 : samples_in_1ms this_iacc = abs(sum(timeData(lambda:end, i3rd*2-1) .* timeData(1:(end-lambda+1), i3rd*2)) / sqrt(IACF)); IACC(i3rd) = max([this_iacc, IACC(i3rd)]); wb.inc; end % shift right channel for lambda = 1 : samples_in_1ms this_iacc = abs(sum(timeData(1:(end-lambda+1), i3rd*2-1) .* timeData(lambda:end, i3rd*2)) / sqrt(IACF)); IACC(i3rd) = max([this_iacc, IACC(i3rd)]); wb.inc; end end wb.close; end %------------------------------------------------------------------ function [Ts] = getCenterTime(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end % calculate schroeder curve edc = obj.getSchroederCurve('nodB', 'nonorm', 'notimecorrect'); if ~iscell(edc) edc = {edc}; end % calculate sound speed c = calculateSoundSpeed(obj.getTemperature(), obj.getHumidity(), obj.getPressure()); for iRec = 1 : numel(edc) if ~isempty(edc{sourceID+1,iRec}) for iFreq = 1 : size(edc{sourceID+1,iRec}, 2) % calculate source to receiver distance and arrival time of direct sound srcPos = obj.getSourcePosition(sourceID); recPos = obj.getReceiverPosition(iRec-1); source_receiver_distance = norm(recPos - srcPos); directSoundTime = source_receiver_distance / c; % create timevector with center times of each time step timeVecDirect = obj.histogram{sourceID+1,iRec}.timevector + (obj.timeSlotLength/1000) / 2; % set t=0 to arrival of the direct sound timeVecDirect = timeVecDirect - directSoundTime; % set negative entries (before direct sound arrival to zero) timeVecDirect(timeVecDirect < 0) = 0; % find timestep AFTER direct sound timestepDirect = floor(directSoundTime / (obj.timeSlotLength/1000)) + 1; % second timestep after direct sound (counting starts at 1 in matlab!) % time in the middle between direct sound arrival and end of the timestep timeFirstTimestep = (timestepDirect * (obj.timeSlotLength/1000) - directSoundTime) / 2; % calculcate numerator tmpIntegral = timeFirstTimestep * obj.histogram{sourceID+1,iRec}.data(timestepDirect, iFreq); % first time step tmpIntegral = tmpIntegral + sum(timeVecDirect((timestepDirect+1):end) .* obj.histogram{sourceID+1,iRec}.data((timestepDirect+1):end, iFreq)); % all remaining time steps % total energy in the impulse response totalEnergy = edc{sourceID+1,iRec}(1, iFreq); % calculate center time if totalEnergy > 0 Ts{iRec}(iFreq) = tmpIntegral / totalEnergy; else Ts{iRec}(iFreq) = 0; end end else Ts{iRec} = []; end end if averageOverReceivers Ts = obj.averageOverReceivers(Ts); end if averageOverFrequencies if nargin < 4 afterDIN = 0; end Ts = obj.averageAfterDIN(Ts, afterDIN); end if (numel(Ts) == 1) && iscell(Ts) Ts = Ts{1}; end end %------------------------------------------------------------------ function [C50, C80] = getClarity(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end % calculate schroeder curve edc = obj.getSchroederCurve('nodB', 'nonorm', 'notimecorrect'); if ~iscell(edc) edc = {edc}; end % calculate sound speed c = calculateSoundSpeed(obj.getTemperature(), obj.getHumidity(), obj.getPressure()); for iRec = 1 : numel(edc) if ~isempty(edc{sourceID+1,iRec}) for iFreq = 1 : size(edc{sourceID+1,iRec}, 2) % calculate source to receiver distance and arrival time of direct sound srcPos = obj.getSourcePosition(sourceID); recPos = obj.getReceiverPosition(iRec-1); source_receiver_distance = norm(recPos - srcPos); directSoundTime = source_receiver_distance / c; % find the exact integration limit in the histogram integral_end_time_50 = (directSoundTime + 0.050); integral_end_time_80 = (directSoundTime + 0.080); last_time_slot_50 = find(obj.histogram{sourceID+1,iRec}.timevector >= integral_end_time_50, 1) - 1; last_time_slot_80 = find(obj.histogram{sourceID+1,iRec}.timevector >= integral_end_time_80, 1) - 1; relativePortionOfLastTimeSlot_50 = rem(integral_end_time_50, obj.timeSlotLength/1000) / (obj.timeSlotLength/1000); relativePortionOfLastTimeSlot_80 = rem(integral_end_time_80, obj.timeSlotLength/1000) / (obj.timeSlotLength/1000); % detect energy in last time slot energyInLastTimeSlot_50 = edc{sourceID+1,iRec}(last_time_slot_50, iFreq) - edc{sourceID+1,iRec}(last_time_slot_50 + 1, iFreq); energyInLastTimeSlot_80 = edc{sourceID+1,iRec}(last_time_slot_80, iFreq) - edc{sourceID+1,iRec}(last_time_slot_80 + 1, iFreq); % total energy in the impulse response totalEnergy = edc{sourceID+1,iRec}(1, iFreq); % energy in the first 50ms/80ms after the direct sound energy50ms = totalEnergy - edc{sourceID+1,iRec}(last_time_slot_50, iFreq) + relativePortionOfLastTimeSlot_50 * energyInLastTimeSlot_50; energy80ms = totalEnergy - edc{sourceID+1,iRec}(last_time_slot_80, iFreq) + relativePortionOfLastTimeSlot_80 * energyInLastTimeSlot_80; % calculate clarity if totalEnergy > energy50ms C50{iRec}(iFreq) = 10 * log10( energy50ms / (totalEnergy-energy50ms) ); else C50{iRec}(iFreq) = inf; end if nargout > 1 if totalEnergy > energy80ms C80{iRec}(iFreq) = 10 * log10( energy80ms / (totalEnergy-energy80ms) ); else C80{iRec}(iFreq) = inf; end end end else C50{iRec} = []; if nargout > 1 C80{iRec} = []; end end end if averageOverReceivers C50 = obj.averageOverReceivers(C50); if nargout > 1 C80 = obj.averageOverReceivers(C80); end end if averageOverFrequencies if nargin < 4 afterDIN = 0; end C50 = obj.averageAfterDIN(C50, afterDIN); if nargout > 1 C80 = obj.averageAfterDIN(C80, afterDIN); end end if (numel(C50) == 1) && iscell(C50) C50 = C50{1}; if nargout > 1 C80 = C80{1}; end end end %------------------------------------------------------------------ function [D50, D80] = getDefinition(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end % get clarity [C50, C80] = obj.getClarity(averageOverReceivers, averageOverFrequencies, afterDIN, sourceID); % calculate definition from clarity if iscell(C50) for iCell = 1 : numel(C50) D50{iCell} = 1 - 1 ./ (10.^(C50{iCell}./10) + 1); if nargout > 1 D80{iCell} = 1 - 1 ./ (10.^(C80{iCell}./10) + 1); end end else D50 = 1 - 1 ./ (10.^(C50./10) + 1); if nargout > 1 D80 = 1 - 1 ./ (10.^(C80./10) + 1); end end if (numel(D50) == 1) && iscell(D50) D50 = D50{1}; if nargout > 1 D80 = D80{1}; end end end %------------------------------------------------------------------ function [LFC, LF, LE] = getLFC(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if ~isfield(obj.histogram{1,1}, 'energyLFC') error('LFC not implemented. Please update your RAVEN installation.'); end if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end % calculate schroeder curve edc = obj.getSchroederCurve('nodB', 'nonorm', 'notimecorrect'); if ~iscell(edc) edc = {edc}; end for iRec = 1 : numel(edc) % check if this receiver looks at the source srcPos = obj.getSourcePosition(sourceID); recPos = obj.getReceiverPosition(iRec-1); recView = obj.getReceiverViewVectors(iRec-1); directionErrorDeg = acosd(dot(srcPos-recPos, recView) / (norm(srcPos - recPos) * norm(recView))); if (directionErrorDeg > 2) warning(['Receiver ' num2str(iRec) ' does not look directly at the source! Re-orient the receiver to calculate LF(C).']); end if ~isempty(edc{sourceID+1,iRec}) for iFreq = 1 : size(edc{sourceID+1,iRec}, 2) % get arrival time of direct sound directSoundTime = obj.histogram{sourceID+1,iRec}.directsounddelay; % find the exact integration limit in the histogram integral_end_time = (directSoundTime + 0.080); last_time_slot = find(obj.histogram{sourceID+1,iRec}.timevector >= integral_end_time, 1) - 1; relativePortionOfLastTimeSlot = rem(integral_end_time, obj.timeSlotLength/1000) / (obj.timeSlotLength/1000); % detect energy in last time slot energyInLastTimeSlot = edc{sourceID+1,iRec}(last_time_slot, iFreq) - edc{sourceID+1,iRec}(last_time_slot + 1, iFreq); % total energy in the impulse response totalEnergy = edc{sourceID+1,iRec}(1, iFreq); % omnidirectional energy in the first 0ms..80ms after the direct sound energy80ms = totalEnergy - edc{sourceID+1,iRec}(last_time_slot, iFreq) + relativePortionOfLastTimeSlot * energyInLastTimeSlot; % calculate lateral fraction [cosine] if energy80ms > realmin LFC{iRec}(iFreq) = obj.histogram{sourceID+1,iRec}.energyLFC(iFreq) / energy80ms; if nargout > 1 LF{iRec}(iFreq) = obj.histogram{sourceID+1,iRec}.energyLF(iFreq) / energy80ms; if nargout > 2 LE{iRec}(iFreq) = obj.histogram{sourceID+1,iRec}.energyLE(iFreq) / energy80ms; end end else LFC{iRec}(iFreq) = 0; if nargout > 1 LF{iRec}(iFreq) = 0; if nargout > 2 LE{iRec}(iFreq) = 0; end end end end else LFC{iRec} = []; if nargout > 1 LF{iRec} = []; if nargout > 2 LE{iRec} = []; end end end end if averageOverReceivers LFC = obj.averageOverReceivers(LFC); if nargout > 1 LF = obj.averageOverReceivers(LF); if nargout > 2 LE = obj.averageOverReceivers(LE); end end end if averageOverFrequencies if nargin < 4 afterDIN = 0; end LFC = obj.averageAfterDINLateral(LFC, afterDIN); if nargout > 1 LF = obj.averageAfterDINLateral(LF, afterDIN); if nargout > 1 LE = obj.averageAfterDINLateral(LE, afterDIN); end end end if (numel(LFC) == 1) && iscell(LFC) LFC = LFC{1}; if nargout > 1 LF = LF{1}; if nargout > 2 LE = LE{1}; end end end end %------------------------------------------------------------------ function G = getStrength(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end edc = obj.getSchroederCurve('nodB', 'nonorm', 'notimecorr'); if ~iscell(edc) edc = {edc}; end %% Important note: This only works if you do not use a source directivity file for iRec = 1 : size(edc,2) % Strength G G{iRec} = 10 * log10( edc{sourceID+1,iRec}(1,:) / (1/100) ); % 1/100 = energy of direct sound per frequency band for 10m distance end if averageOverReceivers G = obj.averageOverReceivers(G); end if averageOverFrequencies if nargin < 4 afterDIN = 0; end G = obj.averageAfterDIN(G, afterDIN); end if (numel(G) == 1) && iscell(G) G = G{1}; end end %------------------------------------------------------------------ function EDT = getEDT(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end EDT = obj.getRT(averageOverReceivers, averageOverFrequencies, afterDIN, 0, -10, sourceID); end %------------------------------------------------------------------ function T20 = getT20(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end T20 = obj.getRT(averageOverReceivers, averageOverFrequencies, afterDIN, -5, -25, sourceID); end %------------------------------------------------------------------ function T30 = getT30(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end T30 = obj.getRT(averageOverReceivers, averageOverFrequencies, afterDIN, -5, -35, sourceID); end %------------------------------------------------------------------ function T60 = getT60(obj, averageOverReceivers, averageOverFrequencies, afterDIN, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 sourceID = 0; end T60 = obj.getRT(averageOverReceivers, averageOverFrequencies, afterDIN, -5, -65, sourceID); end %------------------------------------------------------------------ function RT = getRT(obj, averageOverReceivers, averageOverFrequencies, afterDIN, from_dB, to_dB, sourceID) if nargin < 2 averageOverReceivers = 0; end if nargin < 3 averageOverFrequencies = 0; end if nargin < 4 afterDIN = 0; end if nargin < 5 from_dB = -5; % default to T30 end if nargin < 6 to_dB = -35; % default to T30 end if nargin < 7 sourceID = 0; % default to first source end % get schroeder curve edc = obj.getSchroederCurve('dB', 'norm', 'timecorrect'); if ~iscell(edc) edc = {edc}; end % get number of frequency bands numFrequencyBands = size(obj.histogram{find(~cellfun(@isempty, obj.histogram), 1)}.data, 2); RT = cell(size(obj.histogram,2),1); % T30{1..numReceivers}(1..numFrequencyBands) for iRec = 1 : size(obj.histogram, 2) if ~isempty(obj.histogram{sourceID+1,iRec}) % get slope of schroeder curves slope = ones(numFrequencyBands, 1); for iBand = 1 : numFrequencyBands % linear regression analysis start_index = find(edc{sourceID+1,iRec}(:, iBand) <= from_dB, 1); end_index = find(edc{sourceID+1,iRec}(:, iBand) <= to_dB, 1); if start_index ~= end_index slope_temp = polyfit(obj.histogram{sourceID+1,iRec}.timevector(start_index : end_index), ... edc{sourceID+1,iRec}(start_index : end_index, iBand), 1); slope(iBand) = slope_temp(1); else slope(iBand) = -Inf; end end % reverberation times RT{iRec} = -60 ./ slope; else RT{iRec} = []; end end if averageOverReceivers RT = obj.averageOverReceivers(RT); end if averageOverFrequencies if nargin < 4 afterDIN = 0; end RT = obj.averageAfterDIN(RT, afterDIN); end if (numel(RT) == 1) && iscell(RT) RT = RT{1}; end end %------------------------------------------------------------------ function newReverbTime = adjustAbsorptionToMatchReverbTime(obj, targetReverbTime, roomID, validationsimulation, materialPrefix, materialAppendix, materialIndexVector) %adjustAbsorptionToMatchReverbTime(targetReverbTime, roomID, validationsimulation, materialPrefix, materialAppendix, materialIndexVector) if isempty(obj.modelFileList) return end % check args if nargin < 3 roomID = 0; end % check if simulation was already run if isempty(obj.histogram) if ~obj.exportHistogram obj.setExportHistogram(1); % make sure we get a histogram obj.run(); obj.setExportHistogram(0); else obj.run(); end else disp('Using results from last simulation run.'); end % check back the reverberation times numSources=size(obj.getSourcePosition,1); thisReverbTime = obj.getT30(true,false,false,0); for iSources=2:numSources thisReverbTime = thisReverbTime + obj.getT30(true,false,false,iSources-1); end thisReverbTime = thisReverbTime / numSources; % if desired, rename the materials of the model if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end materialNames = roommodel.getMaterialNames(); % read materials out of the ac3d model if nargin > 4 newMaterialNames = strcat(materialPrefix, materialNames); if nargin > 5 newMaterialNames = strcat(newMaterialNames, materialAppendix); end roommodel.setMaterialNames(newMaterialNames); else newMaterialNames = materialNames; end % calculate new absorption values for the variable model to match the reverberation time of the target model (default values by eyring including air abs) [A, S] = roommodel.getEquivalentAbsorptionArea_eyring(obj.pathMaterials); if (obj.getAirAbsorptionEnabled) airAbscoeffs = determineAirAbsorptionParameter(obj.getTemperature, obj.getPressure, obj.getHumidity); else airAbscoeffs = zeros(1,31); end equivalentAirAbsorptionArea = 4* roommodel.getVolume() * airAbscoeffs; if numel(targetReverbTime) == 10 % go to octave resolution A = mean([A(2:3:end);A(3:3:end); A(4:3:end)]); % A = A(3:3:end); equivalentAirAbsorptionArea = mean([equivalentAirAbsorptionArea(2:3:end);equivalentAirAbsorptionArea(3:3:end); equivalentAirAbsorptionArea(4:3:end)]); % equivalentAirAbsorptionArea = equivalentAirAbsorptionArea(3:3:end); end alphas_default = 1 - exp(((-0.163 * roommodel.getVolume() ./ targetReverbTime) - equivalentAirAbsorptionArea)/(S)); % alphas_alt = min((A+equivalentAirAbsorptionArea)/S,1); alphas_alt = (A+equivalentAirAbsorptionArea)/S; alphas_neu = abs(1 - (1 - alphas_alt(:)).^(thisReverbTime(:) ./ targetReverbTime(:))); absorptionFactors = alphas_neu(:) ./ alphas_alt(:); absorptionFactors(isnan(absorptionFactors)) = 1; absorptionFactors(isinf(absorptionFactors)) = 999; % apply new absorption values matProcessed={}; if nargin < 7 materialIndexVector = 1 : numel(materialNames); else materialIndexVector = materialIndexVector + 1; % C++ starts counting at zero, but MATLAB starts at 1 end for iMat = materialIndexVector if ~sum(strcmp(materialNames{iMat},matProcessed)) matProcessed{iMat}=materialNames{iMat}; [absorp, scatter] = readRavenMaterial(materialNames{iMat}, obj.pathMaterials); if isempty(absorp) absorp = alphas_default; scatter = 0.2; disp(['Material (' materialNames{iMat} ') not found. Using default absorption values by eyring.']); elseif numel(absorptionFactors) == 10 absorp = absorp(3:3:end); % go to octave resolution scatter = scatter(3:3:end); % go to octave resolution end absorp_neu = absorp(:) .* absorptionFactors(:); absorp_neu(absorp_neu > 1) = 1; absorp_neu(absorp_neu <= 0) = 0.0001; obj.setMaterial(newMaterialNames{iMat}, absorp_neu, scatter); end end if (nargin < 4) || validationsimulation disp('repeat simulation for validation'); obj.run(); % check back the reverberation times newReverbTime = obj.getT30(true,false,false,0); for iSources=2:numSources newReverbTime = newReverbTime + obj.getT30(true,false,false,iSources-1); end newReverbTime = newReverbTime / numSources; figure; nRT = numel(targetReverbTime); plot(1:nRT, targetReverbTime, 1:nRT, thisReverbTime, 1:nRT, newReverbTime); legend('Target RT', 'RT before', 'RT after'); ylabel('T30 [s]'); xlabel('Frequency band [Hz]'); ylim([0 1.2*max(targetReverbTime)]); ax = gca; ax.XTickLabel = {32 63 125 250 500 1000 2000 4000 8000 16000}; grid on; else newReverbTime = []; end end %------------------------------------------------------------------ function [ newReverbTime absorp_neu_all ]= adjustAbsorptionToMatchParameter(obj, targetReverbTime, roomID, validationsimulation, materialPrefix, materialAppendix, materialIndexVector) if isempty(obj.modelFileList) return end % check args if nargin < 3 roomID = 0; end % check if simulation was already run if isempty(obj.histogram) if ~obj.exportHistogram obj.setExportHistogram(1); % make sure we get a histogram obj.run(); obj.setExportHistogram(0); else obj.run(); end else disp('Using results from last simulation run.'); end % check back the reverberation times thisReverbTime = obj.getT30(1); % average over receivers % if desired, rename the materials of the model if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end materialNames = roommodel.getMaterialNames(); % read materials out of the ac3d model if nargin > 4 newMaterialNames = strcat(materialPrefix, materialNames); if nargin > 5 newMaterialNames = strcat(newMaterialNames, materialAppendix); end roommodel.setMaterialNames(newMaterialNames); else newMaterialNames = materialNames; end % calculate new absorption values for the variable model to match the reverberation time of the target model (default values by eyring including air abs) [A, S] = roommodel.getEquivalentAbsorptionArea_eyring(obj.pathMaterials); if (obj.getAirAbsorptionEnabled) airAbscoeffs = determineAirAbsorptionParameter(obj.getTemperature, obj.getPressure, obj.getHumidity); else airAbscoeffs = zeros(1,31); end equivalentAirAbsorptionArea = 4* roommodel.getVolume() * airAbscoeffs; if numel(targetReverbTime) == 10 % go to octave resolution A = A(3:3:end); equivalentAirAbsorptionArea = equivalentAirAbsorptionArea(3:3:end); end alphas_default = 1 - exp(((-0.163 * roommodel.getVolume() ./ targetReverbTime) - equivalentAirAbsorptionArea)/(S)); alphas_alt = (A+equivalentAirAbsorptionArea)/S; alphas_neu = 1 - (1 - alphas_alt(:)).^(thisReverbTime(:) ./ targetReverbTime(:)); absorptionFactors = alphas_neu(:) ./ alphas_alt(:); absorptionFactors(isnan(absorptionFactors)) = 1; absorptionFactors(isinf(absorptionFactors)) = 999; % apply new absorption values matProcessed={}; if nargin < 7 materialIndexVector = 1 : numel(materialNames); else materialIndexVector = materialIndexVector + 1; % C++ starts counting at zero, but MATLAB starts at 1 end absorp_neu_all = zeros(length(alphas_neu),numel(materialIndexVector)); for iMat = materialIndexVector if ~sum(strcmp(materialNames{iMat},matProcessed)) matProcessed{iMat}=materialNames{iMat}; [absorp, scatter] = readRavenMaterial(materialNames{iMat}, obj.pathMaterials); if isempty(absorp) absorp = alphas_default; scatter = 0.2; disp(['Material (' materialNames{iMat} ') not found. Using default absorption values by eyring.']); elseif numel(absorptionFactors) == 10 absorp = absorp(3:3:end); % go to octave resolution scatter = scatter(3:3:end); % go to octave resolution end absorp_neu = absorp(:) .* absorptionFactors(:); absorp_neu(absorp_neu > 1) = 1; absorp_neu(absorp_neu <= 0) = 0.0001; obj.setMaterial(newMaterialNames{iMat}, absorp_neu, scatter); absorp_neu_all(:,iMat) = absorp_neu; end end if (nargin < 4) || validationsimulation disp('repeat simulation for validation'); obj.run(); % check back the reverberation times newReverbTime = obj.getT30(1); % figure; % nRT = numel(targetReverbTime); % plot(1:nRT, targetReverbTime, 1:nRT, thisReverbTime, 1:nRT, newReverbTime); % legend('Target RT', 'RT before', 'RT after'); % % ylabel('T30 [s]'); % xlabel('Frequency band [Hz]'); % ylim([0 1.1*max(targetReverbTime)]); % % ax = gca; % ax.XTickLabel = {32 63 125 250 500 1000 2000 4000 8000 16000}; % grid on; else newReverbTime = []; end end %------------------------------------------------------------------ function brir = adjustAbsorptionToMatchReverbTimeOfBinauralImpulseResponse(obj, targetReverbTime, roomID, validationsimulation, materialPrefix, materialAppendix) %adjustAbsorptionToMatchReverbTimeOfBinauralImpulseResponse(targetReverbTime, roomID, validationsimulation) if isempty(obj.modelFileList) return end % check args if nargin < 3 roomID = 0; end % check if simulation was already run if isempty(obj.binauralIR) obj.run(); else disp('Using results from last simulation run.'); end disp('get binaural impulse response'); brir = obj.getBinauralImpulseResponseItaAudio(); % merge IRs of multiple receivers into 1 multichannel itaAudio if iscell(brir) tmp = brir{1}; for i = 2 : numel(brir) tmp = merge(tmp, brir{i}); end brir = tmp; elseif numel(brir) > 1 % multi-instance brir = merge(brir); end % calculate schroeder decay curves and get reverberation times thisReverbTime = ita_roomacoustics(brir, 'T20', 'edcMethod','noCut', 'freqRange', [20 20000], 'bandsPerOctave', 3); thisReverbTime = thisReverbTime.T20; thisReverbTime = mean(thisReverbTime.freqData, 2); % check if target reverb times are values already or if they still have to be calculated from an itaAudio if isa(targetReverbTime, 'itaAudio') targetReverbTime = ita_roomacoustics(targetReverbTime, 'T20', 'edcMethod','noCut', 'freqRange', [20 20000], 'bandsPerOctave', 3); targetReverbTime = targetReverbTime.T20; targetReverbTime = mean(targetReverbTime.freqData, 2); end disp('load room model'); % if desired, rename the materials of the model if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end materialNames = roommodel.getMaterialNames(); % read materials out of the ac3d model if nargin > 4 newMaterialNames = strcat(materialPrefix, materialNames); if nargin > 5 newMaterialNames = strcat(newMaterialNames, materialAppendix); end roommodel.setMaterialNames(newMaterialNames); else newMaterialNames = materialNames; end disp('calculate new absorption values for the variable model to match the reverberation time of the target model'); [A, S] = roommodel.getEquivalentAbsorptionArea(); if numel(targetReverbTime) == 10 A = A(3:3:end); % go to octave resolution end alphas_default = 1 - exp(-0.163 * roommodel.getVolume() ./ (S .* targetReverbTime)); alphas_alt = A/S; alphas_neu = 1 - (1 - alphas_alt(:)).^(thisReverbTime(:) ./ targetReverbTime(:)); absorptionFactors = alphas_neu(:) ./ alphas_alt(:); absorptionFactors(isnan(absorptionFactors)) = 1; absorptionFactors(isinf(absorptionFactors)) = 999; disp('apply new absorption values'); matProcessed={}; for iMat = 1 : numel(materialNames) if ~sum(strcmp(materialNames{iMat},matProcessed)) matProcessed{iMat}=materialNames{iMat}; [absorp, scatter] = readRavenMaterial(materialNames{iMat}, obj.pathMaterials); if isempty(absorp) absorp = alphas_default; scatter = 0.2; disp(['Material (' materialNames{iMat} ') not found. Using default absorption values by eyring.']); end if numel(absorptionFactors) == 10 absorp = absorp(3:3:end); % go to octave resolution scatter = scatter(3:3:end); % go to octave resolution end absorp_neu = absorp(:) .* absorptionFactors(:); absorp_neu(absorp_neu > 1) = 1; absorp_neu(absorp_neu <= 0) = 0.0001; obj.setMaterial(newMaterialNames{iMat}, absorp_neu, scatter); end end if (nargin < 4) || validationsimulation disp('repeat simulation for validation'); obj.run(); disp('generate binaural impulse response'); brir = obj.getBinauralImpulseResponseItaAudio(); % merge IRs of multiple receivers into 1 multichannel itaAudio if numel(brir) > 1 brir = merge(brir); end % calculate schroeder decay curves and get reverberation times newReverbTime = ita_roomacoustics(brir, 'T20', 'edcMethod','noCut', 'freqRange', [20 20000], 'bandsPerOctave', 3); newReverbTime = newReverbTime.T20; newReverbTime = mean(newReverbTime.freqData, 2); figure; nRT = numel(targetReverbTime); plot(1:nRT, targetReverbTime, 1:nRT, thisReverbTime, 1:nRT, newReverbTime); legend('Target RT', 'RT before', 'RT after'); end end %------------------------------------------------------------------ function rir = adjustAbsorptionToMatchReverbTimeOfMonauralImpulseResponses(obj, targetReverbTime, roomID, validationsimulation, materialPrefix, materialAppendix, materialIndexVector) %adjustAbsorptionToMatchReverbTimeOfMonauralImpulseResponses(targetReverbTime, roomID, validationsimulation, materialPrefix, materialAppendix) if isempty(obj.modelFileList) return end % check args if nargin < 3 roomID = 0; end % check if simulation was already run if isempty(obj.monoIR) obj.run(); else disp('Using results from last simulation run.'); end disp('get mono room impulse response'); rir = obj.getImpulseResponseItaAudio(); % merge IRs of multiple receivers into 1 multichannel itaAudio if iscell(rir) tmp = rir{1}; for i = 2 : numel(rir) tmp = merge(tmp, rir{i}); end rir = tmp; elseif numel(rir) > 1 % multi-instance rir = merge(rir); end % calculate schroeder decay curves and get reverberation times if numel(targetReverbTime) == 10 thisReverbTime = ita_roomacoustics(rir, 'T20', 'edcMethod','noCut', 'freqRange', [30 16000], 'bandsPerOctave', 1); else thisReverbTime = ita_roomacoustics(rir, 'T20', 'edcMethod','noCut', 'freqRange', [20 20000], 'bandsPerOctave', 3); end thisReverbTime = mean(thisReverbTime.T20.freqData, 2); % check if target reverb times are values already or if they still have to be calculated from an itaAudio if isa(targetReverbTime, 'itaAudio') targetReverbTime = ita_roomacoustics(targetReverbTime, 'T20', 'edcMethod','noCut', 'freqRange', [20 20000], 'bandsPerOctave', 3); targetReverbTime = mean(targetReverbTime.freqData, 2); end disp('load room model'); % if desired, rename the materials of the model if isempty(obj.model) if iscell(obj.modelFileList) for iRoom = 1 : numel(obj.modelFileList) obj.model{iRoom} = itaAc3dModel(obj.modelFileList{iRoom}); end roommodel = obj.model{roomID + 1}; else obj.model = itaAc3dModel(obj.modelFileList); roommodel = obj.model; end else if iscell(obj.model) roommodel = obj.model{roomID + 1}; else roommodel = obj.model; end end materialNames = roommodel.getMaterialNames(); % read materials out of the ac3d model if nargin > 4 newMaterialNames = strcat(materialPrefix, materialNames); if nargin > 5 newMaterialNames = strcat(newMaterialNames, materialAppendix); end roommodel.setMaterialNames(newMaterialNames); else newMaterialNames = materialNames; end disp('calculate new absorption values for the variable model to match the reverberation time of the target model'); [A, S] = roommodel.getEquivalentAbsorptionArea(); if numel(targetReverbTime) == 10 A = A(3:3:end); % go to octave resolution end alphas_default = 1 - exp(-0.163 * roommodel.getVolume() ./ (S .* targetReverbTime)); alphas_alt = A/S; alphas_neu = 1 - (1 - alphas_alt(:)).^(thisReverbTime(:) ./ targetReverbTime(:)); absorptionFactors = alphas_neu(:) ./ alphas_alt(:); absorptionFactors(isnan(absorptionFactors)) = 1; absorptionFactors(isinf(absorptionFactors)) = 999; disp('apply new absorption values'); matProcessed={}; if nargin < 7 materialIndexVector = 1 : numel(materialNames); else materialIndexVector = materialIndexVector + 1; % C++ starts counting at zero, but MATLAB starts at 1 end for iMat = materialIndexVector if ~sum(strcmp(materialNames{iMat},matProcessed)) matProcessed{iMat}=materialNames{iMat}; [absorp, scatter] = readRavenMaterial(materialNames{iMat}, obj.pathMaterials); if isempty(absorp) absorp = alphas_default; scatter = 0.2; disp(['Material (' materialNames{iMat} ') not found. Using default absorption values by Eyring.']); end if numel(absorptionFactors) == 10 absorp = absorp(3:3:end); % go to octave resolution scatter = scatter(3:3:end); % go to octave resolution end absorp_neu = absorp(:) .* absorptionFactors(:); absorp_neu(absorp_neu > 1) = 1; absorp_neu(absorp_neu <= 0) = 0.0001; obj.setMaterial(newMaterialNames{iMat}, absorp_neu, scatter); end end if (nargin < 4) || validationsimulation disp('repeat simulation for validation'); obj.run(); disp('get monaural impulse response'); rir = obj.getImpulseResponseItaAudio(); % merge IRs of multiple receivers into 1 multichannel itaAudio if numel(rir) > 1 rir = merge(rir); end % calculate schroeder decay curves and get reverberation times if numel(targetReverbTime) == 10 newReverbTime = ita_roomacoustics(rir, 'T20', 'edcMethod','noCut', 'freqRange', [30 16000], 'bandsPerOctave', 1); else newReverbTime = ita_roomacoustics(rir, 'T20', 'edcMethod','noCut', 'freqRange', [20 20000], 'bandsPerOctave', 3); end newReverbTime = mean(newReverbTime.T20.freqData, 2); figure; nRT = numel(targetReverbTime); plot(1:nRT, targetReverbTime, 1:nRT, thisReverbTime, 1:nRT, newReverbTime); legend('Target RT', 'RT before', 'RT after'); end end %------------------------------------------------------------------ function T30 = getT30_fromImpulseResponse(obj) IRs = obj.getImpulseResponseItaAudio(); for i = 1 : numel(IRs) if ~isempty(IRs{i}) T30{i} = ita_roomacoustics_reverberation_time(IRs{i}); else T30{i} = []; end end if numel(T30) == 1 T30 = T30{1}; end end %------------------------------------------------------------------ function schroeder_edc = getSchroederCurve(obj, nodB, norm, timecorrect) if nargin < 3 timecorrect = 'timecorrect'; end if nargin < 3 norm = 'norm'; end if nargin < 2 nodB = 'dummy'; end % load histogram if isempty(obj.histogram) error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end % preallocation schroeder_edc = cell(size(obj.histogram, 1), size(obj.histogram, 2)); for iSrc = 1 : size(obj.histogram, 1) for iRec = 1 : size(obj.histogram, 2) if ~isempty(obj.histogram{iSrc,iRec}) % calculate schroeder function schroeder_edc{iSrc,iRec} = cumsum(obj.histogram{iSrc,iRec}.data(end:-1:1, :), 1); schroeder_edc{iSrc,iRec} = schroeder_edc{iSrc,iRec}(end:-1:1, :); % time correction? (direct sound should be at t = 0) if strcmp(timecorrect, 'timecorrect') directSoundEnergy = max(schroeder_edc{iSrc,iRec}); [timeInd, freqInd] = find(abs(bsxfun(@minus, schroeder_edc{iSrc,iRec}, directSoundEnergy)) < eps, 1, 'last'); schroeder_edc{iSrc,iRec}(1:timeInd, :) = []; end % normalize? if strcmp(norm, 'norm') schroeder_edc{iSrc,iRec} = bsxfun(@rdivide,schroeder_edc{iSrc,iRec},schroeder_edc{iSrc,iRec}(1, :)); end % dB scale? if ~strcmp(nodB, 'nodB') schroeder_edc{iSrc,iRec} = 10 * log10(schroeder_edc{iSrc,iRec} + eps); end else schroeder_edc{iSrc,iRec} = []; end end end if numel(schroeder_edc) == 1 schroeder_edc = schroeder_edc{1,1}; end end %------------------------------------------------------------------ function schroeder_edc = getSchroederCurve_itaResult(obj, nodB, norm, timecorrect) if nargin < 3 timecorrect = 'timecorrect'; end if nargin < 3 norm = 'norm'; end if nargin < 2 nodB = 'nodB'; end schroeder_data = obj.getSchroederCurve(nodB, norm, timecorrect); if ~iscell(schroeder_data) schroeder_data = {schroeder_data}; end schroeder_edc(size(obj.histogram,1), size(obj.histogram,2)) = itaResult(); for iHisto = 1 : numel(obj.histogram) if ~isempty(obj.histogram{iHisto}) schroeder_edc(iHisto).timeData = schroeder_data{iHisto}; schroeder_edc(iHisto).timeVector = obj.histogram{1}.timevector(1:size(schroeder_edc(iHisto).timeData, 1)); if obj.filterResolution == 1 schroeder_edc(iHisto).channelNames = strcat('EDC ', obj.freqLabelOct); else schroeder_edc(iHisto).channelNames = strcat('EDC ', obj.freqLabel3rd); end end end end %------------------------------------------------------------------ function wavedata = loadWaveFile(obj, filename) if isempty(filename) wavedata = []; return; end if ~iscell(filename) filename = {filename}; end % get number of sources and receivers numSources = numel(obj.sourcePositions) / 3; numReceivers = numel(obj.receiverPositions) / 3; % number of impulse responses numFiles = numel(filename); if (numFiles ~= numSources * numReceivers) warning(['There are ' num2str(numSources * numReceivers) ' Source->Receiver combinations, but ' num2str(numFiles) ' impulse response files.']); end info_file = IniConfig(); wavedata = cell(numSources, numReceivers); for i = 1 : numFiles % get source and receiver IDs info_file_name = [filename{i}(1:end-4) '_info.rir']; info_file.ReadFile(info_file_name); sourceID = info_file.GetValues('PrimarySource', 'sourceID', -1); receiverID = info_file.GetValues('Receiver', 'receiverID', -1); if (sourceID < 0) || (receiverID < 0) continue; end % LAS: wavread is not supported by Matlab versions newer than 2015a % [wavedata{sourceID+1, receiverID+1}, obj.sampleRate,~] = wavread(filename{i}); [wavedata{sourceID+1, receiverID+1}, obj.sampleRate] = audioread(filename{i}); sourcename = info_file.GetValuesWithoutDefaults('PrimarySource', 'sourceName'); receivername = info_file.GetValuesWithoutDefaults('Receiver', 'receiverName'); % read performance data if (obj.logPerformance) obj.loadPerformanceData(info_file, i) end % check matlab <-> raven consistency (via name tags) if ~strcmp(sourcename, obj.sourceNames{sourceID+1}) disp(['Warning! Probably wrong source<->IR mapping (source "' obj.sourceNames{sourceID+1} '" vs. IR "' sourcename '").']); end if ~strcmp(receivername, obj.receiverNames{receiverID+1}) disp(['Warning! Probably wrong receiver<->IR mapping (receiver "' obj.receiverNames{receiverID+1} '" vs. IR "' receivername '").']); end end end %------------------------------------------------------------------ function initPerformanceData(obj) % get number of sources and receivers numFiles = (numel(obj.sourcePositions) / 3)*(numel(obj.receiverPositions) / 3); % init variables obj.performance.ISFilterMonaural = zeros(1,numFiles); obj.performance.ISFilterBinaural = zeros(1,numFiles); obj.performance.RTFilterMonaural = zeros(1,numFiles); obj.performance.RTFilterBinaural = zeros(1,numFiles); obj.performance.ISGenerateImageSources = zeros(1,numFiles); obj.performance.ISTransformationMatrix = zeros(1,numFiles); obj.performance.ISAudibilityTest = zeros(1,numFiles); obj.performance.RTTotal = zeros(1,numFiles); if obj.filterResolution==1 nBands = 10; else nBands = 31; end obj.performance.RTBands = zeros(nBands,numFiles); end %------------------------------------------------------------------ function loadPerformanceData(obj, info_file, fileIndex) % NOTE: Only works for BSP (AccelerationType=0) % reads ravens performance data which is stored in the info files i = fileIndex; % infofiles for one simulation may contain redundant information. Only read info % file if no data is present yet (0) if (obj.performance.ISFilterMonaural(i) == 0) obj.performance.ISFilterMonaural(i) = info_file.GetValues('Performance', 'RF_Controller::computeMonauralImpulseResponseImageSources_PS2R',0); end if (obj.performance.ISFilterBinaural(i) == 0) obj.performance.ISFilterBinaural(i) = info_file.GetValues('Performance', 'RF_Controller::computeBinauralImpulseResponseImageSources_PS2R',0); end if (obj.performance.RTFilterMonaural(i) == 0) obj.performance.RTFilterMonaural(i) = info_file.GetValues('Performance', 'RF_Controller::computeMonauralImpulseResponseRayTracer_PS2R',0); end if (obj.performance.RTFilterBinaural(i) == 0) obj.performance.RTFilterBinaural(i) = info_file.GetValues('Performance', 'RF_Controller::computeBinauralImpulseResponseRayTracer_PS2R',0); end if (obj.performance.ISGenerateImageSources(i) == 0) obj.performance.ISGenerateImageSources(i) = info_file.GetValues('Performance', 'RIS_Ops::generate_image_sources_PS',0); end if (obj.performance.ISTransformationMatrix(i) == 0) obj.performance.ISTransformationMatrix(i) = info_file.GetValues('Performance', 'RIS_Ops::calculate_image_sources_transformation_matrix',0); end if (obj.performance.ISAudibilityTest(i) == 0) obj.performance.ISAudibilityTest(i) = info_file.GetValues('Performance', 'RIS_Ops::find_audible_sources_PS2R_BSP',0); end if (obj.performance.RTTotal(i) == 0) obj.performance.RTTotal(i) = info_file.GetValues('Performance', 'RRT_Controller::ExecuteRayTracing_PS2R_BSP',0); end if (sum(obj.performance.RTBands(:,i)) == 0) for l=1:size(obj.performance.RTBands,1) obj.performance.RTBands(l,i) = info_file.GetValues('Performance', ['RRT_Controller::ExecuteRayTracing_PS2R_BSP_band' int2str(l-1)],0); end end end %------------------------------------------------------------------ function histo = readHistogram(obj, histoFiles) % % [histogram, timevector] = readHistogram() % % Returns the data of the given histogram file(s) % if isempty(histoFiles) histo = []; return; end if ~iscell(histoFiles) histoFiles = {histoFiles}; end numHistograms = numel(histoFiles); numSources = numel(obj.sourcePositions) / 3; numReceivers = numel(obj.receiverPositions) / 3; histo = cell(numSources, numReceivers); hst = IniConfig(); for i_histo = 1 : numHistograms % Read INI file hst.ReadFile(histoFiles{i_histo}); % get source ID sourceID = hst.GetValues('Histogram', 'PrimarySourceID', -1); if sourceID < 0 continue; end % get receiver ID receiverID = hst.GetValues('Histogram', 'ReceiverID', -1); if receiverID < 0 continue; end histo{sourceID+1,receiverID+1}.sourcename = hst.GetValuesWithoutDefaults('Histogram', 'SourceName'); histo{sourceID+1,receiverID+1}.receivername = hst.GetValuesWithoutDefaults('Histogram', 'ReceiverName'); if ~strcmp(histo{sourceID+1,receiverID+1}.receivername, obj.receiverNames{receiverID+1}) disp(['Warning! Probably wrong receiver-histogram mapping (receiver "' obj.receiverNames{receiverID+1} '" vs. histogram "' histo{sourceID+1,receiverID+1}.receivername '").']); end % read direct sound time delay histo{sourceID+1,receiverID+1}.directsounddelay = hst.GetValuesWithoutDefaults('Histogram', 'directSoundTimeDelay'); % read lateral energy sums histo{sourceID+1,receiverID+1}.energyLE = hst.GetValuesWithoutDefaults('Histogram', 'energyLE'); histo{sourceID+1,receiverID+1}.energyLF = hst.GetValuesWithoutDefaults('Histogram', 'energyLF'); histo{sourceID+1,receiverID+1}.energyLFC = hst.GetValuesWithoutDefaults('Histogram', 'energyLFC'); % read time steps numTimeSteps = hst.GetValuesWithoutDefaults('Histogram', 'Timesteps'); % numFreqBands = hst.GetValues('Histogram', 'Bands'); % histo{receiverID+1}.data = zeros(numTimeSteps, numFreqBands); % histo{receiverID+1}.timevector = zeros(1, numTimeSteps); % for i = 1 : numTimeSteps % read_line = hst.GetValues('Histogram', num2str(i-1)); % histo{receiverID+1}.timevector(i) = read_line(1); % histo{receiverID+1}.data(i,:) = read_line(3:end); % end all_time_steps = cellstr(int2str((0:numTimeSteps-1)')); read_line = cell2mat(hst.GetFloatVector('Histogram', all_time_steps)); histo{sourceID+1,receiverID+1}.timevector = read_line(:, 1); histo{sourceID+1,receiverID+1}.data = read_line(:, 3:end); end end %------------------------------------------------------------------ function wall_log = readWallHitLog_IS(obj, wallFiles) % % wall_log = readWallHitLog_IS(obj, wallFiles) % % Returns the list of image source reflections with wall hit list % if isempty(wallFiles) wall_log = []; return; end if ~iscell(wallFiles) wallFiles = {wallFiles}; end numWallLogs = numel(wallFiles); % get number of sources and receivers if ~isempty(obj.sourcePositions) && numWallLogs > 1 numSources = numel(obj.sourcePositions) / 3; else numSources = 1; end if ~isempty(obj.receiverPositions) && numWallLogs > 1 numReceivers = numel(obj.receiverPositions) / 3; else numReceivers = 1; end % prealloc wall_log = cell(numSources,numReceivers); log = IniConfig(); for i_log = 1 : numWallLogs % Read INI file log.ReadFile(wallFiles{i_log}); % get source and receiver ID sourceID = log.GetValues('WallHitLog', 'sourceID', -1); receiverID = log.GetValues('WallHitLog', 'receiverID', -1); if (sourceID < 0) || (receiverID < 0) continue; end numReflections = log.GetValues('WallHitLog', 'numReflections', -1); % numFreqBands = log.GetValues('WallHitLog', 'numFreqBands', -1); % prealloc wall_log{sourceID+1,receiverID+1}(numReflections) = struct; for i = 1 : numReflections wall_log{sourceID+1,receiverID+1}(i).time = log.GetValuesWithoutDefaults(num2str(i-1), 'time'); wall_log{sourceID+1,receiverID+1}(i).timeslot = log.GetValuesWithoutDefaults(num2str(i-1), 'timeslot'); wall_log{sourceID+1,receiverID+1}(i).azimuth_PS = log.GetValuesWithoutDefaults(num2str(i-1), 'azimuth_PS'); wall_log{sourceID+1,receiverID+1}(i).elevation_PS = log.GetValuesWithoutDefaults(num2str(i-1), 'elevation_PS'); wall_log{sourceID+1,receiverID+1}(i).azimuth_Rec = log.GetValuesWithoutDefaults(num2str(i-1), 'azimuth_Rec'); wall_log{sourceID+1,receiverID+1}(i).elevation_Rec = log.GetValuesWithoutDefaults(num2str(i-1), 'elevation_Rec'); wall_log{sourceID+1,receiverID+1}(i).spectrum = log.GetValuesWithoutDefaults(num2str(i-1), 'spectrum'); wall_log{sourceID+1,receiverID+1}(i).materials = log.GetValuesWithoutDefaults(num2str(i-1), 'walls'); end end end %------------------------------------------------------------------ function [wall_log, initial_particle_energy] = readWallHitLog_RT(obj, wallFiles) % % wall_log = readWallHitLog_IS(obj, wallFiles) % % Returns the list of ray tracing reflections with wall hit list % if isempty(wallFiles) wall_log = []; return; end if ~iscell(wallFiles) wallFiles = {wallFiles}; end numWallLogs = numel(wallFiles); % get number of sources and receivers if ~isempty(obj.sourcePositions) && numWallLogs > 1 numSources = numel(obj.sourcePositions) / 3; else numSources = 1; end if ~isempty(obj.receiverPositions) && numWallLogs > 1 numReceivers = numel(obj.receiverPositions) / 3; else numReceivers = 1; end % get number of bands if obj.filterResolution == 1 numBands = 10; else numBands = 31; end % prealloc wall_log = cell(numSources,numReceivers,numBands); if nargout > 1 initial_particle_energy = cell(numSources,numBands); end for i_log = 1 : numWallLogs % Read source/receiver IDs and frequency band [logspath, fname, ~] = fileparts(wallFiles{i_log}); IDs = sscanf(fname, 'WallHitLog_PrimarySource[%i]_Receiver[%i]_Band[%i]_RT.log'); % get source and receiver ID sourceID = IDs(1); receiverID = IDs(2); bandID = IDs(3); % read data try wall_log{sourceID+1,receiverID+1,bandID+1} = dlmread(wallFiles{i_log}); catch disp(['Wall Log File empty for source ' num2str(sourceID) ' receiver ' num2str(receiverID) ' in band ' num2str(bandID)]); end wall_log{sourceID+1,receiverID+1,bandID+1} = sortrows(wall_log{sourceID+1,receiverID+1,bandID+1}); % %check for inital particle energy file % if nargout > 1 % initPartFileName = ['InitialParticleEnergy_PrimarySource[' num2str(sourceID) ']_Band[' num2str(bandID) ']_RT.log']; % if exist(fullfile(logspath, initPartFileName), 'file') % try % initial_particle_energy{sourceID+1,bandID+1} = dlmread(fullfile(logspath, initPartFileName)); % catch % disp(['Initial particle energy file empty or not found for source ' num2str(sourceID) ' receiver ' num2str(receiverID) ' in band ' num2str(bandID)]); % end % else % initial_particle_energy{sourceID+1,bandID+1} = []; % end % end % fid = fopen(wallFiles{i_log}); % reflectionIndex = 1; % while ~feof(fid) % readline = fgetl(fid); % data = sscanf(readline, '%f,'); % % % prealloc % %wall_log{sourceID+1,receiverID+1,bandID+1}(numReflections) = struct; % % wall_log{sourceID+1,receiverID+1,bandID+1}(reflectionIndex).initialEnergy = data(1); % wall_log{sourceID+1,receiverID+1,bandID+1}(reflectionIndex).time = data(2); % wall_log{sourceID+1,receiverID+1,bandID+1}(reflectionIndex).timeslot = floor(data(2) / (obj.timeSlotLength/1000)); % wall_log{sourceID+1,receiverID+1,bandID+1}(reflectionIndex).lambertIntegral = data(3); % wall_log{sourceID+1,receiverID+1,bandID+1}(reflectionIndex).energy = data(4); % wall_log{sourceID+1,receiverID+1,bandID+1}(reflectionIndex).materials = data(5:end); % % reflectionIndex = reflectionIndex + 1; % end % fclose(fid); end end %------------------------------------------------------------------ function planeWaveList = readPlaneWaveList(obj, planeWaveFiles) % % planeWaveList = readPlaneWaveList_RT(planeWaveFiles) % % Returns the list of reflections % if isempty(planeWaveFiles) planeWaveList = []; return; end if ~iscell(planeWaveFiles) planeWaveFiles = {planeWaveFiles}; end numPWLists = numel(planeWaveFiles); % get number of sources and receivers if ~isempty(obj.sourcePositions) && numPWLists > 1 numSources = numel(obj.sourcePositions) / 3; else numSources = 1; end if ~isempty(obj.receiverPositions) && numPWLists > 1 numReceivers = numel(obj.receiverPositions) / 3; else numReceivers = 1; end % get number of bands if obj.filterResolution == 1 numBands = 10; else numBands = 31; end % prealloc planeWaveList = cell(numSources, numReceivers); for i_log = 1 : numPWLists % Read source/receiver IDs and frequency band [logspath, fname, ~] = fileparts(planeWaveFiles{i_log}); IDs = sscanf(fname, 'PlaneWaves_PrimarySource[%i]_Receiver[%i]_%c%c.txt'); % get source and receiver ID sourceID = IDs(1); receiverID = IDs(2); % read data try tempData = dlmread(planeWaveFiles{i_log}); catch disp(['Plane Wave List empty for source ' num2str(sourceID) ' receiver ' num2str(receiverID)]); end % sort by time stemp tempData = sortrows(tempData); % get names planeWaveList{sourceID+1,receiverID+1}.sourcename = obj.sourceNames{sourceID+1}; planeWaveList{sourceID+1,receiverID+1}.receivername = obj.receiverNames{receiverID+1}; % arrival time planeWaveList{sourceID+1,receiverID+1}.time = tempData(:, 1); % angles planeWaveList{sourceID+1,receiverID+1}.azimuth = tempData(:, 2); planeWaveList{sourceID+1,receiverID+1}.elevation = tempData(:, 3); % frequency data planeWaveList{sourceID+1,receiverID+1}.freqData = tempData(:, 4:end); end end %------------------------------------------------------------------ function plotSphereEnergy(obj, sourceID) if nargin < 2 sourceID = 0; end if isempty(obj.histogram) if ~isempty(obj.monoIR) disp('Calculating energy from impulse response.'); sphereEnergy = zeros(size(obj.monoIR,2), 1); for i = 1 : size(obj.monoIR,2) sphereEnergy(i) = 10*log10( sum(obj.monoIR{sourceID+1,i}.^2) ); end else error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end else sphereEnergy = zeros(size(obj.histogram,2), 1); for i = 1 : size(obj.histogram,2) if isempty(obj.histogram{sourceID+1,i}) sphereEnergy(i) = -63; % -63 dB is the lower energy limit for particles in raven else sphereEnergy(i) = 10*log10(sum(sum(obj.histogram{sourceID+1,i}.data)) ); end end end % x = obj.receiverPositions(1:3:end); % y = obj.receiverPositions(2:3:end); % z = obj.receiverPositions(3:3:end); % scatter3(x, y, z, obj.radiusSphere * 1000, sphereEnergy, 'filled'); figure; for i = 1 : numel(sphereEnergy) [x,y,z] = sphere(16); x = x * obj.radiusSphere + -obj.receiverPositions(i, 3); y = y * obj.radiusSphere + obj.receiverPositions(i, 2); z = z * obj.radiusSphere + obj.receiverPositions(i, 1); c = ones(size(z,1), size(z,2)) * sphereEnergy(i); surf(z, x, y, c); hold on; end shading flat; set(gca, 'CameraViewAngleMode', 'manual'); axis equal; colorbar; hold on; obj.plotModel(gca); set(gca, 'CameraViewAngle', 10); end %------------------------------------------------------------------ function soundspeed = getSoundSpeed(obj) % Give temperature in degree celsius, relative humidity in percent and % pressure in Pa. Defaults are 20 degree celcius, 50% humidity and 101325 Pa % The exact solution is valid from 0 to 30 degrees celcius. % Taken from: % http:%resource.npl.co.uk/acoustics/techguides/speedair/: % "The calculator presented here computes the zero-frequency speed of sound % in humid air according to Cramer (J. Acoust. Soc. Am., 93, p2510, 1993), % with saturation vapour pressure taken from Davis, Metrologia, 29, p67, 1992, % and a mole fraction of carbon dioxide of 0.0004. % [...] temperature = obj.getTemperature(); humidity = obj.getHumidity(); pressure = obj.getPressure(); Kelvin = 273.15; % Measured ambient temp T_kel = Kelvin + temperature; % Molecular concentration of water vapour calculated from Rh using Giacomos method by Davis (1991) as implemented in DTU report 11b-1997 ENH = 3.14 * .00000001 * pressure + 1.00062 + temperature*temperature * 5.6 * .0000001; % These commented lines correspond to values used in Cramer (Appendix) % PSV1 = sqr(T_kel)*1.2811805*Math.pow(10,-5)-1.9509874*Math.pow(10,-2)*T_kel ; % PSV2 = 34.04926034-6.3536311*Math.pow(10,3)/T_kel; PSV1 = ( T_kel*T_kel * 1.2378847 * 0.00001 ) - ( 1.9121316 * 0.01 * T_kel ); PSV2 = 33.93711047 - 6.3431645 * 1000 / T_kel; PSV = exp(PSV1) * exp(PSV2); H = humidity * ENH * PSV / pressure; Xw = H / 100.0; % Xc = 314.0 * 10^-6; Xc = 400.0 * .000001; % Speed calculated using the method of Cramer from JASA vol 93 pg 2510 C1 = 0.603055 * temperature + 331.5024 - temperature*temperature * 5.28 / 10000 ... + (0.1495874 * temperature + 51.471935 - temperature*temperature * 7.82 / 10000) * Xw; C2 =(-1.82 / 10000000 + 3.73 / 100000000 * temperature - temperature*temperature * 2.93 / 10000000000) * pressure ... + (-85.20931 - 0.228525 * temperature + temperature*temperature * 5.91 / 100000 ) * Xc; C3 = Xw*Xw * 2.835149 + pressure*pressure * 2.15 / 10000000000000 ... - Xc*Xc * 29.179762 - 4.86 / 10000 * Xw * pressure * Xc; soundspeed = C1 + C2 - C3; end %------------------------------------------------------------------ function N = getNumberOfParticlesRecommendation(obj, stdDev_dB, roomID) % N = getNumberOfParticlesRecommendation(stdDev_dB, roomID) if nargin < 3 roomID = 0; end if nargin < 2 stdDev_dB = 1; % default 1dB as assumed JND end N = obj.getRoomVolume(roomID) ./ ((stdDev_dB/4.34)^2 * pi*obj.radiusSphere^2 * obj.getSoundSpeed() * obj.timeSlotLength/1000); end %------------------------------------------------------------------ function plotSurfaceEnergy(obj, sourceID) if nargin < 2 sourceID = 0; end if isempty(obj.uniformReceiverGridX) disp('For surface plots a uniform receiver grid has to be created first. (Or rebuilt using "rebuildReceiverGrid")'); disp('Trying to restore the grid from the given receiver positions...'); obj.restoreReceiverGrid(); end if isempty(obj.histogram) if ~isempty(obj.monoIR) disp('Calculating energy from impulse response.'); surfaceEnergy = zeros(size(obj.monoIR,2), 1); for i = 1 : size(obj.monoIR,2) surfaceEnergy(i) = 10*log10( sum(obj.monoIR{sourceID+1,i}.^2) ); end else error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end else sphereEnergy = zeros(size(obj.histogram,2), 1); for i = 1 : size(obj.histogram,2) if isempty(obj.histogram{sourceID+1,i}) sphereEnergy(i) = -63; % -63 dB is the lower energy limit for particles in raven else sphereEnergy(i) = 10*log10( sum(sum(obj.histogram{sourceID+1,i}.data)) ); end end end c = reshape(sphereEnergy, size(obj.uniformReceiverGridX,1), size(obj.uniformReceiverGridX,2)); surf(obj.uniformReceiverGridZ, obj.uniformReceiverGridX, obj.uniformReceiverGridY, c); shading flat; axis equal; colorbar; end %------------------------------------------------------------------ function plotSurfaceEnergyUniform(obj, resolution, sourceID) if nargin < 3 sourceID = 0; end if isempty(obj.histogram) if ~isempty(obj.monoIR) disp('Calculating energy from impulse response.'); surfaceEnergy = zeros(size(obj.monoIR,2), 1); for i = 1 : size(obj.monoIR,2) surfaceEnergy(i) = sum(obj.monoIR{sourceID+1,i}.^2); end else error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end else sphereEnergy = zeros(size(obj.histogram,2), 1); for i = 1 : size(obj.histogram,2) if isempty(obj.histogram{sourceID+1,i}) sphereEnergy(i) = -63; % -63 dB is the lower energy limit for particles in raven else sphereEnergy(i) = sum(sum(obj.histogram{sourceID+1,i}.data)); end end end if nargin < 2 resolution = 1; end x = obj.receiverPositions(:, 1); y = obj.receiverPositions(:, 2); z = obj.receiverPositions(:, 3); minx = min(x); maxx = max(x); % miny = min(y); % maxy = max(y); minz = min(z); maxz = max(z); c = reshape(sphereEnergy, size(x,1), size(x,2)); [uniformz, uniformx] = meshgrid(minz : resolution : maxz, minx : resolution : maxx); uniformy = uniformx; uniformy(:,:) = mean(mean(y)); uniformc = griddata(z, x, c, uniformz, uniformx); surf(uniformz, uniformx, uniformy, uniformc); shading flat; axis equal; colorbar; end %------------------------------------------------------------------ function plotSphereEnergyAnimation(obj, sourceID) if nargin < 2 sourceID = 0; end numTimeSteps = 0; if isempty(obj.histogram) if ~isempty(obj.monoIR) disp('Calculating energy from impulse response.'); numTimeSteps = size(obj.monoIR{find(~cellfun(@isempty, obj.monoIR), 1)}, 1); sphereEnergy = zeros(size(obj.monoIR,2), numTimeSteps); for i_sphere = 1 : size(obj.monoIR,2) for i_sample = 1 : numTimeSteps sphereEnergy(i_sphere, i_sample) = 10*log10(obj.monoIR{sourceID+1,i}(i_sample)^2 + eps); end end else error(sprintf('No histogram present.\nHINT: This error message might occur due to a corrupted SketchUp export.\nIf you are sure that the histogram is set, you could try to export your room model and raven project file again.')); %#ok end else numTimeSteps = numel(obj.histogram{find(~cellfun(@isempty, obj.histogram), 1)}.timevector); sphereEnergy = zeros(size(obj.histogram,2), numTimeSteps); for i_sphere = 1 : size(obj.histogram,2) if isempty(obj.histogram{sourceID+1,i_sphere}) sphereEnergy(i_sphere, :) = -63; else for i_timestep = 1 : numTimeSteps sphereEnergy(i_sphere, i_timestep) = 10*log10(sum(obj.histogram{sourceID+1,i_sphere}.data(i_timestep, :)) + eps); end end end end if numTimeSteps > 0 fig = figure; for i_sphere = 1 : size(sphereEnergy, 1) [x,y,z] = sphere(16); x = x * obj.radiusSphere + -obj.receiverPositions(i_sphere, 3); y = y * obj.radiusSphere + obj.receiverPositions(i_sphere, 2); z = z * obj.radiusSphere + obj.receiverPositions(i_sphere, 1); c = zeros(size(z,1), size(z,2)); sphere_handle(i_sphere) = surf(z, x, y, c); hold on; end shading flat; set(gca, 'CameraViewAngleMode', 'manual'); axis equal; caxis([-60 max(max(sphereEnergy))]); colorbar; hold on; obj.plotModel(gca); set(gca, 'CameraViewAngle', 10); for i_rep = 1 : 1 for i_time = 1 : numTimeSteps for i_sphere = 1 : size(sphereEnergy, 1) sphere_color = ones(size(x,1), size(x,2)) * sphereEnergy(i_sphere, i_time); set(sphere_handle(i_sphere), 'CData', sphere_color); end % show the current time in the impulse response if exist('h', 'var') delete(h); end h = text(0,0, [num2str(i_time/numTimeSteps * obj.timeSlotLength) 's']); pause(1 / i_time); % lighting phong; % set(gcf,'Renderer','zbuffer'); % frame(i_time) = getframe; end end close(fig); % figure; % movie(frame, 1, 10); end end %------------------------------------------------------------------ function getWallHitLogBand(obj,iBand) if obj.logPerformance obj.initPerformanceData; end disp('gatherResultsBand is called.'); if ~isnumeric(iBand) error('Your argument has to be numeric and positive scalar.'); end if obj.exportWallHitLog wall_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'WallHitLog_', ['_IS.log']}, '.log'); wall_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'WallHitLog_', ['[' num2str(iBand-1) ']_RT.log']}, '.log'); % read the wall hit logs back from disk obj.wallHitLog_IS = obj.readWallHitLog_IS(wall_files_IS); [obj.wallHitLog_RT, obj.initialParticleEnergy] = obj.readWallHitLog_RT(wall_files_RT); else obj.wallHitLog_IS = []; obj.wallHitLog_RT = []; obj.initialParticleEnergy = []; end if obj.exportPlaneWaveList planewave_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'PlaneWaves_', ['_IS.txt']}, '.txt'); planewave_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'PlaneWaves_', ['[' num2str(iBand-1) ']_RT.txt']}, '.txt'); % read the plane wave lists back from disk obj.planeWaveList_IS = obj.readPlaneWaveList(planewave_files_IS); obj.planeWaveList_RT = obj.readPlaneWaveList(planewave_files_RT); else obj.planeWaveList_IS = []; obj.planeWaveList_RT = []; end disp(['Successfully gathered results from band ' num2str(iBand)]); end end % public methods methods(Hidden = true) function gatherResultsFromSpecificFolder(obj, path) %Gathers the results from a specific folder. % Only use this if you know what you are doing. oldPath = obj.pathResults; obj.setOutputPath(path); obj.gatherResults(); obj.setOutputPath(oldPath); end %------------------------------------------------------------------ function monoIR = getMonauralImpulseResponse(obj) if isempty(obj.monoIR) error('No mono room impulse response present.'); end if numel(obj.monoIR) > 1 monoIR = obj.monoIR; else monoIR = obj.monoIR{1,1}; end end %------------------------------------------------------------------ function monoIRitaAudio = getMonauralImpulseResponseItaAudio(obj) if isempty(obj.monoIR) error('No mono room impulse response present.'); end monoIRitaAudio = itaAudio(size(obj.monoIR)); for iSrc = 1 : size(obj.monoIR, 1) for iRec = 1 : size(obj.monoIR, 2) monoIRitaAudio(iSrc,iRec).timeData = obj.monoIR{iSrc,iRec}; monoIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; monoIRitaAudio(iSrc,iRec).signalType = 'energy'; monoIRitaAudio(iSrc,iRec).channelNames{1} = obj.receiverNames{iRec}; %Make sure number of samples is even if mod(monoIRitaAudio(iSrc,iRec).nSamples,2)==1 monoIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function monoIR_IS = getMonauralImpulseResponseImageSources(obj) if isempty(obj.monoIR_IS) error('No mono image source impulse response present.'); end if numel(obj.monoIR_IS) > 1 monoIR_IS = obj.monoIR_IS; else monoIR_IS = obj.monoIR_IS{1,1}; end end %------------------------------------------------------------------ function monoIR_RT = getMonauralImpulseResponseRayTracing(obj) if isempty(obj.monoIR_RT) error('No mono ray tracing impulse response present.'); end if numel(obj.monoIR_RT) > 1 monoIR_RT = obj.monoIR_RT; else monoIR_RT = obj.monoIR_RT{1,1}; end end %------------------------------------------------------------------ function monoIRitaAudio = getMonauralImpulseResponseImageSourcesItaAudio(obj) if isempty(obj.monoIR_IS) error('No mono image source impulse response present.'); end monoIRitaAudio = itaAudio(size(obj.monoIR_IS)); for iSrc = 1 : size(obj.monoIR_IS, 1) for iRec = 1 : size(obj.monoIR_IS, 2) monoIRitaAudio(iSrc,iRec).timeData = obj.monoIR_IS{iSrc,iRec}; monoIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; monoIRitaAudio(iSrc,iRec).signalType = 'energy'; monoIRitaAudio(iSrc,iRec).channelNames{1} = obj.receiverNames{iRec}; %Make sure number of samples is even if mod(monoIRitaAudio(iSrc,iRec).nSamples,2)==1 monoIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end %------------------------------------------------------------------ function monoIRitaAudio = getMonauralImpulseResponseRayTracingItaAudio(obj) if isempty(obj.monoIR_RT) error('No mono ray tracing impulse response present.'); end monoIRitaAudio = itaAudio(size(obj.monoIR_RT)); for iSrc = 1 : size(obj.monoIR_RT, 1) for iRec = 1 : size(obj.monoIR_RT, 2) monoIRitaAudio(iSrc,iRec).timeData = obj.monoIR_RT{iSrc,iRec}; monoIRitaAudio(iSrc,iRec).samplingRate = obj.sampleRate; monoIRitaAudio(iSrc,iRec).signalType = 'energy'; monoIRitaAudio(iSrc,iRec).channelNames{1} = obj.receiverNames{iRec}; %Make sure number of samples is even if mod(monoIRitaAudio(iSrc,iRec).nSamples,2)==1 monoIRitaAudio(iSrc,iRec).timeData(end, :) = []; end end end end end %---------------------- PRIVATE METHODS ------------------------------% methods (Access = 'private') %------------------------------------------------------------------ function gatherResults(obj) if obj.logPerformance obj.initPerformanceData; end if obj.generateRIR ir_files = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\RIR_Combined', 'PrimarySource', 'Receiver'}, '.wav'); ir_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\RIR_IS', 'PrimarySource', 'Receiver'}, '.wav'); ir_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\RIR_RT', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.monoIR = obj.loadWaveFile(ir_files); obj.monoIR_IS = obj.loadWaveFile(ir_files_IS); obj.monoIR_RT = obj.loadWaveFile(ir_files_RT); else obj.monoIR = []; obj.monoIR_IS = []; obj.monoIR_RT = []; end if obj.generateBRIR ir_files = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\BRIR_Combined', 'PrimarySource', 'Receiver'}, '.wav'); ir_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\BRIR_IS', 'PrimarySource', 'Receiver'}, '.wav'); ir_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\BRIR_RT', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.binauralIR = obj.loadWaveFile(ir_files); obj.binauralIR_IS = obj.loadWaveFile(ir_files_IS); obj.binauralIR_RT = obj.loadWaveFile(ir_files_RT); else obj.binauralIR = []; obj.binauralIR_IS = []; obj.binauralIR_RT = []; end if obj.generateISHOA ir_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\HOA_IS', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.ambisonicsIR_IS = obj.loadWaveFile(ir_files_IS); % check for combined results with ray tracing if obj.generateRTHOA ir_files_combined = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\HOA_Combined', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.ambisonicsIR = obj.loadWaveFile(ir_files_combined); else obj.ambisonicsIR = []; end else obj.ambisonicsIR_IS = []; end if obj.generateRTHOA ir_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\HOA_RT', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.ambisonicsIR_RT = obj.loadWaveFile(ir_files_RT); else obj.ambisonicsIR_RT = []; end if obj.generateISVBAP ir_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\VBAP_IS', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.vbapIR_IS = obj.loadWaveFile(ir_files_IS); % check for combined results with ray tracing if obj.generateRTVBAP ir_files_combined = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\VBAP_Combined', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.vbapIR = obj.loadWaveFile(ir_files_combined); else obj.vbapIR = []; end else obj.vbapIR_IS = []; end if obj.generateRTVBAP ir_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'\VBAP_RT', 'PrimarySource', 'Receiver'}, '.wav'); % read the wave files back from disk obj.vbapIR_RT = obj.loadWaveFile(ir_files_RT); else obj.vbapIR_RT = []; end if obj.exportHistogram histo_files = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), '_Hybrid', '.hst'); raytracing_histo_files = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), '_RT', '.hst'); % read the histogram files back from disk obj.histogram = obj.readHistogram(histo_files); obj.histogramRT = obj.readHistogram(raytracing_histo_files); else obj.histogram = []; obj.histogramRT = []; end if obj.exportWallHitLog wall_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'WallHitLog_', '_IS.log'}, '.log'); wall_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'WallHitLog_', '_RT.log'}, '.log'); % read the wall hit logs back from disk obj.wallHitLog_IS = obj.readWallHitLog_IS(wall_files_IS); [obj.wallHitLog_RT, obj.initialParticleEnergy] = obj.readWallHitLog_RT(wall_files_RT); else obj.wallHitLog_IS = []; obj.wallHitLog_RT = []; obj.initialParticleEnergy = []; end if obj.exportPlaneWaveList planewave_files_IS = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'PlaneWaves_', '_IS.txt'}, '.txt'); planewave_files_RT = obj.scan_output_folder(fullfile(obj.pathResults, obj.projectTag), {'PlaneWaves_', '_RT.txt'}, '.txt'); % read the plane wave lists back from disk obj.planeWaveList_IS = obj.readPlaneWaveList(planewave_files_IS); obj.planeWaveList_RT = obj.readPlaneWaveList(planewave_files_RT); else obj.planeWaveList_IS = []; obj.planeWaveList_RT = []; end end %------------------------------------------------------------------ function deleteResultsInRavenFolder(obj) if exist(fullfile(obj.pathResults, obj.projectTag), 'dir') rmdir(fullfile(obj.pathResults, obj.projectTag), 's'); end end %------------------------------------------------------------------ function restoreReceiverGrid(obj) numReceivers = numel(obj.receiverPositions) / 3; %xyz if numReceivers < 1 error('No receiver position data present.'); end x = obj.receiverPositions(:, 1); y = obj.receiverPositions(:, 2); z = obj.receiverPositions(:, 3); [x, new_index_order] = sort(x); y = y(new_index_order); z = z(new_index_order); uniform_x = []; uniform_y = []; uniform_z = []; row_size = 0; i = 1; while i <= numel(x) uniform_row_indices = (x==x(i)); row_x = x(uniform_row_indices); row_y = y(uniform_row_indices); row_z = z(uniform_row_indices); [row_z, new_row_order] = sort(row_z); row_x = row_x(new_row_order); row_y = row_y(new_row_order); uniform_x = [uniform_x row_x]; uniform_y = [uniform_y row_y]; uniform_z = [uniform_z row_z]; row_size = numel(row_x); i = i + numel(row_x); end if row_size > 0 if rem(numel(x), row_size) == 0 col_size = numel(x) / row_size; obj.uniformReceiverGridX = reshape(uniform_x, row_size, col_size); obj.uniformReceiverGridY = reshape(uniform_y, row_size, col_size); obj.uniformReceiverGridZ = reshape(uniform_z, row_size, col_size); end end end end % private methods %---------------------- STATIC METHODS -------------------------------% methods(Static) %------------------------------------------------------------------ function rvn_pos = pSU2RVN(su_pos) if size(su_pos, 2) ~= 3 error('Input positions must be given in Nx3 matrix.'); end rvn_pos = su_pos; rvn_pos(:, 1) = sign(RavenProject.COORD_TRAFO_SKETCHUP2RAVEN(1)) .* su_pos(:, abs(RavenProject.COORD_TRAFO_SKETCHUP2RAVEN(1))); rvn_pos(:, 2) = sign(RavenProject.COORD_TRAFO_SKETCHUP2RAVEN(2)) .* su_pos(:, abs(RavenProject.COORD_TRAFO_SKETCHUP2RAVEN(2))); rvn_pos(:, 3) = sign(RavenProject.COORD_TRAFO_SKETCHUP2RAVEN(3)) .* su_pos(:, abs(RavenProject.COORD_TRAFO_SKETCHUP2RAVEN(3))); end %------------------------------------------------------------------ function su_pos = pRVN2SU(rvn_pos) if size(rvn_pos, 2) ~= 3 error('Input positions must be given in Nx3 matrix.'); end su_pos = rvn_pos; su_pos(:, 1) = sign(RavenProject.COORD_TRAFO_RAVEN2SKETCHUP(1)) .* rvn_pos(:, abs(RavenProject.COORD_TRAFO_RAVEN2SKETCHUP(1))); su_pos(:, 2) = sign(RavenProject.COORD_TRAFO_RAVEN2SKETCHUP(2)) .* rvn_pos(:, abs(RavenProject.COORD_TRAFO_RAVEN2SKETCHUP(2))); su_pos(:, 3) = sign(RavenProject.COORD_TRAFO_RAVEN2SKETCHUP(3)) .* rvn_pos(:, abs(RavenProject.COORD_TRAFO_RAVEN2SKETCHUP(3))); end %------------------------------------------------------------------ function outputFilename = convolve(soundFilename, IRitaAudioOrFile, outputFilename, overwrite, amplification_dB) %outputFilename = convolve(soundFilename, IRitaAudio, outputFilename, overwrite, amplification_dB) if ~exist(soundFilename, 'file') error('Sound file not found.'); end % check if convolution program is available in current directory if ~exist('FIRFilter.exe', 'file') error('Convolution engine (''FIRFilter.exe'') not found in current directory.'); end if exist(outputFilename, 'file') % ask user if file should be overwritten if nargin < 4 answ = questdlg('Do you want to overwrite the target file?', ... 'File already exists. Overwrite?', ... 'Overwrite', 'Cancel', 'Overwrite'); switch answ, case 'Overwrite', overwrite = 1; case 'Cancel', outputFilename = ''; return; end % switch end if ~overwrite warning('File already exists and was not allowed to be overwritten.'); return; end end if isa(IRitaAudioOrFile, 'itaAudio') tmpFileName = 'RavenProjectTempImpulseResponse1234567890.wav'; ita_write(IRitaAudioOrFile, tmpFileName, 'overwrite', 'nbits', 32); elseif ischar(IRitaAudioOrFile) tmpFileName = IRitaAudioOrFile; else error('Couldn''t read impulse response.'); end executeCommand = ['FIRFilter.exe ']; if exist('amplification_dB', 'var') && isnumeric(amplification_dB) && (amplification_dB ~= 0) executeCommand = [executeCommand '-g ' num2str(amplification_dB) ' ']; end executeCommand = [executeCommand '"' soundFilename '" "' tmpFileName '" "' outputFilename '"']; disp(executeCommand); system(executeCommand); if isa(IRitaAudioOrFile, 'itaAudio') delete(tmpFileName); end end %------------------------------------------------------------------ function outfile = mix(outfile, varargin) % outfile = mix(outfile, infile1, [infile2], [infile3], ...) % outfile = mix(outfile, {infile1, [infile2], [infile3], ...}) if nargin < 2 error('At least two parameters expected.'); end % check if mixing program is available in current directory if ~exist('CLmix.exe', 'file') error('Mixing engine (''CLmix.exe'') not found in current directory.'); end % make command executeCommand = ['CLmix.exe "' outfile]; % append input file names if iscell(varargin{1}) for iFile = 1 : numel(varargin{1}) executeCommand = [executeCommand '" "' varargin{1}{iFile}]; end else for iFile = 1 : nargin-1 executeCommand = [executeCommand '" "' varargin{iFile}]; end end executeCommand = [executeCommand '"']; % execute command disp(executeCommand); system(executeCommand); end %------------------------------------------------------------------ function lspSignals = ambisonicsDecoding(BformatItaAudio, decodingMatrix, nfcFilter) % get parameters nSmpls = BformatItaAudio.nSamples; nLS = size(decodingMatrix, 2); % pre-allocate memory lspSignals = itaAudio; lspSignals.timeData = zeros(nSmpls, nLS); % partition data to avoid memory problems blockSize = 1000; numCycles = ceil(nSmpls / blockSize); startSample = 1; endSample = startSample + blockSize; % signal too short? (nSmpls < blockSize) if numCycles < 2 endSample = nSmpls; end % NFC if nargin > 2 % 1:1 channel by channel convolution BformatItaAudio = ita_convolve(BformatItaAudio, nfcFilter); % TODO: undo time shift disp('NFC filters did apply a time shift which should be reversed. TODO.'); end % DECODING h = waitbar(0, 'Decoding...'); for i = 1 : numCycles lspSignals.timeData(startSample : endSample, :) = real(BformatItaAudio.timeData(startSample : endSample, 1:size(decodingMatrix, 1)) * decodingMatrix); startSample = startSample + blockSize; endSample = endSample + blockSize; if endSample > size(BformatItaAudio.timeData,1) endSample = size(BformatItaAudio.timeData,1); end if ishandle(h) waitbar(i / numCycles, h); end end delete(h); end %------------------------------------------------------------------ function energySpectrum = powerSpectrum2energySpectrum(powerSpectrum) RF_DEF_thirdOctaveBoundary = [ ... 22.4, ... 27.9, 35.1, 44.2, ... 55.6, 70.1, 88.4, ... 111.4, 140.3, 176.8, ... 222.8, 280.7, 353.6, ... 445.4, 561.2, 707.1, ... 890.9, 1122.5, 1414.2, ... 1781.8, 2244.9, 2828.4, ... 3563.6, 4489.9, 5656.9, ... 7127.3, 8979.8, 11313.8, ... 14255.0, 17960.0, 22627.6 ]; RF_DEF_octaveBoundary = [ ... 44.2, 88.4, 176.8, ... 353.6, 707.1, 1414.2, ... 2828.4, 5656.9, 11313.8, ... 22627.6 ]; numFrequencyBands = size(powerSpectrum, 2); numTimeSlots = size(powerSpectrum, 1); energySpectrum = zeros(numTimeSlots, numFrequencyBands); % OCTAVE RESOLUTION if (numFrequencyBands == 10) % get total bandwidth in Hz of all octave bands bandwidthTotal = RF_DEF_octaveBoundary(10); % RF_DEF_octaveBoundary array defines the upper band limits (see RD_Defines.h) for frequencystep = 1 : numFrequencyBands % define the upper and lower band borders of the current band upperBandLimit = RF_DEF_octaveBoundary(frequencystep); if (frequencystep == 1) lowerBandLimit = 0; else lowerBandLimit = RF_DEF_octaveBoundary(frequencystep - 1); end % get bandwidth in Hz of the current band bandwidthCurrentBand = upperBandLimit - lowerBandLimit; % multiply the power level at the center frequency with the band's width to get the energy per band energySpectrum(:, frequencystep) = (bandwidthCurrentBand / bandwidthTotal) * powerSpectrum(:, frequencystep); end % THIRD-OCTAVE RESOLUTION elseif (numFrequencyBands == 31) % get total bandwidth in Hz of all third-octave bands bandwidthTotal = RF_DEF_thirdOctaveBoundary(31); % RF_DEF_octaveBoundary array defines the upper band limits (see RD_Defines.h) for frequencystep = 1 : numFrequencyBands % define the upper and lower band borders of the current band upperBandLimit = RF_DEF_thirdOctaveBoundary(frequencystep); if (frequencystep == 1) lowerBandLimit = 0; else lowerBandLimit = RF_DEF_thirdOctaveBoundary(frequencystep - 1); end % get bandwidth in Hz of the current band bandwidthCurrentBand = upperBandLimit - lowerBandLimit; % multiply the power level at the center frequency with the band's width to get the energy per band energySpectrum(:, frequencystep) = (bandwidthCurrentBand / bandwidthTotal) * powerSpectrum(:, frequencystep); end % UNKNOWN RESOLUTION - ERROR else disp('Error in RavenProject.powerSpectrum2energySpectrum: Unknown frequency resolution.'); energySpectrum = []; end end %------------------------------------------------------------------ function thisReverbTime = getT30_fromBinauralImpulseResponse(brir) if ~iscell(brir) brir = {brir}; end thisReverbTime = cell(numel(brir), 1); for i = 1 : numel(brir) if ~isempty(brir{i}) disp('third-octave filtering'); brir{i} = ita_mpb_filter(brir{i}, '3-oct'); disp('schroeder decay curve and reverberation time'); thisReverbTime{i} = ita_roomacoustics_reverberation_time(brir{i}, 'edc_type','normal'); disp('average left and right ear''s reverberation time'); if numel(thisReverbTime{i}) > 1 thisReverbTime{i} = mean(thisReverbTime{i}); end thisReverbTime{i} = mean([thisReverbTime{i}.freqData(1:2:60) thisReverbTime{i}.freqData(2:2:60)], 2); else thisReverbTime{i} = []; end end if numel(thisReverbTime) == 1 thisReverbTime = thisReverbTime{1}; end end %------------------------------------------------------------------ function value_out = averageAfterDIN(values_in, afterDIN) if nargin < 2 afterDIN = 1; % default -> average after DIN, otherwise just return mean end if ~iscell(values_in) values_in = {values_in}; end for iCell = 1 : numel(values_in) % averaging a la DIN (average of 500Hz and 1kHz for % octaves or 400, 500, 630, 800, 1000, 1250 Hz for % third-octave resolution if (numel(values_in{iCell}) == 10) && afterDIN value_out{iCell} = (values_in{iCell}(5) + values_in{iCell}(6)) / 2; % 500 and 1000 Hz elseif (numel(values_in{iCell}) == 31) && afterDIN value_out{iCell} = ... (values_in{iCell}(14) + ... % 400 Hz values_in{iCell}(15) + ... % 500 Hz values_in{iCell}(16) + ... % 630 Hz values_in{iCell}(17) + ... % 800 Hz values_in{iCell}(18) + ... % 1000 Hz values_in{iCell}(19)) / 6; % 1250 Hz else % neither octaves nor thirds, average over all bands value_out{iCell} = mean(values_in{iCell}); end end if numel(value_out) == 1 value_out = value_out{1}; end end %------------------------------------------------------------------ function value_out = averageAfterDINLateral(values_in, afterDIN) if nargin < 2 afterDIN = 1; % default -> average after DIN, otherwise just return mean end if ~iscell(values_in) values_in = {values_in}; end for iCell = 1 : numel(values_in) % averaging a la DIN (average of 125Hz ro 1kHz for % octaves or 100 to 1250 Hz for third-octave resolution if (numel(values_in{iCell}) == 10) && afterDIN value_out{iCell} = ... (values_in{iCell}(3) + ... % 125 Hz values_in{iCell}(4) + ... % 250 Hz values_in{iCell}(5) + ... % 500 Hz values_in{iCell}(6)) / 4; % 1000 Hz elseif (numel(values_in{iCell}) == 31) && afterDIN value_out{iCell} = ... (values_in{iCell}(8) + ... % 100 Hz values_in{iCell}(9) + ... % 125 Hz values_in{iCell}(10) + ... % 160 Hz values_in{iCell}(11) + ... % 200 Hz values_in{iCell}(12) + ... % 250 Hz values_in{iCell}(13) + ... % 315 Hz values_in{iCell}(14) + ... % 400 Hz values_in{iCell}(15) + ... % 500 Hz values_in{iCell}(16) + ... % 630 Hz values_in{iCell}(17) + ... % 800 Hz values_in{iCell}(18) + ... % 1000 Hz values_in{iCell}(19)) / 12; % 1250 Hz else % neither octaves nor thirds, average over all bands value_out{iCell} = mean(values_in{iCell}); end end if numel(value_out) == 1 value_out = value_out{1}; end end %------------------------------------------------------------------ function value_out = averageOverReceivers(rec_data) if ~iscell(rec_data) rec_data = {rec_data}; end summe = []; validReceiverCount = 0; for iRec = 1 : numel(rec_data) if ~isempty(rec_data{iRec}) if isempty(summe) summe = rec_data{iRec}; else summe = summe + rec_data{iRec}; end validReceiverCount = validReceiverCount + 1; end end if validReceiverCount > 0 % average over all receivers value_out = summe / validReceiverCount; else value_out = []; end end %------------------------------------------------------------------ function histo = applyTemporalLowPassToLatePartOfHistogram(histo, min_block_size, max_block_size) numTimeSteps = size(histo, 1); if nargin < 2 min_block_size = 0.9; % just slightly below 1 to avoid 2-slot-windows at the beginning (1/2 = 0.5 -> rounded to 1) end if nargin < 2 max_block_size = numTimeSteps * 0.1; % at the end moving average windows will have have the size of 10% of the impulse response end startTimeStep = ceil(min_block_size / 2) + 1; endTimeStep = numTimeSteps - ceil(max_block_size / 2); currentWindowLength = min_block_size; stepFactor = (max_block_size/min_block_size)^(1/numTimeSteps); for timestep = startTimeStep : endTimeStep windowBegin = round(timestep - currentWindowLength/2); windowEnd = round(timestep + currentWindowLength/2); windowLength = windowEnd - windowBegin + 1; histo(timestep, :) = sum(histo(windowBegin:windowEnd, :), 1) / windowLength; currentWindowLength = currentWindowLength * stepFactor; end for timestep = endTimeStep : numTimeSteps windowBegin = round(timestep - currentWindowLength/2); windowEnd = numTimeSteps; windowLength = windowEnd - windowBegin + 1; histo(timestep, :) = sum(histo(windowBegin:windowEnd, :), 1) / windowLength; end end %------------------------------------------------------------------ function wav_files = scan_output_folder(outputDir, searchTag, extension) % % scan_output_folder(outputDir, searchTag, extension) % % Reads all impulse responses in the given output directory that contain % the seachTag(s) in their filename. SearchTag might be e.g. "BRIR_IS_", % "BRIR_RT_" or "BRIR_Combined_". SearchTag can also be an array of % cells: {'BRIR_IS_', 'portal'}. If no Extension is given, the function % will return all .wav files by default. % if nargin < 3 extension = '.wav'; end if (nargin > 1) && ~isempty(searchTag) if ~iscell(searchTag) searchTag = {searchTag}; end searchTag = lower(searchTag); wav_files = []; all_ir_files = dirrec(outputDir, extension); for i_ir = 1 : numel(all_ir_files) % look for searchTag(s) in current file name tags_found = true; for i_tag = 1 : numel(searchTag) if isempty(strfind(lower(all_ir_files{i_ir}), searchTag{i_tag})) tags_found = false; end end if tags_found wav_files{end+1} = all_ir_files{i_ir}; end end else % no search tags given, output all files wav_files = dirrec(outputDir, extension); end end %------------------------------------------------------------------ function out_string = cat_cell_of_strings(cell_of_strings) if iscell(cell_of_strings) out_string = ''; for i = 1 : numel(cell_of_strings) out_string = [out_string cell_of_strings{i} ',']; end out_string(end) = []; % delete last "," else out_string = cell_of_strings; end end %------------------------------------------------------------------ function out_string = make_proper_string(in_matrix) if iscell(in_matrix) out_string = ''; for i = 1 : numel(in_matrix) out_string = [out_string num2str(in_matrix{i}, '%.3f,')]; end out_string(end) = []; % delete last "," else if numel(in_matrix) > 3 temp = num2str(in_matrix(:,:), '%.3f,'); out_string(1:numel(temp)) = deal(temp'); out_string(end) = []; % delete last ", " elseif numel(in_matrix) > 0 % 1 frame out_string = num2str(in_matrix(:)', '%.3f,'); out_string(end) = []; % delete last ", " else out_string = []; end end end %------------------------------------------------------------------ function out_string = makeNumberedNames(name, xtimes) out_string = sprintf([name '%05d,'], 0 : xtimes-1); out_string(end) = []; % delete last "," end %------------------------------------------------------------------ function out_string = writeXtimes(str, xtimes) [cell_of_strings{1:xtimes, 1}] = deal(str); out_string = RavenProject.cat_cell_of_strings(cell_of_strings); end %------------------------------------------------------------------ function out_string = writeXtimes_num(num, xtimes) str = num2str(num); out_string = RavenProject.writeXtimes(str, xtimes); end end % static methods end % classdef %--------------------------------------------------------------------------