Commit 1ed57954 authored by Lukas Aspöck's avatar Lukas Aspöck

Merge branch 'master' of https://git.rwth-aachen.de/ita/toolbox

parents 385d9751 bd0b3942
......@@ -368,16 +368,29 @@ classdef itaMSTF < itaMSPlaybackRecord
function sweepRate = sweepRate(this,value)
% get the sweep rate of the excitation
%% sweep rate from analytic calculation, only using sweep parameters / PDI
nSamples = ita_nSamples( this.fftDegree );
% MMT: use nSamples-1 here to be conform with sweep calculation
% based on timeVector and chirp function
finalExcitationLength = (nSamples-1)/this.samplingRate - this.stopMargin;
sweepRate(1) = log2(this.finalFreqRange(2)/this.finalFreqRange(1))/finalExcitationLength;
%% sweep rate of analysis of excitation signal
sweepRate(2) = ita_sweep_rate(this.raw_excitation,[2000 this.samplingRate/3]);
sweepRate = nan(2,1);
switch this.type
case 'exp'
% equation is: f(t) = f0*2^(sweeprate*t)
% sweep rate from analytic calculation, only using sweep parameters
nSamples = ita_nSamples( this.fftDegree );
% use nSamples-1 here to be conform with sweep calculation
% based on timeVector and chirp function
finalExcitationLength = (nSamples-1)/this.samplingRate - this.stopMargin;
sweepRate(1) = log2(this.finalFreqRange(2)/this.finalFreqRange(1))/finalExcitationLength;
case 'lin'
% equation is: f(t) = f0*(1 + sweeprate*t)
% sweep rate from analytic calculation, only using sweep parameters
nSamples = ita_nSamples( this.fftDegree );
% use nSamples-1 here to be conform with sweep calculation
% based on timeVector and chirp function
finalExcitationLength = (nSamples-1)/this.samplingRate - this.stopMargin;
sweepRate(1) = (this.finalFreqRange(2)/this.finalFreqRange(1) - 1)/finalExcitationLength;
otherwise
end
sweepRate(2) = ita_sweep_rate(this.raw_excitation,'freqRange',[max(this.finalFreqRange(1),2000) min(this.finalFreqRange(2),this.samplingRate/3)],'type',this.type,'f0',this.finalFreqRange(1));
if exist('value','var')
sweepRate = sweepRate(value);
end
......
......@@ -174,7 +174,7 @@ clipping = false;
if record
% determine clipping limit from NI session information
clippingLimit = Inf;
for iChannel = numel(niSession.Channels)
for iChannel = 1:numel(niSession.Channels)
isInput = ~isempty(strfind(niSession.Channels(iChannel).ID,'ai'));
if isInput
clippingLimit = min(clippingLimit,max(abs(double(niSession.Channels(iChannel).Range))));
......
......@@ -15,7 +15,7 @@ sArgs = struct('pos1_data','anything','bandsperoctave',3,'freqVector',[],'create
%% reference curves
if strcmpi(sArgs.type,'iso') % Reference curve and frequencies according to ISO 717-2
outputStr = 'Ln_w (C_I)';
outputStr = 'Ln_w (C_I, C_{I50-2500})';
roundingFactor = 0.1;
deficiencyLimit = Inf;
if sArgs.bandsperoctave == 1
......@@ -81,7 +81,13 @@ end
%% adaptation term for ISO
if strcmpi(sArgs.type,'iso')
C = round(round((10.*log10(sum(10.^(NISPL./10))))/0.1)*0.1 - 15 - impactInsulationClass);
C = round(round((10.*log10(sum(data.freq2value(100,2500).^2))+93.98)/0.1)*0.1 - 15 - impactInsulationClass);
if min(data.freqVector) <= 50 % if we have data
C_50 = round(round((10.*log10(sum(data.freq2value(50,2500).^2))+93.98)/0.1)*0.1 - 15 - impactInsulationClass);
else
C_50 = nan;
end
C = [C C_50];
else
C = 0;
end
......@@ -100,7 +106,7 @@ if sArgs.createPlot
bar(gca,deficiencies.freqVector,deficiencies.freq,'hist');
[maxDef,maxIdx] = max(deficiencies.freq);
if strcmpi(sArgs.type,'iso')
singleNumberString = [outputStr ' = ' num2str(impactInsulationClass) ' (' num2str(C) ') dB'];
singleNumberString = [outputStr ' = ' num2str(impactInsulationClass) ' (' num2str(C(1)) ',' num2str(C(2)) ') dB'];
else
singleNumberString = [outputStr ' = ' num2str(impactInsulationClass) 'dB'];
end
......
......@@ -207,6 +207,11 @@ classdef (Abstract)itaAbstract3DModelVisualizer < handle
axis(this.mAxes, 'off');
axis(this.mAxes, 'equal');
end
function ClearPlot(this)
%Clears all plot items
this.clearPlotItems();
end
end
......
......@@ -39,7 +39,7 @@ classdef itaAc3dVisualizer < itaAbstract3DModelVisualizer
elseif isa(input, 'itaAc3dModel') && isscalar(input)
obj.mModel = input;
else
error('Input must be either a valid .ac filename or a AC3D object (itaAc3dModel).')
error('Input must be either a valid .ac filename or an AC3D object (itaAc3dModel).')
end
obj.clearPlotItems();
......
......@@ -80,8 +80,8 @@ classdef itaComsolPhysics < itaComsolNode
impedanceNodes = obj.getFeatureNodesByType(obj.mActiveNode, 'Impedance');
end
function impedanceNodeOfBoundary = ImpedanceNodeByBoundaryGroupName(obj, boundaryGroupName)
%Returns the impedance node of the active physics node that is
%connected to given boundary group (= Comsol selection). Returns
%Returns the first impedance node of the active physics node that
%is connected to given boundary group (= Comsol selection). Returns
%an empty object if no impedance is connected to the given selection.
assert(ischar(boundaryGroupName) && isrow(boundaryGroupName), 'Input must be a char row vector')
impedanceNodeOfBoundary = [];
......@@ -138,6 +138,50 @@ classdef itaComsolPhysics < itaComsolNode
end
end
%% Poroacoustics
methods
function poroacousticsNodes = PoroacousticsNodes(obj)
%Returns the poroacoustics nodes of the active physics node.
poroacousticsNodes = obj.getFeatureNodesByType(obj.mActiveNode, 'PoroacousticsModel');
end
function poroacousticsNodeOfDomain = PoroacousticsNodeByVolumeGroupName(obj, volumeGroupName)
%Returns the first poroacoustics node of the active physics node
%that is connected to given volume group (= Comsol domain selection).
%Returns an empty object if no poroacoustics model is connected
%to the given selection.
assert(ischar(volumeGroupName) && isrow(volumeGroupName), 'Input must be a char row vector')
poroacousticsNodeOfDomain = [];
volumeGroup = obj.mModel.selection.VolumeGroup(volumeGroupName);
if isempty(volumeGroup); return; end
poroacousticsNodes = obj.PoroacousticsNodes();
for idxPoroacoustics = 1:numel(poroacousticsNodes)
selectionTag = poroacousticsNodes{idxPoroacoustics}.selection.named;
if strcmp(selectionTag, volumeGroup.tag)
poroacousticsNodeOfDomain = poroacousticsNodes{idxPoroacoustics};
return
end
end
end
function poroacousticsNode = CreatePoroacoustics(obj, poroacousticsTag, selectionTag)
%Creates a poroacoustics model for the active physics node
% Inputs:
% poroacousticsTag Tag for poroacoustics node [char row vector]
% selectionTag Tag of the selection the poroacoustics model is linked to [char row vector]
assert(ischar(poroacousticsTag) && isrow(poroacousticsTag), 'First input must be a char row vector')
assert(ischar(selectionTag) && isrow(selectionTag), 'Second input must be a char row vector')
physics = obj.activeNode;
if ~obj.hasFeatureNode(physics, poroacousticsTag)
physics.create(poroacousticsTag, 'PoroacousticsModel', 3);
end
poroacousticsNode = physics.feature(poroacousticsTag);
poroacousticsNode.selection.named(selectionTag);
end
end
%% Source related
methods
function normalVelocityNode = CreateNormalVelocity(obj, sourceTag, selectionTag, normalVelocityExpression)
......
......@@ -81,8 +81,9 @@ classdef itaComsolResult < handle
end
function [res, metaData] = ByExpression(obj, expression, evalPos)
%For the active study, evaluates the given expression either at
%the mesh nodes or at the given itaCoordinates. Also returns a
%metadata struct that contains information for parametric sweeps.
%the mesh nodes (FEM only) or at the given itaCoordinates. Also
%returns a metadata struct that contains information for#
%parametric sweeps.
% Inputs:
% expression: char row-vector with expression to be evaluated
% evalPos (optional): itaCoordinates with coordinates for evaluation
......@@ -93,15 +94,12 @@ classdef itaComsolResult < handle
% metaData: struct with info on parametric sweep
assert(ischar(expression) && isrow(expression), 'First input must be char row vector with the expression of the result')
%NOTE - PSC:
%Since the Comsol functions mphinterp and mpheval work
%significantly different for the pabe physics, getting pabe
%results is disabled for now. I am waiting for the Comsol
%Support to help me with this.
if contains(expression, 'pabe')
error('Evaluating BEM (pabe) results is not supported at the moment');
end
if nargin == 2
%NOTE - PSC:
%Since the Comsol function mpheval seems to have a bug if
%using the pabe physics, getting pabe results at the mesh
%nodes is disabled.
assert( ~contains(expression, 'pabe'), 'You have to specify coordinates if evaluating a BEM (pabe) result' )
[res, metaData] = obj.getResultAtMeshNodes(expression);
else
if isa(evalPos, 'itaReceiver')
......@@ -116,16 +114,15 @@ classdef itaComsolResult < handle
%% Extract results
methods(Access = private)
function [res, metaData] = getResultAtMeshNodes(obj, expression)
assert( ~contains(expression, 'pabe'), 'Cannot return expression' )
datasetTag = obj.getDirectDatasetTag();
metaData = obj.createMetaDataStruct(datasetTag);
if ~isempty(metaData.parameterValues) %Param sweep
res = itaResult([1 metaData.nSimulations]);
for idxParam = 1:metaData.nSimulations
%Note that 'outersolnum' = 'all' does not work if the
%mesh is parameter dependent since then the number of
%data points differs between parametric solutions.
%NOTE - PSC:
%'outersolnum' = 'all' does not work if the mesh is
%parameter-dependent since then the number of data
%points differs between parametric solutions.
%Thus, we have to get the results one by one.
data = mpheval(obj.mModel.modelNode, expression, 'dataset',datasetTag,'outersolnum',idxParam);
itaCoords = itaCoordinates(data.p.');
......@@ -142,10 +139,24 @@ classdef itaComsolResult < handle
end
end
function [res, metaData] = getResultAtCoords(obj, expression, itaCoords)
datasetTag = obj.getDatasetTag(expression);
freqData = mphinterp(obj.mModel.modelNode, expression, 'coord', itaCoords.cart.', 'dataset', datasetTag,'outersolnum','all');
freqVector = mphglobal(obj.mModel.modelNode,'freq','dataset',datasetTag,'outersolnum','all');
metaData = obj.createMetaDataStruct(datasetTag);
freqDatasetTag = obj.getResultAtCoordsDatasetTag(expression);
directDatasetTag = obj.getDirectDatasetTag();
metaData = obj.createMetaDataStruct(directDatasetTag);
if contains(expression, 'pabe')
disp('**ITA-COMSOL** Starting evaluating BEM result at user-defined coordinates.')
disp(' This might take a while...')
end
%NOTE - PSC:
%Appearently, 'outersolnum' = 'all' does not work for pabe at
%the moment due to a bug. Thus, the solver range is set
%manually using the metaData struct.
freqData = mphinterp(obj.mModel.modelNode, expression, 'coord', itaCoords.cart.', 'dataset', freqDatasetTag,'outersolnum',1:metaData.nSimulations);
freqVector = mphglobal(obj.mModel.modelNode,'freq','dataset',directDatasetTag,'outersolnum',1:metaData.nSimulations);
if contains(expression, 'pabe')
disp('**ITA-COMSOL** Done!')
end
if numel(size(freqData))==3 %Param sweep
res = obj.createParametricItaResult(freqData, freqVector, itaCoords);
else
......@@ -203,12 +214,17 @@ classdef itaComsolResult < handle
%% Getting dataset from study
methods(Access = private)
function datasetTag = getDatasetTag(obj, expression)
function datasetTag = getResultAtCoordsDatasetTag(obj, expression)
%Returns the dataset used to extract results at user-defined
%coordinates
% In case of acpr, datasetTag is the same as for
% getDirectDatasetTag() but in case of pabe, datasetTag
% refers to a grid dataset.
if contains(expression, 'pabe')
datasetTag = obj.getBemDatasetTag();
else
if ~contains(expression, 'acpr')
warning('Expression does not contain physics-tag. Assuming acpr physics.')
warning('Expression does not contain a known physics-tag. Assuming acpr physics.')
end
datasetTag = obj.getDirectDatasetTag();
end
......@@ -236,6 +252,8 @@ classdef itaComsolResult < handle
solTag = '';
if isempty(allSolvers); return; end
if obj.mModel.study.IsParametric && numel(allSolvers)>= 2
%For a parametric sweep, the first solver is the one for a
%single simulation, while the second contains all simulations.
solTag = char(allSolvers(2));
else
solTag = char(allSolvers(1));
......
......@@ -28,6 +28,9 @@ itaComsolModelObj.study.Run(showProgress);
%% --- Read results ---
%% at mesh nodes
%IMPORTANT NOTE:
%Evaluating at mesh nodes only works if using FEM (acpr). It does not work
%if using BEM (pabe).
resultAtMeshNodes = itaComsolModelObj.result.Pressure();
%% at user defined points within mesh
......
......@@ -115,11 +115,14 @@ classdef itaMaterial < itaSimulationInputItem
function this = set.scattering(this, scattering)
if isnumeric(scattering) && isempty(scattering)
this.mScattering = [];
return;
elseif isnumeric(scattering) && isscalar(scattering)
assert(scattering >= 0 && scattering <= 1, 'When applying a constant scattering, the value must be between 0 and 1.')
freqData = scattering*ones( size(this.gaThirdOctavefreqs) );
this.mScattering = itaResult(freqData, this.gaThirdOctavefreqs, 'freq');
else
this.checkDataTypeForFreqData(scattering)
this.mScattering = scattering;
end
this.checkDataTypeForFreqData(scattering)
this.mScattering = scattering;
end
function this = set.rho0Air(this, rho0)
......
......@@ -50,11 +50,24 @@ else
end
if max(size(asData)) > 1
for idx = 1:numel(asData)
result(idx) = ita_extract_dat(asData(idx),varargin(2:end));
varargout{1} = result;
return;
end
for idx = 1:numel(asData)
% MMT: this seems to be the easiest way to make this work
if isnumeric(varargin{2})
optionStr = num2str(varargin{2});
else
optionStr = ['''' varargin{2} ''''];
end
for iVar = 3:numel(varargin)
if isnumeric(varargin{iVar})
optionStr = [optionStr ',' num2str(varargin{iVar})]; %#ok<*AGROW>
else
optionStr = [optionStr ',''' varargin{iVar} ''''];
end
end
eval(['result(idx) = ita_extract_dat(asData(idx),' optionStr ');']);
end
varargout{1} = result;
return;
end
......
......@@ -41,7 +41,7 @@ end
%% Initialization and Input Parsing
sArgs = struct('pos1_a','itaAudio');
[Data, sArgs] = ita_parse_arguments(sArgs,varargin(1)); %#ok<NASGU>
[Data, sArgs] = ita_parse_arguments(sArgs,varargin(1)); %#ok<ASGLU>
if nargin == 1
NewSamplingRate = 44100;
......@@ -50,6 +50,16 @@ else
NewSamplingRate = varargin{2};
end
if numel(Data) > 1
ita_verbose_info('Calling for all instances.',1)
result = itaAudio(size(Data));
for idx = 1:numel(Data)
result(idx) = ita_resample(Data(idx),NewSamplingRate);
end
varargout{1} = result;
return
end
OldSamplingRate = Data.samplingRate;
if OldSamplingRate == NewSamplingRate
ita_verbose_info('ITA_RESAMPLE:The Sampling Rate is already fine.',1);
......
......@@ -9,47 +9,46 @@ function varargout = ita_sweep_rate(varargin)
% Note: A too short sweep or a too small frequency range may be problematic.
% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved.
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder.
% This file is part of the ITA-Toolbox. Some rights reserved.
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder.
% </ITA-Toolbox>
% Author: Alexandre Bleus -- alexandre.bleus@akustik.rwth-aachen.de
% Created: 12-March-2010
%% Get ITA Toolbox preferences
thisFuncStr = [upper(mfilename) ':'];
%% Initialisation
narginchk(1,2);
if nargin == 1
varargin{2} = [1000 4000];
elseif nargin~=2
error([thisFuncStr 'Two input arguments are required.'])
end
if ~isa(varargin{1},'itaAudio')
error([thisFuncStr 'Oh Lord! An audio object (Exponential Sweep) is required.'])
else
inSweep = varargin{1};
freq_vec=varargin{2};
end
sArgs = struct('pos1_inSweep','itaAudio', 'freqRange',[1000 4000],'type','exp','f0',[]);
[inSweep,sArgs] = ita_parse_arguments(sArgs,varargin);
%% Calculate the sweep rate for some points
freqjump = round((freq_vec(2)-freq_vec(1))/10);
freq_gdel = freq_vec(1) + [(1:8).' (2:9).'].*freqjump;
freqjump = round((sArgs.freqRange(2)-sArgs.freqRange(1))/10);
freq_gdel = sArgs.freqRange(1) + [(1:8).' (2:9).'].*freqjump;
gdel = ita_groupdelay_ita(inSweep);
sweep_rate = log(freq_gdel(:,2)./freq_gdel(:,1)) ./ log(2) ./ (gdel.freq2value(freq_gdel(:,2)) - gdel.freq2value(freq_gdel(:,1)));
sweep_rate(sweep_rate < 0) = sweep_rate(end);
switch sArgs.type
case 'exp'
sweep_rate = log2(freq_gdel(:,2)./freq_gdel(:,1)) ./ (gdel.freq2value(freq_gdel(:,2)) - gdel.freq2value(freq_gdel(:,1)));
case 'lin'
sweep_rate = (freq_gdel(:,2) - freq_gdel(:,1)) ./ (gdel.freq2value(freq_gdel(:,2)) - gdel.freq2value(freq_gdel(:,1)));
if isempty(sArgs.f0)
error('Need f0 for linear sweep rate');
else
sweep_rate = sweep_rate./sArgs.f0;
end
otherwise
error('Type can only be exp or lin');
end
sweep_rate(sweep_rate < 0) = nan;
%% Verify that the sweep rates are similar
meanRate = round(10*mean(sweep_rate));
meanRate = round(10*mean(sweep_rate,'omitnan'));
if any(abs(round(sweep_rate*10) - meanRate) > 5)
ita_verbose_info('Oh Lord! Are you sure that you are using an Exponential Sweep?',0);
end
%% Set Output
inSweep.userData.sweep_rate = mean(sweep_rate);
inSweep.userData.sweep_rate = mean(sweep_rate,'omitnan');
varargout(1) = {inSweep.userData.sweep_rate};
end
\ No newline at end of file
This diff is collapsed.
......@@ -451,6 +451,12 @@ if ~iscell(plot_type)
plot_type = {plot_type};
end
if numel(plot_type) == 1 && Nplots > 1
plot_type = repmat(plot_type,[Nplots 1]);
end
errup = cell(Nplots,1);
errdown = cell(Nplots,1);
for i = 1:Nplots
if strcmpi(plot_type{i},'bar')
if size(x_data{i},1) > 1
......@@ -461,17 +467,15 @@ for i = 1:Nplots
end
end
errup = cell(Nplots,1);
errdown = cell(Nplots,1);
if strcmpi(plot_type{i},'errorbar')
errup = get(chdr(i),'UData');
errdown = get(chdr(i),'LData');
if ~iscell(errup)
errup = {errup};
end
if ~iscell(errdown)
errdown = {errdown};
end
errup{i} = get(chdr(i),'UData');
errdown{i} = get(chdr(i),'LData');
% if ~iscell(errup)
% errup = {errup};
% end
% if ~iscell(errdown)
% errdown = {errdown};
% end
end
end
......
function varargout = ita_plot_variation_area(input, sArgs)
%ITA_PLOT_VARIATION - This is a private helper function for
%ita_plot_variation that creates the area plots
% Do not use this function directly! Returns the legend string for the
% created area plot.
%
% Syntax:
% audioObjOut = ita_plot_variation_area(audioObjIn, sArgs)
%
% sArgs is the argument struct given by ita_plot_variation.
%
% Example:
% legendString = ita_plot_variation_area(audioObjIn)
%
% See also:
% ita_plot_variation
%
% Reference page in Help browser
% <a href="matlab:doc ita_plot_variation_area">doc ita_plot_variation_area</a>
% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved.
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder.
% </ITA-Toolbox>
% Author: Philipp Schfer -- Email: psc@akustik.rwth-aachen.de
% Created: 22-Oct-2019
%% Initialization
domain = sArgs.domain;
xScaleLog = sArgs.logXscale;
axh = sArgs.axes_handle;
% area color
if numel(sArgs.areaColor) == 3
areaColor = sArgs.areaColor;
else
error('wrong area color format')
end
% edge color
if ischar(sArgs.edgeColor) && strcmp(sArgs.edgeColor, 'same')
edgeColor = areaColor;
else
edgeColor = sArgs.edgeColor;
end
%% x-values
xVec = input.([domain 'Vector']);
%% Raw y-data
if strcmpi(domain, 'freq') % erstmal so, bis sich da jemand was schlaues einfallen laesst
if input.allowDBPlot
inputData = input.freqData_dB;
else
inputData = input.freqData;
end
else
inputData = input.([domain 'Data']);
end
%% Remove x=0
if xScaleLog % there is no zero in log scale => delete zeros
idx2del = xVec == 0;
xVec(idx2del) = [];
inputData(idx2del,:) = [];
end
%% Calculate y-values
dataContainsNan = any(isnan(inputData(:)));
switch lower(sArgs.areaMethod)
case 'std'
if dataContainsNan
error('Data contains NaNs. Use option ita_plot_variation(..., ''areaMethod'', ''nanStd'')')
end
yMeanData = mean(inputData, 2);
yStdData = std(inputData, 1, 2);
yAreaData = [yMeanData-yStdData yMeanData+yStdData];
legendStr = 'std';
case 'nanstd'
yMeanData = nanmean(inputData, 2);
yStdData = nanstd(inputData, 1, 2);
yAreaData = [yMeanData-yStdData yMeanData+yStdData];
legendStr = 'std (without NaNs)';
case 'percentile'
percentileValues = sArgs.percValues;
assert(isnumeric(percentileValues) && numel(percentileValues) == 2, 'Percentile values must be given as 2-element vector.')
percentileValues = sort(percentileValues);
yAreaData = prctile(inputData, percentileValues, 2);
legendStr = sprintf('percentiles [%d %d]%%', percentileValues(1), percentileValues(2));
case 'minmax'
yAreaData = [min(inputData,[],2) max(inputData,[],2)];
legendStr = 'min / max';
case 'directinput'
yAreaData = inputData;
legendStr = '';
otherwise
error('unknown areaMethod')
end
%% Remove unwanted data
if sArgs.justTakeFiniteValues
idx2del = any(~isfinite(yAreaData),2);
xVec(idx2del) = [];
yAreaData(idx2del,:) = [];
end
%% plot
if size(xVec,1) ~= size(yAreaData,1)
error('Number samples in x- and y-data does not match')
end
patch(axh, [xVec(:); xVec(end:-1:1)], [yAreaData(:,1); flipud(yAreaData(:,2)) ],...
areaColor, 'FaceAlpha', sArgs.faceTransparency,...
'EdgeColor', edgeColor, 'EdgeAlpha', sArgs.edgeTransparency);
%% Output
if nargout
varargout{1} = legendStr;
end
%end function