diff --git a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/interp.m b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/interp.m index 0f9fc79c0c6a50bcb9874cc127e902d5df2cf5f3..5054c6bf2c07503f117bb99eefe9ec2dd1b9b12f 100644 --- a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/interp.m +++ b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/interp.m @@ -106,7 +106,7 @@ ear_d = [-0.07 0.07]; W = sparse(diag(w)); % diagonal matrix containing weights D = sparse(diag(dweights_rep)); % decomposition order-dependent Tikhonov regularization -Y = ita_sph_base(this.dirCoord,Nmax,'orthonormal',false); % calculate real-valued SHs using the measurement grid +Y = ita_sph_base(this.dirCoord,Nmax,'real'); % calculate real-valued SHs using the measurement grid %% Calculate HRTF data for field points if Nmax > 25 @@ -131,7 +131,7 @@ for ear=1:2 % % calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization % % freqData_temp = data.freqData; -% Y = ita_sph_base(data.channelCoordinates,Nmax,'orthonormal',false); +% Y = ita_sph_base(data.channelCoordinates,Nmax,'real'); a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.'; %a0 = (Y.'*W*Y + epsilon*D) \ Y.' * @@ -142,10 +142,10 @@ for ear=1:2 % calculate range-extrapolated HRTFs a1 = a0 .* hankel_rep.'; - Yest = ita_sph_base(field,N,'orthonormal',false); % use real-valued SH's + Yest = ita_sph_base(field,N,'real'); % use real-valued SH's hrtf_arbi(:,ear:2:end) = (Yest*a1).'; % interpolated + range-extrapolated HRTFs else - Yest = ita_sph_base(field,Nmax,'orthonormal',false); % use real-valued SH's + Yest = ita_sph_base(field,Nmax,'real'); % use real-valued SH's hrtf_arbi(:,ear:2:end) = (Yest*a0).'; % interpolated HRTFs end end diff --git a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m index 7e3c65ea91089dc6f57d6d9194ace9b78289ba6a..6ab573b245135156bf6e3478422983707da35c69 100644 --- a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m +++ b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/itaHRTF.m @@ -1110,7 +1110,7 @@ classdef itaHRTF < itaAudio W = diag(vWeights); % diagonal matrix containing weights D = diag(regweights_rep); % decomposition order-dependent Tikhonov regularization - Y = ita_sph_base(this.dirCoord,Nmeas,'orthonormal',false); % calculate real-valued SHs using the measurement grid (high SH-order) + Y = ita_sph_base(this.dirCoord,Nmeas,'real'); % calculate real-valued SHs using the measurement grid (high SH-order) %% Calculate spatially smoothed HRTF data set hrtf_smoo_wo_ITD = zeros(this.nBins,2*this.dirCoord.nPoints); % init.: columns: LRLRLR... @@ -1302,11 +1302,11 @@ classdef itaHRTF < itaAudio earSidePlot = sArgs.earSide; if numel(phiC_deg)>1, xData = phiC_deg; - strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '°']; + strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '�']; strXlabel = '\phi in Degree'; else xData = thetaC_deg; - strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '°']; + strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '�']; strXlabel = '\theta in Degree'; end @@ -1394,11 +1394,11 @@ classdef itaHRTF < itaAudio earSidePlot = sArgs.earSide; if numel(phiC_deg)>1, xData = phiC_deg; - strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '°']; + strTitle =[ earSidePlot ' ear, \theta = ' num2str(round(thetaC_deg)) '�']; strXlabel = '\phi in Degree'; else xData = thetaC_deg; - strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '°']; + strTitle =[earSidePlot ' ear, \phi = ' num2str(round(phiC_deg)) '�']; strXlabel = '\theta in Degree'; end diff --git a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/private/dfitaHRIRDAFFDataFunc.m b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/private/dfitaHRIRDAFFDataFunc.m index d2fffb293784038c0757dba42642576f5a635025..80031655366ddd53d8c06e793343365debcb645c 100644 --- a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/private/dfitaHRIRDAFFDataFunc.m +++ b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/private/dfitaHRIRDAFFDataFunc.m @@ -6,7 +6,7 @@ function [ data, samplerate, metadata ] = dfitaHRIRDAFFDataFunc( alpha, beta, it % DAFF requires data alignment by multiple of 4 nResidual = mod( hrtf.nSamples, 4 ); - data = [ hrtf.timeData', zeros( hrtf.nChannels, 4 - nResidual ) ]; + data = [ hrtf.timeData', zeros( hrtf.nChannels, mod(4 - nResidual,4) ) ]; metadata = []; end diff --git a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/writeDAFFFile.m b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/writeDAFFFile.m index 4c393d7b91278c110e09220f4faf48d13c26faf8..1f71b73febba30b5c6980911c8ea49ae85725d49 100644 --- a/applications/Binaural-HRTF/HRTF_class/@itaHRTF/writeDAFFFile.m +++ b/applications/Binaural-HRTF/HRTF_class/@itaHRTF/writeDAFFFile.m @@ -44,26 +44,33 @@ end theta_start_deg = rad2deg( min( this.channelCoordinates.theta ) ); theta_end_deg = rad2deg( max( this.channelCoordinates.theta ) ); -theta_num_elements = size( unique( this.channelCoordinates.theta ), 1 ); +theta_num_elements = size( uniquetol( this.channelCoordinates.theta ), 1 ); -phi_start_deg = rad2deg( min( mod( this.channelCoordinates.phi, 2*pi ) ) ); -phi_end_deg = rad2deg( max( mod( this.channelCoordinates.phi, 2*pi ) ) ); -phi_num_elements = size( unique( this.channelCoordinates.phi ), 1 ); +phi_start_deg = rad2deg( min( mod( this.channelCoordinates.phi, 2 * pi ) ) ); +phi_end_deg = rad2deg( max( mod( this.channelCoordinates.phi, 2 * pi ) ) ); +phi_num_elements = size( uniquetol( this.channelCoordinates.phi ), 1 ); assert( phi_num_elements ~= 0 ); alphares = ( phi_end_deg - phi_start_deg ) / phi_num_elements; % phi end does not cover entire circle in this case alphares_full_circle = ( phi_end_deg - phi_start_deg ) / ( phi_num_elements - 1 ); % phi end does not cover entire circle in this case if phi_end_deg + alphares_full_circle >= 360.0 - alpharange = [ phi_start_deg ( phi_end_deg + alphares_full_circle ) ]; % Account for full circle + alpharange = [ phi_start_deg 360 ]; % Account for full circle and force end of range to 360 deg alphares = alphares_full_circle; else alpharange = [ phi_start_deg phi_end_deg ]; end +assert( alpharange( 1 ) >= 0.0 ) +assert( alpharange( 2 ) <= 360.0 ) + assert( theta_num_elements ~= 0 ); betares = ( theta_end_deg - theta_start_deg ) / ( theta_num_elements - 1 ); % phi end does not cover entire circle betarange = 180 - [ theta_start_deg theta_end_deg ]; % Flip poles (DAFF starts at south pole) +assert( betarange( 2 ) >= 0.0 ) +assert( betarange( 1 ) <= 180.0 ) + + %% Assemble metadata metadata = daffv17_add_metadata( metadata, 'Generation script', 'String', 'writeDAFFFile.m' ); diff --git a/applications/Binaural-HRTF/HpTF_class/itaHpTF_Audio.m b/applications/Binaural-HRTF/HpTF_class/itaHpTF_Audio.m index 17056ee71d74f9b6a1708a52b73c0a0d7c9dd55a..909333ca8233a9191667180c1eb6995961d1752a 100644 --- a/applications/Binaural-HRTF/HpTF_class/itaHpTF_Audio.m +++ b/applications/Binaural-HRTF/HpTF_class/itaHpTF_Audio.m @@ -14,6 +14,7 @@ classdef itaHpTF_Audio < itaHpTF % Audio Engineering Society Convention 130, May 2011) % normalized (Normalization of HpTF) % smoothing (Smoothing bandwidth has to be defined in fractions of octaves, see also ita_smooth) + % mic (measured microphone transfer funcions - not inverted: length does not matter; will be adapted) % % itaHpTF_Audio Methods (private): % HP_equalization If this.mic is not set, it will be unused @@ -30,7 +31,9 @@ classdef itaHpTF_Audio < itaHpTF m_fUpper = 18000; mMethod = 'mSTD'; mNormalized = true; + mGainComp = true; mSmoothing = 1/6; + mSelectMeas = 1:8; % select the channels which should be used for the HpTF end properties(Dependent = true, Hidden = false) @@ -38,8 +41,10 @@ classdef itaHpTF_Audio < itaHpTF fLower = 100; fUpper = 18000; method = 'mSTD'; - normalized = true; + normalized = true; % normalization of the signal + gainComp = true; % microphone compensation smoothing = 1/6; + selectMeas = 1:8; % select specific measurements end properties (Dependent = true, SetAccess = private) @@ -122,6 +127,13 @@ classdef itaHpTF_Audio < itaHpTF function normalized = get.normalized(this) normalized = this.mNormalized; end + function comp = get.gainComp(this) + comp = this.mGainComp; end + + function sel = get.selectMeas(this) + sel = this.mSelectMeas; end + + function smoothing = get.smoothing(this) smoothing = this.mSmoothing; end %% ................................................................ @@ -151,6 +163,16 @@ classdef itaHpTF_Audio < itaHpTF this = HP_equalization(this); end + function this = set.selectMeas(this,ch) + this.mSelectMeas = ch; + this = HP_equalization(this); + end + + function this = set.gainComp(this,comp) + this.mGainComp = comp; + this = HP_equalization(this); + end + function this = set.smoothing(this,smoothing) this.mSmoothing = smoothing; this = HP_equalization(this); @@ -162,11 +184,12 @@ classdef itaHpTF_Audio < itaHpTF % SET PRIVATE %.................................................................. function this = set.init(this,var) - this.nameHP = var.nameHP; + this.nameHP = var.nameHP; this.nameMic = var.nameMic; this.nameSubj = var.nameSubj; this.repeat = var.repeat; this.mTF.time = var.time; + this.selectMeas = 1:var.repeat; this = HP_equalization(this); end @@ -175,11 +198,8 @@ classdef itaHpTF_Audio < itaHpTF % Audio Engineering Society Convention 130, May 2011 tWin = this.TF.trackLength; % crop HPTF - if this.mic.dimensions == 2 % mic min phase + extension - minMic = ita_minimumphase(this.mic); - [mic, TF] = ita_extend_dat(minMic,this.TF); - else TF = this.TF; - end + measSel = sort([2*this.selectMeas-1 2*this.selectMeas]); + TF = this.TF.ch(measSel); %init thisEQ = this; @@ -200,14 +220,19 @@ classdef itaHpTF_Audio < itaHpTF otherwise error('Unknown type'); end - if this.mic.dimensions == 2 - Rec = ita_multiply_spk(Rec,mic); + + if this.mic.dimensions == 2 % Compensation of mic + minMic = ita_minimumphase(this.mic); + RecM = ita_convolve(Rec,minMic.ch(cdx)); + RecM.fftDegree = TF.fftDegree; + else + RecM = Rec; end - %% Short Filter with no correction for low freqs and minimun phase - aux = max(abs(Rec.freqData),[],2); + %% Short Filter with no correction for low freqs and minimum phase + aux = max(abs(RecM.freqData),[],2); % find first maximum and truncate low freq correction at this point - idxDF = Rec.freq2index([this.fLower this.fLower*1.5 ]); + idxDF = RecM.freq2index([this.fLower this.fLower*1.5 ]); d_aux = diff(aux(idxDF(1):idxDF(2))); idx = find(diff(sign(d_aux)) ~= 0,1,'first'); % Bruno style... aux(1:idxDF(1)+idx+1) = aux(idxDF(1)+idx+2); @@ -227,10 +252,18 @@ classdef itaHpTF_Audio < itaHpTF this_min = ita_minimumphase(ita_time_shift(HpTF,tWin/2)); this_win = ita_time_window(this_min,[tWin*0.99,tWin],'time','crop'); - if this.normalized, - this_norm = ita_normalize_dat(this_win); + if this.gainComp % compensate gains from mics (equalize) + idxDF = this_win.freq2index(this.fLower); + chGain = 20*log10(abs(this_win.freqData(idxDF,:))); + this_comp = this_win; + this_comp.freqData = bsxfun(@times,this_win.freqData,10.^(-chGain/20)); + else, this_comp = this_win; + end + + if this.normalized + this_norm = ita_normalize_dat(this_comp); thisEQ.timeData = this_norm.timeData; - else thisEQ.timeData = this_win.timeData; + else, thisEQ.timeData = this_comp.timeData; end if ~isempty(this.savePath) diff --git a/applications/Binaural-HRTF/HpTF_class/itaHpTF_MS.m b/applications/Binaural-HRTF/HpTF_class/itaHpTF_MS.m index 1fedae95f223937cdacb958ed08e5c2bad7320a2..57c586d6eb292f21be381c0ef40841e9b37d25e0 100644 --- a/applications/Binaural-HRTF/HpTF_class/itaHpTF_MS.m +++ b/applications/Binaural-HRTF/HpTF_class/itaHpTF_MS.m @@ -88,6 +88,7 @@ classdef itaHpTF_MS < itaHpTF % Measure Right side MS.inputChannels = chIn(2); MS.outputChannels = chOut(2); + commandwindow resultR = MS.run; resultR.channelNames{1} = strR; @@ -108,7 +109,7 @@ classdef itaHpTF_MS < itaHpTF disp(['Measurement ',num2str(idxM),'/',num2str(this.repeat),' finished.']) disp('Put headphones off and on again. Press any key to continue.') pause - pause(2) + pause(1) else commandwindow fprintf('\nMEASUREMENT DONE!\n') @@ -118,6 +119,7 @@ classdef itaHpTF_MS < itaHpTF MS.inputChannels = chIn; HpTF = itaHpTF_Audio(this); + HpTF.selectMeas = 1:this.repeat; assignin('base', ['HpTF_' this.nameSubj], HpTF) end diff --git a/applications/HTMLhelp/ita_generate_documentation.m b/applications/HTMLhelp/ita_generate_documentation.m index 954e7322732d552b811350856489ccbd381d3957..7ce7fd863fc1b6cb90fca63a6b22ad589ca346b3 100644 --- a/applications/HTMLhelp/ita_generate_documentation.m +++ b/applications/HTMLhelp/ita_generate_documentation.m @@ -40,7 +40,8 @@ ignoreList = {'.svn', ... 'doc', ... 'GuiCallbacks', ... 'external_packages', ... - 'ExternalPackages'}; + 'ExternalPackages', ... + 'helpers'}; pathStr = genpath(sArgs.rootpath); prefixToolbox = fliplr(strtok(fliplr(sArgs.rootpath),filesep)); %get Toolbox folder name diff --git a/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorControlNanotec.m b/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorControlNanotec.m index 68d238cf689b7c2832ea39172b039103e3afaaa8..ca58b8460ccd38f50098a3c646acdf79cec04888 100644 --- a/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorControlNanotec.m +++ b/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorControlNanotec.m @@ -142,7 +142,7 @@ classdef itaMotorControlNanotec < itaMotorControl % parse the notfound options for wait [controlOptions, notFound] = ita_parse_arguments(controlOptions,notFound); - + this.wait = controlOptions.wait; for index = 1:length(this.motorList) motorposition = sArgs.(this.motorList{index}.getMotorName()); this.started(index) = this.motorList{index}.prepareMove(motorposition,notFound{:}); @@ -222,6 +222,15 @@ classdef itaMotorControlNanotec < itaMotorControl ita_verbose_info('Finished preparing',2) end + function startMove(this) + for index = 1:length(this.motorList) + if this.started(index) + this.motorList{index}.startMoveToPosition(); + end + end + this.wait4everything; + end + function startContinuousMoveNow(this) % start moves for index = 1:length(this.motorList) @@ -235,6 +244,7 @@ classdef itaMotorControlNanotec < itaMotorControl this.mIsInitialized = false; error(sprintf('Motor %s is not responding!',this.motorList{index}.getMotorName)); end + this.wait4everything this.preparedList = []; end @@ -447,7 +457,7 @@ classdef itaMotorControlNanotec < itaMotorControl else ita_verbose_info('Position NOT reached! - Check for errors!', 0); this.send_commandlist(this.failed_command_repititions); % mpo: bugfix: send_commandlist needs argument - this.isReferenced = false; +% this.isReferenced = false; end this.clear_receivedlist; this.started(1:length(this.motorIDList)) = false; diff --git a/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorNanotec_HRTFarc.m b/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorNanotec_HRTFarc.m index c1e728bc3c1b17fc920c2811a3cd46853ab51b66..13620922bf8d207b69ed899428211d56c4e97017 100644 --- a/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorNanotec_HRTFarc.m +++ b/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/itaMotorNanotec_HRTFarc.m @@ -14,18 +14,25 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec end properties(Constant, Hidden = true) + sArgs_default_motor = struct( ... - 'wait', true, ... - 'speed', 2, ... - 'VST', 'adaptiv', ... - 'limit', true, ... - 'continuous', false, ... - 'absolut', true, ... - 'closed_loop', false, ... - 'acceleration_ramp', 1, ... - 'gear_ratio', 5, ... - 'current', 100, ... - 'ramp_mode', 2 ); + 'wait', true, ... + 'speed', 2, ... + 'VST', 'adaptiv', ... + 'limit', true, ... + 'continuous', false, ... + 'absolut', true, ... + 'closed_loop', false, ... + 'acceleration_ramp',1, ... + 'gear_ratio', 5, ... + 'current', 100, ... + 'ramp_mode', 2, ... + 'P', 400, ... + 'I', 2.0, ... + 'D', 700, ... + 'P_nenner', 3, ... + 'I_nenner', 5,... + 'D_nenner', 3); end methods @@ -87,11 +94,11 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec % phasenstrom im stillstand motorControl.add_to_commandlist(sprintf('#%dr=0\r' , this.motorID)); % fehlerkorrekturmodus - motorControl.add_to_commandlist(sprintf('#%dU=1\r' , this.motorID)); + motorControl.add_to_commandlist(sprintf('#%dU=0\r' , this.motorID)); % ausschwingzeit motorControl.add_to_commandlist(sprintf('#%dO=1\r' , this.motorID)); % umkehrspiel - motorControl.add_to_commandlist(sprintf('#%dz=0\r' , this.motorID)); + motorControl.add_to_commandlist(sprintf('#%dz=5\r' , this.motorID)); % automatisches senden des status deaktivieren motorControl.add_to_commandlist(sprintf('#%dJ=0\r' , this.motorID)); @@ -125,7 +132,7 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec motorControl.add_to_commandlist(sprintf('#%dd=1\r' , this.motorID)); % Calculate and set lower speed: % stepspersecond = round(this.sArgs_default_motor.speed/0.9*this.sArgs_default_motor.gear_ratio); - motorControl.add_to_commandlist(sprintf('#%du=%d\r' , this.motorID, 3)); + motorControl.add_to_commandlist(sprintf('#%du=%d\r' , this.motorID, 1)); % Calculate and set upper speed: stepspersecond = round(this.sArgs_default_motor.speed/0.9*this.sArgs_default_motor.gear_ratio); motorControl.add_to_commandlist(sprintf('#%do=%d\r' , this.motorID, stepspersecond)); @@ -149,7 +156,7 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec sArgs = ita_parse_arguments(sArgs,varargin); end if sArgs.continuous - ret = this.prepare_move(position, 'speed', sArgs.speed,'continuous', true); + ret = this.prepare_move(position, sArgs); started = true; else % if only the phi angle is given @@ -300,19 +307,25 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec % Activate CL during movement: motorControl.add_to_commandlist(sprintf('#%d:CL_enable=2\r' , this.motorID)); % Nice values for the speed closed loop control: - pos = [0.5 1 2 3 4 8 12 16 25 32 40 50]; - vecP = [0.5 1.5 2.5 3.5 4.5 4.5 5.5 2.5 2.0 1.3 1.3 1.3]; - vecI = [0.05 0.1 0.2 0.3 0.4 0.8 1.2 1.6 2.0 2.5 2.5 2.5]; - vecD = [9 6 4 3 2 1 1 3 6 10 10 10]; - pP = polyfit(pos,vecP,5); - pI = polyfit(pos,vecI,5); - pD = polyfit(pos,vecD,5); - P = polyval(pP,this.sArgs_motor.speed); - I = polyval(pI,this.sArgs_motor.speed); - D = polyval(pD,this.sArgs_motor.speed); - P_nenner = 5; - I_nenner = 5; - D_nenner = 5; +% pos = [0.5 1 2 3 4 8 12 16 25 32 40 50]; +% vecP = [0.5 1.5 2.5 3.5 4.5 4.5 5.5 2.5 2.0 1.3 1.3 1.3]; +% vecI = [0.05 0.1 0.2 0.3 0.4 0.8 1.2 1.6 2.0 2.5 2.5 2.5]; +% vecD = [9 6 4 3 2 1 1 3 6 10 10 10]; +% pP = polyfit(pos,vecP,5); +% pI = polyfit(pos,vecI,5); +% pD = polyfit(pos,vecD,5); +% P = polyval(pP,this.sArgs_motor.speed); +% I = polyval(pI,this.sArgs_motor.speed); +% D = polyval(pD,this.sArgs_motor.speed); +% P_nenner = 5; +% I_nenner = 5; +% D_nenner = 5; + P = this.sArgs_motor.P;% (400 = default) + I = this.sArgs_motor.I;% (2 = default) + D = this.sArgs_motor.D;% (700 = default) + P_nenner = this.sArgs_motor.P_nenner; + I_nenner = this.sArgs_motor.I_nenner; + D_nenner = this.sArgs_motor.D_nenner; P_zaehler = round(P*2^P_nenner); I_zaehler = round(I*2^I_nenner); D_zaehler = round(D*2^D_nenner); @@ -322,13 +335,14 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec motorControl.add_to_commandlist(sprintf('#%d:CL_KI_v_N=%d\r' , this.motorID, I_nenner)); motorControl.add_to_commandlist(sprintf('#%d:CL_KD_v_Z=%d\r' , this.motorID, D_zaehler)); motorControl.add_to_commandlist(sprintf('#%d:CL_KD_v_N=%d\r' , this.motorID, D_nenner)); - % Nice values for the positioning closed loop control: - P = 200;% (400 = default) - I = 1.0;% (2 = default) - D = 300;% (700 = default) - P_nenner = 3; - I_nenner = 5; - D_nenner = 3; +% % Nice values for the positioning closed loop control: + + P = this.sArgs_motor.P;% (400 = default) + I = this.sArgs_motor.I;% (2 = default) + D = this.sArgs_motor.D;% (700 = default) + P_nenner = this.sArgs_motor.P_nenner; + I_nenner = this.sArgs_motor.I_nenner; + D_nenner = this.sArgs_motor.D_nenner; P_zaehler = round(P*2^P_nenner); I_zaehler = round(I*2^I_nenner); D_zaehler = round(D*2^D_nenner); @@ -338,6 +352,7 @@ classdef itaMotorNanotec_HRTFarc < itaMotorNanotec motorControl.add_to_commandlist(sprintf('#%d:CL_KI_s_N=%d\r' , this.motorID, I_nenner)); motorControl.add_to_commandlist(sprintf('#%d:CL_KD_s_Z=%d\r' , this.motorID, D_zaehler)); motorControl.add_to_commandlist(sprintf('#%d:CL_KD_s_N=%d\r' , this.motorID, D_nenner)); + % Kask V-Regler: P = 1.2, I = 0.85, D = 0.7 % Kask- P-Regler: P = 400 (default), I = 2 (default), D = 700 (default) % TODO: Send values to the motor... or we just skip the diff --git a/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/test_itaEimarMotorControl.m b/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/test_itaEimarMotorControl.m index b0e48383f145d0f0136a0c19dbaafae17cf518f6..7dda8d0213d201792334fa370d00b5551a349ae8 100644 --- a/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/test_itaEimarMotorControl.m +++ b/applications/Hardware/ITAMotorControl/ClassStuff/itaMotorControl/test_itaEimarMotorControl.m @@ -204,7 +204,7 @@ classdef test_itaEimarMotorControl < itaMeasurementTasksScan % calculate the pre angle % the pre angle is needed because the measurement setup will % not start recording imidiately - numRepetitions = this.measurementSetup.repititions; + numRepetitions = this.measurementSetup.repetitions; timePerRepetition = this.measurementSetup.twait*length(this.measurementSetup.outputChannels); speed = 360/(numRepetitions*timePerRepetition); @@ -215,7 +215,7 @@ classdef test_itaEimarMotorControl < itaMeasurementTasksScan preAngle = min(preAngle,15); preAngle = max(preAngle,8); numTotalRepetitions = numRepetitions+ceil(preAngleTime/(timePerRepetition))+9; - this.measurementSetup.repititions = numTotalRepetitions; + this.measurementSetup.repetitions = numTotalRepetitions; %prepare motors for continuous measurement this.mMotorControl.prepareForContinuousMeasurement('speed',speed,'preAngle',preAngle); diff --git a/applications/Hardware/Tracking/Optitrack/@itaOptitrack/itaOptitrack.m b/applications/Hardware/Tracking/Optitrack/@itaOptitrack/itaOptitrack.m index bea23b373d97149e33bcaa54a3f45f290f33d229..58baec899344fb0b2d193baaa5290fdaeae95fd4 100644 --- a/applications/Hardware/Tracking/Optitrack/@itaOptitrack/itaOptitrack.m +++ b/applications/Hardware/Tracking/Optitrack/@itaOptitrack/itaOptitrack.m @@ -553,7 +553,7 @@ classdef itaOptitrack < handle if Optitrack_obj.singleShot Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'ExecutionMode','singleShot'); % execute timer callback once else - Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,'ExecutionMode','fixedRate','BusyMode','drop'); + Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,'ExecutionMode','fixedSpacing','BusyMode','drop'); end else Optitrack_obj.correctRowIdx = 1; % reset row counter @@ -561,7 +561,7 @@ classdef itaOptitrack < handle Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'ExecutionMode','singleShot'); % execute timer callback once else Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,... - 'ExecutionMode','fixedRate','BusyMode','drop','TasksToExecute',32*Optitrack_obj.numFrames); % execute timer 16*Optitrack_obj.numFrames times but abort as soon as Optitrack_obj.rigidBodyLogData.data is filled + 'ExecutionMode','fixedSpacing','BusyMode','drop','TasksToExecute',32*Optitrack_obj.numFrames); % execute timer 16*Optitrack_obj.numFrames times but abort as soon as Optitrack_obj.rigidBodyLogData.data is filled end end diff --git a/applications/Hardware/Tracking/Optitrack/ita_vu2quat.m b/applications/Hardware/Tracking/Optitrack/ita_vu2quat.m index 206803bcb21f77ad715a71d4343a583d9a67eab0..9c8ab76d426264abca5bdf4089fc17fa4d220320 100644 --- a/applications/Hardware/Tracking/Optitrack/ita_vu2quat.m +++ b/applications/Hardware/Tracking/Optitrack/ita_vu2quat.m @@ -52,8 +52,15 @@ if abs(dot(v,u,2)) > 1e-5 end % normalize view/up vectors -v = normr(v); -u = normr(u); +[~,colv]=size(v); +[~,colu]=size(u); +if (colv == 1) + v(~isnan(v(:,1)),:) = v(~isnan(v(:,1)),:) ./ abs(v(~isnan(v(:,1)),:)); + u(~isnan(u(:,1)),:) = u(~isnan(u(:,1)),:) ./ abs(u(~isnan(u(:,1)),:)); +else + v = sqrt( ones ./ (sum((v(~isnan(v(:,1)),:).*v(~isnan(v(:,1)),:))')) )' * ones(1,colv).*v(~isnan(v(:,1)),:); + u = sqrt( ones ./ (sum((u(~isnan(u(:,1)),:).*u(~isnan(u(:,1)),:))')) )' * ones(1,colu).*u(~isnan(u(:,1)),:); +end % calculate side vector s = cross(v, u); @@ -65,7 +72,7 @@ qx = qw; qy = qw; qz = qw; -for idx = 1:vec_ent; +for idx = 1:vec_ent % build rotation matrix R = [s(idx,:); u(idx,:); -v(idx,:)]; diff --git a/applications/Hardware/Tracking/Optitrack/ita_vu2rpy.m b/applications/Hardware/Tracking/Optitrack/ita_vu2rpy.m index b11f7f3665cba756cadf87971b978ced5571feed..249c01971f8694c5e0cb4428d8b7e92a9b065c06 100644 --- a/applications/Hardware/Tracking/Optitrack/ita_vu2rpy.m +++ b/applications/Hardware/Tracking/Optitrack/ita_vu2rpy.m @@ -63,8 +63,15 @@ if abs(dot(v,u,2)) > 1e-5 end % normalize view/up vectors -v(~isnan(v(:,1)),:) = normr(v(~isnan(v(:,1)),:)); -u(~isnan(u(:,1)),:) = normr(u(~isnan(u(:,1)),:)); +[~,colv]=size(v); +[~,colu]=size(u); +if (colv == 1) + v(~isnan(v(:,1)),:) = v(~isnan(v(:,1)),:) ./ abs(v(~isnan(v(:,1)),:)); + u(~isnan(u(:,1)),:) = u(~isnan(u(:,1)),:) ./ abs(u(~isnan(u(:,1)),:)); +else + v = sqrt( ones ./ (sum((v(~isnan(v(:,1)),:).*v(~isnan(v(:,1)),:))')) )' * ones(1,colv).*v(~isnan(v(:,1)),:); + u = sqrt( ones ./ (sum((u(~isnan(u(:,1)),:).*u(~isnan(u(:,1)),:))')) )' * ones(1,colu).*u(~isnan(u(:,1)),:); +end % init. y = NaN(size(v,1),1); diff --git a/applications/Measurement/ClassStuff/@itaMSTF/itaMSTF.m b/applications/Measurement/ClassStuff/@itaMSTF/itaMSTF.m index 8d41395473e2de8a8de6736e4e22b150e86ecc32..9e07fffbeb2b6cfb0e0e9c5b4fcc760b45510481 100644 --- a/applications/Measurement/ClassStuff/@itaMSTF/itaMSTF.m +++ b/applications/Measurement/ClassStuff/@itaMSTF/itaMSTF.m @@ -114,7 +114,7 @@ classdef itaMSTF < itaMSPlaybackRecord fieldName = fieldnames(rmfield(this.saveobj,'dateSaved')); end - for ind = 1:numel(fieldName); + for ind = 1:numel(fieldName) try this.(fieldName{ind}) = varargin{1}.(fieldName{ind}); catch errmsg @@ -371,11 +371,10 @@ classdef itaMSTF < itaMSPlaybackRecord %% sweep rate from analytic calculation, only using sweep parameters / PDI nSamples = ita_nSamples( this.fftDegree ); - finalFreqRange = this.finalFreqRange; % 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(finalFreqRange(2)/finalFreqRange(1))/finalExcitationLength; + 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]); @@ -563,7 +562,7 @@ classdef itaMSTF < itaMSPlaybackRecord if this.minimumphasedeconvolution % get minimumphase part of deconvolution, neglect % all-pass component - [this.mCompensation, allpass_component] = ita_invert_spk_regularization(ita_extend_dat(this.raw_excitation,this.final_excitation.nSamples*factor),[min(this.freqRange(:)) max(this.freqRange(:))],'filter',this.filter); + [this.mCompensation, allpass_component] = ita_invert_spk_regularization(ita_extend_dat(this.raw_excitation,this.final_excitation.nSamples*factor),[min(this.freqRange(:)) max(this.freqRange(:))],'filter',this.filter); %#ok else this.mCompensation = ita_invert_spk_regularization(ita_extend_dat(this.raw_excitation,this.final_excitation.nSamples*factor),[min(this.freqRange(:)) max(this.freqRange(:))],'filter',this.filter); end @@ -619,7 +618,7 @@ classdef itaMSTF < itaMSPlaybackRecord token = this.(list{idx}); if isempty(token) || isa(token,'itaSuper') continue; - end; + end if ischar(token) token = ['''' token '''']; @@ -630,7 +629,7 @@ classdef itaMSTF < itaMSPlaybackRecord end else error([upper(mfilename) '.commandline: What kind of field value is this?']); - end; + end str = [str '''' list{idx} '''' ',' token ]; if idx < numel(list) @@ -643,7 +642,7 @@ classdef itaMSTF < itaMSPlaybackRecord %% Hidden methods methods(Hidden = true) - function display(this) + function display(this) %#ok this.excitation; %pre-init % Begin Display Start Line diff --git a/applications/Measurement/ImpedanceTube/Protocol/KundtGermanTemplate.tex b/applications/Measurement/ImpedanceTube/Protocol/KundtGermanTemplate.tex index 64b02dfe345dcd608062eebd2550d4a2bbeb2eea..a841b521df306339c231b8320310dc64a6b92825 100644 --- a/applications/Measurement/ImpedanceTube/Protocol/KundtGermanTemplate.tex +++ b/applications/Measurement/ImpedanceTube/Protocol/KundtGermanTemplate.tex @@ -109,7 +109,7 @@ \begin{tabular}{|p{0.43\textwidth}p{0.51\textwidth}|}\hline \multicolumn{2}{|p{0.99\textwidth}|}{ \underline{Beschreibung der Prüfeinrichtung:} \newline - Impedanzmessrohr mit 2 Zoll Bohrung aus Aluminium. Probenhalter mit verstellbarem Rohrabschluss durch massiven, abgedichteten Aluminiumkolben. Messrohr mit vier Mikrofonbohrungen (Distanzen zur ersten Position: 17~mm, 110~mm, 510~mm). Sondenmikrofon mit Abnahme des Schallfelds im Rohrmittelpunkt. Auswertbarer Frequenzbereich ca.~100 Hz bis 8 kHz durch Ausblendung der ersten zwei Rohrmoden. Breitbandige Anregung mittels Sinus-Sweeps. Auswertung nach ISO 10534-2 mit Hilfe einer im ITA entwickelten FFT-Messtechnik.}\\ + Impedanzmessrohr mit 2 Zoll Bohrung aus Aluminium. Probenhalter mit verstellbarem Rohrabschluss durch massiven, abgedichteten Aluminiumkolben. Messrohr mit vier Mikrofonbohrungen (Distanzen zur ersten Position: 17~mm, 110~mm, 510~mm). Sondenmikrofon zur Abnahme des Schallfelds im Rohrmittelpunkt. Auswertbarer Frequenzbereich ca.~100 Hz bis 8 kHz durch Ausblendung der ersten zwei Rohrmoden. Breitbandige Anregung mittels Sinus-Sweeps. Auswertung nach ISO 10534-2 mit Hilfe einer im ITA entwickelten FFT-Messtechnik.}\\ &\\ Prüfdatum: & \datumDerMessung \\ Prüfer: & \nameDesPruefers\\ diff --git a/applications/Measurement/NI_DAQ/itaMSTFni.m b/applications/Measurement/NI_DAQ/itaMSTFni.m new file mode 100644 index 0000000000000000000000000000000000000000..4c112fcb95c4b990d7f38896858ce8bf4b86aec9 --- /dev/null +++ b/applications/Measurement/NI_DAQ/itaMSTFni.m @@ -0,0 +1,351 @@ +classdef itaMSTFni < itaMSTF + % This is a class for Transfer Function or Impulse Response + % measurements with a National Instruments (NI) DAC using the DAQ toolbox. + % It supports everything that the regular itaMSTF class does, so see that for + % more info. + % + % Specific to this class: the NI setup is done with the Data + % Acquisition Toolbox, and changes to the channel setup have to be + % monitored and the NI session has to be updated accordingly. + % + % See also: itaMSTF + + % Author: Markus Mueller-Trapet 2017 - markus.mueller-trapet@nrc.ca + + properties(Access = public, Hidden = true) % internal variables + niSession = []; % to store information about NI card setup + end + + methods + + %% CONSTRUCT / INIT / EDIT / COMMANDLINE + + function this = itaMSTFni(varargin) + % itaMSTFni - Constructs an itaMSTFni object. + if nargin == 0 + + % For the creation of itaMSTFni objects from commandline strings + % like the ones created with the commandline method of this + % class, 2 or more input arguments have to be allowed. All + % desired properties have to be given in pairs of two, the + % first element being an identifying string which will be used + % as field name for the property, and the value of the + % specified property. + elseif nargin >= 2 + if ~isnatural(nargin/2) + error('Even number of input arguments expected!'); + end + + % For all given pairs of two, use the first element as + % field name, the second one as value. The validity of the + % field names will NOT be checked. + for idx = 1:2:nargin + this.(varargin{idx}) = varargin{idx+1}; + end + + % Only one input argument is required for the creation of an + % itaMSTFni class object from a struct, created by the saveobj + % method, or as a copy of an already existing itaMSTFni class + % object. In the latter case, only the properties contained in + % the list of saved properties will be copied. + elseif isstruct(varargin{1}) || isa(varargin{1},'itaMSTF') + % Check type of given argument and obtain the list of saved + % properties accordingly. + if isa(varargin{1},'itaMSTF') + %The save struct is obtained by using the saveobj + % method, as in the case in which a struct is given + % from the start (see if-case above). + if isa(varargin{1},'itaMSTFni') + deleteDateSaved = true; + else + deleteDateSaved = false; + end + varargin{1} = saveobj(varargin{1}); + % have to delete the dateSaved field to make clear it + % might be from an inherited class + if deleteDateSaved + varargin{1} = rmfield(varargin{1},'dateSaved'); + end + end + if isfield(varargin{1},'dateSaved') + varargin{1} = rmfield(varargin{1},'dateSaved'); + fieldName = fieldnames(varargin{1}); + else %we have a class instance here, maybe a child + fieldName = fieldnames(rmfield(this.saveobj,'dateSaved')); + end + + for ind = 1:numel(fieldName) + try + this.(fieldName{ind}) = varargin{1}.(fieldName{ind}); + catch errmsg + disp(errmsg); + end + end + else + error('itaMSTFni::wrong input arguments given to the constructor'); + end + + % Define listeners to automatically call the init function of + % this class in case of a change in the below specified + % properties. + addlistener(this,'inputChannels','PostSet',@this.init); + addlistener(this,'outputChannels','PostSet',@this.init); + end + + function this = edit(this) + % edit - Start GUI. + % + % This function calls the itaMSTFni GUI. + + this = ita_mstfni_gui(this); + end + + function this = init(this,varargin) + % init - Initialize the itaMSTFni class object. + init@itaMSTF(this); + % if this is an object loaded from disk (niSession is empty), + % we do not need to initialize the card + if ~isempty(this.niSession) + this.niSession = init_NI_card(this); + end + end + + function checkready(this) + %check if the instance is ready for measurement run and ask for + %missing entries + if isempty(this.inputChannels) || isempty(this.outputChannels) + this.edit; + end + % has the NI session been initialized + if isempty(this.niSession) || isempty(this.niSession.Channels) + this.niSession = init_NI_card(this); + end + % has the samplingRate been changed (NI rate is no exact) + if abs(this.niSession.Rate - this.samplingRate) > 1 + this.niSession = init_NI_card(this); + end + end + + function MS = calibrationMS(this) + this.checkready; + % call the parent function + MS = calibrationMS@itaMSTF(this); + % convert to instance of this class + MS = itaMSTFni(MS); + % release NI hardware to enable measurement with calibrationMS + this.niSession.release; + end + + function [result, max_rec_lvl] = run_raw(this) + % run_raw - Run measurement + this.checkready; + singleprecision = strcmpi(this.precision,'single'); % Bool for single precision for portaudio. + + result = ita_NI_daq_run(this.final_excitation,this.niSession,'InputChannels',this.inputChannels, ... + 'OutputChannels', this.outputChannels,'repeats',1,... + 'latencysamples',this.latencysamples,'singleprecision',singleprecision); + + if this.outputVoltage ~= 1 % only if output is calibrated + result.comment = [result.comment ' @' num2str(round(this.outputVoltage*1000)/1000) 'Vrms']; + end + max_rec_lvl = max(abs(result.timeData),[],1); + end + + function [result, max_rec_lvl] = run_latency(this) + % call parent function + [result, max_rec_lvl] = run_latency@itaMSTF(this); + end + + function this = calibrate_input(this,elementIds) + % have to do this here because of different run function + % do only specific elements (e.g. only AD) + this.checkready + if ~exist('elementIds','var') + elementIds = 1:3; + else + elementIds = unique(min(3,max(1,elementIds))); + end + % and only active channels + inputChannels = this.inputChannels; + imcIdx = zeros(numel(inputChannels),1); + for chIdx = 1:numel(inputChannels) + imcIdx(chIdx) = find(this.inputMeasurementChain.hw_ch == inputChannels(chIdx)); + end + tmpChain = this.inputMeasurementChain(imcIdx); + % element by element + for iElement = elementIds + for iCh = 1:numel(imcIdx) + this.inputChannels = inputChannels(iCh); + if numel(tmpChain(iCh).elements) >= iElement + hw_ch = tmpChain(iCh).hardware_channel; + disp(['Calibration of sound card channel ' num2str(hw_ch)]) + % go thru all elements of the chain and calibrate + if tmpChain(iCh).elements(iElement).calibrated ~= -1 % only calibratable devices + disp([' Calibration of ' upper(tmpChain(iCh).elements(iElement).type) ' ' tmpChain(iCh).elements(iElement).name]) + this.inputChannels = inputChannels(iCh); + [tmpChain(iCh).elements(iElement).sensitivity] = measurement_chain_elements_calibration_ni(this.niSession,tmpChain(iCh),iElement); %calibrate each element + end + end + end + end + this.inputMeasurementChain(imcIdx) = tmpChain; + this.inputChannels = inputChannels; + disp('****************************** FINISHED *********************************') + end + + function this = calibrate_output(this,input_chain_number) + % have to do this here because of different run function + % Calibrates all output chains, using only the first + % (hopefully calibrated) input chain. Input chain calibration + this.checkready + if ~exist('input_chain_number','var') + input_chain_number = find(this.inputMeasurementChain.hw_ch == this.inputChannels(1)); + end + ita_verbose_info(['Calibrating using input channel ' num2str(this.inputMeasurementChain(input_chain_number).hardware_channel)],1); + + MS = this.calibrationMS; % Get new simple Measurement Setup for calibration. See above. + MS.inputChannels = MS.inputChannels(input_chain_number); + mco = this.outputMeasurementChain; % Get all output measurement chains. + outChannels = this.outputChannels; % Get all output channels. + + % The calibration of the multiple output measurement chains / + % outout channels will be executed one-by-one. + for outIdx = 1:numel(outChannels) + chIdx = find(mco.hw_ch == outChannels(outIdx)); % Return single index of entry in 'mco', equal to the out channel, which is to be calibrated. + MS.outputMeasurementChain = mco(chIdx); % Set Measurement Setup's single output chain to match the one, which is to be calibrated. + MS.outputChannels = outChannels(outIdx); % Set Measurement Setup's single output channel to match the one, which is to be calibrated. + + % Execute calibration for every single element in the + % current output measurement chain. + % 'ita_mstfoutput_calibration' determines if the object can + % be calibrated at all. + for ele_idx = 1:length(mco(chIdx).elements) + MS = measurement_chain_output_calibration_ni(MS,input_chain_number,ele_idx); + end + + % if there was no latency info before, copy it from the + % calibrationMS because latency was measured in the output + % calibration routine + if this.latencysamples == 0 + this.latencysamples = MS.latencysamples; + end + % Put the current calibrated measurement chain back into + % its appropriate position in the list of all output + % measurment chains. + mco(chIdx) = MS.outputMeasurementChain; + end + this.outputMeasurementChain = mco; % Copy over the list of all calibrated output measurement chains into the real Measurement Setup. + % release hardware for the standard object + if ~isempty(MS.niSession) % only if something was measured + MS.niSession.release; + end + end + + function [niSession,inputChannels,outputChannels] = init_NI_card(this) + % uses Christoph Hoellers's (hoellerc@nrc.ca) code for initilaization of NI session + % for now only as simple DAC, so only Voltage type + + % Initialization + [inputChannels,outputChannels,niDevices,rateLimits] = ita_get_ni_deviceinfo(); + if isempty(niDevices) % no card attached + error('Cannot initialize NI card, maybe not connected?') + end + % create session (will be stored) + niSession = daq.createSession('ni'); + % turn off useless warning + warning('off','daq:Session:clockedOnlyChannelsAdded') + + if this.samplingRate < rateLimits(1) + warning(['Device does not support a sampling rate of ' num2str(this.samplingRate) ', changing to lower limit of ' num2str(rateLimits(1))]); + this.samplingRate = rateLimits(1); + elseif this.samplingRate > rateLimits(2) + warning(['Device does not support a sampling rate of ' num2str(this.samplingRate) ', changing to upper limit of ' num2str(rateLimits(2))]); + this.samplingRate = rateLimits(2); + end + niSession.Rate = this.samplingRate; + + % INPUT + % set channel data from MS + if any(this.inputChannels > numel(inputChannels.name)) + error(['Your device does not have ' num2str(max(this.inputChannels)) ' input channels!']); + else + inputChannels.isActive = ismember(1:numel(inputChannels.name),this.inputChannels); + end + + % Add analog input channels + for iChannel = 1:numel(inputChannels.name) + if inputChannels.isActive(iChannel) + iDevice = inputChannels.mapping{iChannel}(1); + iDeviceChannel = inputChannels.mapping{iChannel}(2); + niSession.addAnalogInputChannel(get(niDevices(iDevice),'ID'),iDeviceChannel-1,inputChannels.type{iChannel}); + niSession.Channels(end).Name = inputChannels.name{iChannel}; + % set to AC coupling to get rid of large DC offset + niSession.Channels(end).Coupling = 'AC'; + % if any(strcmpi(Channels.type{iDevice,iChannel},{'Accelerometer' 'Microphone'})) + % niSession.Channels(end).Sensitivity = Channels.sensitivity(iDevice,iChannel); + % end + end + end + + % OUTPUT + % set channel data from MS + if any(this.outputChannels > numel(outputChannels.name)) + error(['Your device does not have ' num2str(max(this.outputChannels)) ' output channels!']); + else + outputChannels.isActive = ismember(1:numel(outputChannels.name),this.outputChannels); + end + + % Add analog output channels + for iChannel = 1:numel(outputChannels.name) + if outputChannels.isActive(iChannel) + iDevice = outputChannels.mapping{iChannel}(1); + iDeviceChannel = inputChannels.mapping{iChannel}(2); + niSession.addAnalogOutputChannel(get(niDevices(iDevice),'ID'),iDeviceChannel-1,outputChannels.type{iChannel}); + niSession.Channels(end).Name = inputChannels.name{iChannel}; + end + end + + end % function + + function sObj = saveobj(this) + % saveobj - Saves the important properties of the current + % measurement setup to a struct. + + sObj = saveobj@itaMSTF(this); + % Get list of properties to be saved for this measurement + % class. + propertylist = itaMSTFni.propertiesSaved; + + % Write the content of every item in the list of the to be saved + % properties into its own field in the save struct. + for idx = 1:numel(propertylist) + sObj.(propertylist{idx}) = this.(propertylist{idx}); + end + end + end + + methods(Static, Hidden = true) + + function result = propertiesSaved + % propertiesSaved - Creates a list of all the properties to be + % saved of the current measurement setup. + % + % This function gets the list of all + % properties to be saved during the saving process. + + % Get list of saved properties for this class. + result = {}; + end + + function this = loadobj(sObj) + % loadobj - Creates a new measurement setup and loads the + % properties of a save struct into it. + % + % This function creates a new measurement setup by calling the + % class constructor and passes it the specified save struct. + + this = itaMSTFni(sObj); + end + end + +end % classdef diff --git a/applications/Measurement/NI_DAQ/ita_NI_daq_run.m b/applications/Measurement/NI_DAQ/ita_NI_daq_run.m new file mode 100644 index 0000000000000000000000000000000000000000..b529beaaa19cc5fbcd8a7b63f9411b327c23acc7 --- /dev/null +++ b/applications/Measurement/NI_DAQ/ita_NI_daq_run.m @@ -0,0 +1,203 @@ +function varargout = ita_NI_daq_run(varargin) +% ITA_NI_DAQ_RUN play and record sound with NI card and DAQ toolbox (adapted from ita_portaudio_run) +% see ita_portaudio_run for options and help + +% Autor: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca +% Created: 13-Apr-2017 + +%% ITA Toolbox preferences for verbose level +verboseMode = ita_preferences('verboseMode'); + +%% Input checks +if isa(varargin{1},'itaAudio') + % if a signal is the first input, then playback that signal + sArgs.pos1_data = 'itaAudioTime'; + playback = true; +else + % otherwise just record for the given number of samples + sArgs.pos1_data = 'numeric'; + playback = false; +end + +% second argument has to be the DAQ NI session +sArgs.pos2_niSession = 'anything'; + +% only do recording if requested +if nargout > 0 + record = true; +else + record = false; +end + +%% Init Settings +sArgs.normalizeinput = false; +sArgs.normalizeoutput = false; + +sArgs.inputchannels = []; +sArgs.outputchannels = []; +sArgs.recsamples = -1; +sArgs.repeats = 1; +sArgs.samplingRate = -1; % will always be taken from NI session +% cancel button and monitor not supported currently, but kept for compatibility +sArgs.cancelbutton = ita_preferences('portAudioMonitor'); % -1:automatic; 0:off; 1:on +sArgs.latencysamples = 0; +sArgs.singleprecision = false; + +[data, niSession, sArgs] = ita_parse_arguments(sArgs,varargin); + +if ~playback && sArgs.recsamples == -1 + sArgs.recsamples = data; +end + +sArgs.samplingRate = round(niSession.Rate); % NI rate is not exact + +in_channel_vec = sArgs.inputchannels; +out_channel_vec = sArgs.outputchannels; +normalize_input = sArgs.normalizeinput; +normalize_output = sArgs.normalizeoutput; + +if ~playback && ~record + error('ITA_NI_DAQ_RUN:What do you want? Play or Record, or both? You need to specify an input and/or an output!') +end + +if playback + % Extract channels to play + % are there enough channels to play + if data.nChannels == 1 && (length(out_channel_vec) > 1) %playback the same on alle channels activated + ita_verbose_info('Oh Lord. I will playback the same on all channels.', 2) + nChannelsToPlay = length(out_channel_vec); + data = data.ch(ones(1,nChannelsToPlay)); %duplicated data to all channels + elseif data.nChannels < length(out_channel_vec) %too many output channels for this input data + ita_verbose_info('Not enough channels in data to play.',0) + out_channel_vec = out_channel_vec(1:data.nChannels); + elseif data.nChannels > length(out_channel_vec) + ita_verbose_info('Too many channels in data file, I will only play the first ones',1); + data = ita_split(data,1:length(out_channel_vec)); + end + + % Show channel data if requested + if verboseMode==2 + ita_metainfo_show_channelnames(data); + end + + % Check levels - Normalizing + peak_value = max(max(abs(data.timeData))); + if (peak_value > 1) || (normalize_output) + ita_verbose_info('Oh Lord! Levels too high for playback. Normalizing...',0) + data = ita_normalize_dat(data); + end +end + +%% Extend excitation to compensate soundcard latency +if playback && record && ~isempty(sArgs.latencysamples) && (sArgs.latencysamples > 0) + latencysamples = sArgs.latencysamples; + data = ita_extend_dat(data,data.nSamples+latencysamples,'forcesamples'); +end + +% record as many samples as are in the playback signal +if playback + sArgs.recsamples = data.nSamples; +end + +% Full (double) precision +recordDatadat = zeros(sArgs.recsamples,numel(in_channel_vec),sArgs.repeats); +if sArgs.singleprecision % only single precision + recordDatadat = single(recordDatadat); +end + +% run measurement, possibly repeated +for idrep = 1:sArgs.repeats + if playback && record + niSession.queueOutputData(double(data.time)); + ita_verbose_info('start playback and record',1) + elseif record + niSession.queueOutputData(zeros(sArgs.recsamples,1)); + ita_verbose_info('start record',1) + elseif playback + niSession.queueOutputData(double(data.time)); + ita_verbose_info('start playback',1) + else + error('ITA_NI_DAQ_run:No input and no output, what should I do?') + end + pause(0.01) + % do the measurement + if ~sArgs.singleprecision % Full (double) precision + recordDatadat(:,:,idrep) = niSession.startForeground(); + else + recordDatadat(:,:,idrep) = single(niSession.startForeground()); + end +end % loop for repeats + +ita_verbose_info('playback/record finished ',1); + +recordData = itaAudio(); +recordData.dataType = class(recordDatadat); +recordData.dataTypeOutput = class(recordDatadat); +% Check if we need to average multiple measurements: +if size(recordDatadat,3) > 1 + % average: + recordData.timeData = mean(recordDatadat,3); +else + % no average: (This saves memory!) + recordData.timeData = recordDatadat; +end +recordData.samplingRate = sArgs.samplingRate; + +for idx = 1:numel(in_channel_vec) + recordData.channelNames{idx} = ['Ch ' int2str(in_channel_vec(idx))]; +end + +%% Remove Latency +if playback && record && ~isempty(sArgs.latencysamples) && (sArgs.latencysamples > 0) + recordData = ita_extract_dat(recordData,recordData.nSamples-latencysamples,'firstsample',latencysamples+1); +end + +%% Check for clipping and other errors +clipping = false; +if record + if any(any(abs(recordData.timeData)>=10)) % Check for clipping (NI card actually handles up to 10Vpk) + ita_verbose_info('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!',0); + ita_verbose_info('!ITA_NI_DAQ_RUN:Careful, Clipping!',0); + ita_verbose_info('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!',0); + clipping = true; + end + % This check for singularities needs a lot of memory! TODO: find a + % better solution!: + if any(isnan(recordData.timeData(:))) || any(isinf(recordData.timeData(:))) %Check for singularities + ita_verbose_info('There are singularities in the audio signal! You''d better check your settings!',0) + clipping = true; + end + if any(all(recordData.timeData == 0,1)) && record + ita_verbose_info('There are empty channels in the audio signal! You''d better check your settings!',0) + end + % maximum for each channel + maxData = max(abs(recordData.timeData),[],1); + [channelMax, indexMax] = max(maxData); + + % jri: to detect non working microphones etc, the minimum of the + % maximums on the channels is also outputted + if length(in_channel_vec) > 1 + [channelMin, indexMin] = min(maxData); + ita_verbose_info(['Minimum digital level: ' int2str(20*log10(channelMin)) ' dBFS on channel: ' int2str(in_channel_vec(indexMin))],0); + end + ita_verbose_info(['Maximum digital level: ' int2str(20*log10(channelMax)) ' dBFS on channel: ' int2str(in_channel_vec(indexMax))],0); +end + +%% Add history line +infosforhistory = struct('PlayDevice','NI','Play_Channels',out_channel_vec,'RecDevice','NI','Rec_Channels',in_channel_vec,'Sampling_Rate',niSession.Rate,'Normalize_Input',normalize_input,'Normalize_Output',0,'Rec_Samples',sArgs.recsamples,'Block',1,'Repeats',sArgs.repeats); +recordData = ita_metainfo_add_historyline(recordData,'ita_NI_daq_run',[{data}; ita_struct2arguments(infosforhistory)],'withsubs'); + +if clipping + recordData = ita_metainfo_add_historyline(recordData,'!!!ITA_NI_DAQ_RUN:Careful, clipping or something else went wrong!!!'); + recordData = ita_errorlog_add(recordData,'!!!ITA_NI_DAQ_RUN:Careful, clipping or something else went wrong!!!'); +end + +%% Find output parameters +if nargout ~= 0 + if normalize_input + ita_normalize_dat(recordData) + end + varargout{1} = recordData; +end + +end % function \ No newline at end of file diff --git a/applications/Measurement/NI_DAQ/ita_channelselect_gui_ni.m b/applications/Measurement/NI_DAQ/ita_channelselect_gui_ni.m new file mode 100644 index 0000000000000000000000000000000000000000..2662b4ef003d850f6b1e662a3f8f0819e0a42263 --- /dev/null +++ b/applications/Measurement/NI_DAQ/ita_channelselect_gui_ni.m @@ -0,0 +1,212 @@ +function varargout = ita_channelselect_gui_ni( varargin ) +% Show gui to select channels for input and output (NI hardware) +% +% Syntax: [in_ch out_ch] = ita_channelselect_gui_ni() +% [in_ch out_ch] = ita_channelselect_gui_ni(in_ch_preselct, out_ch_preselect) +% + +% +% 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. +% + +% Author: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca +% Created: 10-May-2017 + +%% Get maximum input and output channels +% Initialization +[inputChannels,outputChannels] = ita_get_ni_deviceinfo(); + +out = numel(outputChannels.name); +out = max(out,0); %Filter '-1' in case that no device is selected +in = numel(inputChannels.name); +in = max(in,0); %Filter '-1' in case that no device is selected +if nargin == 0 + varargin{1} = 1:in; varargin{2} = 1:out; +end +sArgs = struct('pos1_a','int','pos2_b','int','onlyinput','false','onlyoutput','false'); +[in_ch,out_ch,sArgs] = ita_parse_arguments(sArgs,varargin); + +if sArgs.onlyinput + out = 0; +elseif sArgs.onlyoutput + in = 0; +end + +%% Draw - 8 channels in a row +if ~mod(out,8) + out_height = out/8; +else + out_height = (out - mod(out,8))/8 + 1; +end +if ~mod(in,8) + in_height = in/8; +else + in_height = (in - mod(in,8))/8 + 1; +end + +%figure +height = (2 + in_height + 1 + out_height + 1) * 40; +ch_width = 70; +width = (8+3)*(ch_width - 15); +persistent hFig +if ~isempty(hFig) && ishandle(hFig) && strcmpi(get(hFig,'Name'),'PortAudio IO Channels') + close(hFig) +end +hFig = figure('Name','PortAudio IO Channels',... + 'Position',[300 200 width height],... + 'MenuBar','none', ... + 'Toolbar','none', ... + 'HandleVisibility','on',... + 'NumberTitle','off', ... + 'Color', [0.8 0.8 0.8]); + +%ITA toolbox logo with grey background +a_im = importdata(which('ita_toolbox_logo.png')); +image(a_im);axis off +set(gca,'Units','pixel', 'Position', [20 20 200 35]); %TODO: later set correctly the position + +%% Headline Texts +if in + uicontrol(... + 'Parent', hFig,... + 'Position',[20 (1+out_height+1+in_height+1)*40-20 200 20],... + 'HorizontalAlignment','left',... + 'String','Input Channels',... + 'FontSize',13,... + 'FontWeight','bold',... + 'Style', 'text',... + 'ForegroundColor', [0 0 0],... + 'BackgroundColor', [0.8 0.8 0.8]); +end +if out + uicontrol(... + 'Parent', hFig,... + 'Position',[20 (1+out_height+1)*40-20 200 20],... + 'HorizontalAlignment','left',... + 'String','Output Channels',... + 'FontSize',13,... + 'FontWeight','bold',... + 'Style', 'text',... + 'ForegroundColor', [0 0 0],... + 'BackgroundColor', [0.8 0.8 0.8]); +end + +%% Input Channels +for jdx = 1:in + uicontrol(... + 'Parent', hFig,... + 'Position',[20+(mod(jdx-1,8))*ch_width (1+out_height+1+in_height-(jdx-mod(jdx-1,8)+1)/8)*40+2 30 15],... + 'HorizontalAlignment','right',... + 'String', ['Ch' num2str(jdx)],... + 'Style', 'text',... + 'ForegroundColor', [1 .1 .1],... + 'BackgroundColor', [0.8 0.8 0.8]); +end +if isempty(in_ch) + in_ch = 1:in; +end +for jdx = 1:in + h(jdx).in = uicontrol(... + 'Parent', hFig,... + 'Position',[50 + (mod(jdx-1,8))*ch_width (1+out_height+1+in_height-(jdx-mod(jdx-1,8)+1)/8)*40-5 20 30],... + 'Style', 'checkbox',... + 'Value',ismember(jdx,in_ch),... + 'ForegroundColor', [0 0 .7],... + 'BackgroundColor', [0.8 0.8 0.8]); +end + + +%% Output Channels +for jdx = 1:out + uicontrol(... + 'Parent', hFig,... + 'Position',[20+(mod(jdx-1,8))*ch_width (1+out_height-(jdx-mod(jdx-1,8)+1)/8)*40+2 30 15],... + 'HorizontalAlignment','right',... + 'String', ['Ch' num2str(jdx)],... + 'Style', 'text',... + 'ForegroundColor', [.1 0.6 .1],... + 'BackgroundColor', [0.8 0.8 0.8]); +end +if isempty(out_ch) + out_ch = 1:out; +end +for jdx = 1:out + + h(jdx).out = uicontrol(... + 'Parent', hFig,... + 'Position',[50+(mod(jdx-1,8))*ch_width (1+out_height-(jdx-mod(jdx-1,8)+1)/8)*40-5 20 30],... + 'Style', 'checkbox',... + 'Value',ismember(jdx,out_ch),... + 'ForegroundColor', [0 0 .7],... + 'BackgroundColor', [0.8 0.8 0.8]); +end + +%% Buttons +% Cancel Button +uicontrol(... + 'Parent', hFig, ... + 'Position',[330 20 80 30],... + 'String', 'Cancel',... + 'Style', 'pushbutton',... + 'Callback', @CancelButtonCallback,... + 'BackgroundColor', [0.7 0.7 0.7]); + +% Ok Button +uicontrol(... + 'Parent', hFig, ... + 'Position',[430 20 80 30],... + 'String', 'OK',... + 'Style', 'pushbutton',... + 'Callback', @OkayButtonCallback,... + 'BackgroundColor', [0.7 0.7 0.7]); + +uiwait(hFig); + +%% Callbacks + function CancelButtonCallback(hObject,eventdata) %#ok + uiresume(gcf); + close(hFig) + %choose output arguments + if sArgs.onlyinput + varargout{1} = in_ch; + elseif sArgs.onlyoutput + varargout{1} = out_ch; + else + varargout{1} = in_ch; + varargout{2} = out_ch; + end + + return; + end + + function OkayButtonCallback(hObject,eventdata) %#ok + in_ch = []; + for jdx2 = 1:in + if get(h(jdx2).in,'Value') == get(h(jdx2).in,'Max') + in_ch = [in_ch jdx2]; %#ok<*AGROW> + end + end + out_ch = []; + for jdx2 = 1:out + if get(h(jdx2).out,'Value') == get(h(jdx2).out,'Max') + out_ch = [out_ch jdx2]; + end + end + uiresume(gcf); + close(hFig) + + %choose output arguments + if sArgs.onlyinput + varargout{1} = in_ch; + elseif sArgs.onlyoutput + varargout{1} = out_ch; + else + varargout{1} = in_ch; + varargout{2} = out_ch; + end + + return; + end + +end diff --git a/applications/Measurement/NI_DAQ/ita_get_ni_deviceinfo.m b/applications/Measurement/NI_DAQ/ita_get_ni_deviceinfo.m new file mode 100644 index 0000000000000000000000000000000000000000..de3363ad78a445dfd188f1838f05b8ef9bed063d --- /dev/null +++ b/applications/Measurement/NI_DAQ/ita_get_ni_deviceinfo.m @@ -0,0 +1,58 @@ +function [inputChannels,outputChannels,niDevices,rateLimits] = ita_get_ni_deviceinfo(inputSens) +% determine NI hardware setup for maximum input and output capabilities + +% Author: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca +% Created: 10-May-2017 + +if ~nargin + inputSens = 0.01; +end + +%% gather device info +niDevices = daq.getDevices(); % returns an object with all the NI cards +nDevices = length(niDevices); + +% samplingRate limits +rateLimits = [1 192000]; % set max defaults + +%% input channels +inputChannelIdx = 0; +for iDevice = 1:nDevices + for iSub = 1:numel(niDevices(iDevice).Subsystems) + % only have to do this here, as we go through all subsystems + tmpLimit = round(niDevices(iDevice).Subsystems(iSub).RateLimit); + rateLimits(1) = max(rateLimits(1),min(tmpLimit)); + rateLimits(2) = min(rateLimits(2),max(tmpLimit)); + if ~isempty(strfind(lower(niDevices(iDevice).Subsystems(iSub).SubsystemType),'input')) + for iChannel = 1:niDevices(iDevice).Subsystems(iSub).NumberOfChannelsAvailable + inputChannelIdx = inputChannelIdx + 1; + inputChannels.mapping(inputChannelIdx) = {[iDevice,iChannel]}; + inputChannels.name{inputChannelIdx} = [niDevices(iDevice).ID '_' niDevices(iDevice).Subsystems(iSub).ChannelNames{iChannel,1}]; + inputChannels.type{inputChannelIdx} = niDevices(iDevice).Subsystems(iSub).DefaultMeasurementType; + inputChannels.sensitivity(inputChannelIdx) = inputSens; + end + end + end +end +% default active state is zero +inputChannels.isActive = zeros(inputChannelIdx,1); + +%% output channels +outputChannelIdx = 0; +for iDevice = 1:nDevices + for iSub = 1:numel(niDevices(iDevice).Subsystems) + if ~isempty(strfind(lower(niDevices(iDevice).Subsystems(iSub).SubsystemType),'output')) + for iChannel = 1:niDevices(iDevice).Subsystems(iSub).NumberOfChannelsAvailable + outputChannelIdx = outputChannelIdx + 1; + outputChannels.mapping(outputChannelIdx) = {[iDevice,iChannel]}; + outputChannels.name{outputChannelIdx} = [niDevices(iDevice).ID '_' niDevices(iDevice).Subsystems(iSub).ChannelNames{iChannel,1}]; + outputChannels.type{outputChannelIdx} = niDevices(iDevice).Subsystems(iSub).DefaultMeasurementType; + end + end + end +end +% default active state is zero +outputChannels.isActive = zeros(outputChannelIdx,1); + + +end % function \ No newline at end of file diff --git a/applications/Measurement/NI_DAQ/ita_mstfni_gui.m b/applications/Measurement/NI_DAQ/ita_mstfni_gui.m new file mode 100644 index 0000000000000000000000000000000000000000..f871908789e48a3407b5026b9c772ebbb27995be --- /dev/null +++ b/applications/Measurement/NI_DAQ/ita_mstfni_gui.m @@ -0,0 +1,186 @@ +function varargout = ita_mstfni_gui(varargin) +%ITA_MSTFNI_GUI - Edit a measurement setup for transfer functions (NI hardware) + +% +% This file is part of the application Measurement for the ITA-Toolbox. All rights reserved. +% You can find the license for this m-file in the application folder. +% + + +% Author: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca +% Created: 10-May-2017 + + +pList = []; +argList = []; +%GUI Init + +if nargin == 1 + MS = varargin{1}; +else + MS = itaMSTFni; +end + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +ele = numel(pList)+1; +pList{ele}.datatype = 'text'; +pList{ele}.description = 'Basic settings'; +pList{ele}.color = 'black'; + +ele = numel(pList)+1; +pList{ele}.description = 'Preferences'; +pList{ele}.helptext = 'Call ita_preferences()'; +pList{ele}.datatype = 'simple_button'; +pList{ele}.default = ''; +pList{ele}.callback = 'ita_preferences();'; + +ele = numel(pList)+1; +pList{ele}.description = 'Input Channels'; +pList{ele}.helptext = 'Vector with the input channel numbers. The order specified here is respected!'; +pList{ele}.datatype = 'int_result_button'; +pList{ele}.default = MS.inputChannels; +if isempty(pList{ele}.default) + pList{ele}.default = 1; +end +pList{ele}.callback = 'ita_channelselect_gui_ni([$$],[],''onlyinput'')'; +argList = [argList {'inputChannels'}]; + + +ele = numel(pList)+1; +pList{ele}.description = 'Output Channels'; +pList{ele}.helptext = 'Vector with the output channel numbers. The order specified here is respected!'; +pList{ele}.datatype = 'int_result_button'; +pList{ele}.default = MS.outputChannels; +if isempty(pList{ele}.default) + pList{ele}.default = 1; +end +pList{ele}.callback = 'ita_channelselect_gui_ni([],[$$],''onlyoutput'')'; +argList = [argList {'outputChannels'}]; + + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +ele = numel(pList)+1; +pList{ele}.datatype = 'text'; +pList{ele}.description = 'Signal Specifications'; + +ele = numel(pList)+1; +pList{ele}.description = 'FFT Degree'; +pList{ele}.helptext = 'Length of the signal (2^fft_degree samples)'; +pList{ele}.datatype = 'int'; +pList{ele}.default = MS.fftDegree; +argList = [argList {'fftDegree'}]; + + +ele = numel(pList)+1; +pList{ele}.description = 'Frequency Limits [Hz]'; +pList{ele}.helptext = 'The sweep will start at low frequency and rise up to high frequency'; +pList{ele}.datatype = 'int'; +pList{ele}.default = MS.freqRange; +argList = [argList {'freqRange'}]; + +ele = numel(pList)+1; +pList{ele}.description = 'Signal Type'; +pList{ele}.helptext = 'Exponential/logarithmic sweeps and linear sweeps can be choosen.)'; +pList{ele}.datatype = 'char_popup'; +pList{ele}.default = MS.type; +pList{ele}.list = 'exp|lin|noise'; +argList = [argList {'type'}]; + +ele = numel(pList)+1; +pList{ele}.description = 'Stop Margin [s]'; +pList{ele}.helptext = 'This is the time of silence in the end of the sweep. It should be longer than the reverberation time at the highest frequency.'; +pList{ele}.datatype = 'int'; +pList{ele}.default = MS.stopMargin; +argList = [argList {'stopMargin'}]; + + +ele = numel(pList)+1; +pList{ele}.description = 'Output Amplification [dBFS]'; +pList{ele}.helptext = 'Attenuate your sweep by this value in dB Full Scale. Compensation will follow, we will take care of this.'; +pList{ele}.datatype = 'char'; +pList{ele}.default = MS.outputamplification; +argList = [argList {'outputamplification'}]; + + +ele = numel(pList)+1; +pList{ele}.description = 'Comment'; +pList{ele}.helptext = 'Give your child a name'; +pList{ele}.datatype = 'char_long'; +pList{ele}.default = MS.comment; +argList = [argList {'comment'}]; + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +ele = numel(pList)+1; +pList{ele}.datatype = 'text'; +pList{ele}.description = 'Advanced settings'; +pList{ele}.color = 'red'; + +ele = numel(pList)+1; +pList{ele}.description = 'Pause before measurements'; +pList{ele}.helptext = 'Time in seconds the routine waits before each measurement'; +pList{ele}.datatype = 'int'; +pList{ele}.default = MS.pause; +argList = [argList {'pause'}]; + + +ele = numel(pList)+1; +pList{ele}.description = 'Number of Averages'; +pList{ele}.helptext = 'How many measurements should be averaged for the final results?'; +pList{ele}.datatype = 'int'; +pList{ele}.default = MS.averages; +argList = [argList {'averages'}]; + + +ele = numel(pList)+1; +pList{ele}.description = 'Linear Deconvolution'; +pList{ele}.helptext = 'Standard is cyclic deconvolution. To seperate distortion from the noise tail use linear deconvolution. Signal length is therefore doubled.'; +pList{ele}.datatype = 'bool'; +pList{ele}.default = MS.lineardeconvolution; +argList = [argList {'lineardeconvolution'}]; + +ele = numel(pList)+1; +pList{ele}.description = 'Output Equalization'; +pList{ele}.helptext = 'Do a broadband output equalization (only with calibrated output measurement chain)'; +pList{ele}.datatype = 'bool'; +pList{ele}.default = MS.outputEqualization; +argList = [argList {'outputEqualization'}]; + +ele = numel(pList)+1; +pList{ele}.description = 'Measurement Chain'; +pList{ele}.helptext = 'Whether to use the measurement chain functionality (calibration)'; +pList{ele}.datatype = 'bool'; +pList{ele}.default = MS.useMeasurementChain; +argList = [argList {'useMeasurementChain'}]; + + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +%call gui +pList = ita_parametric_GUI(pList,[mfilename ' - Modify an itaMSTF']); +pause(0.02); %wait for GUI to close first + +% Check output amplifiction +if isempty(pList) %user cancelled + varargout{1} = []; + return; +end + +%% settings to MSTF +%reorder first, useMeasurementChain does not work otherwise +pList = [pList(end) pList(1:end-1)]; +argList = [argList(end) argList(1:end-1)]; + +for idx = 1:numel(pList) + MS.(argList{idx}) = pList{idx}; +end + +varargout{1} = MS; + + diff --git a/applications/Measurement/NI_DAQ/measurement_chain_elements_calibration_ni.m b/applications/Measurement/NI_DAQ/measurement_chain_elements_calibration_ni.m new file mode 100644 index 0000000000000000000000000000000000000000..f53718a727fcb86c2f09606fda06d09194e7b313 --- /dev/null +++ b/applications/Measurement/NI_DAQ/measurement_chain_elements_calibration_ni.m @@ -0,0 +1,234 @@ +function varargout = measurement_chain_elements_calibration_ni(niSession,MC,element_idx,oldSens,oldReference) +% calibration routine for input elements using the NI hardware (ita_NI_daq_run) +% otherwise identical to ita_measurement_chain_elements_calibration (except ModulITA, Aurelio etc.) + +% Author: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca +% Created: 10-May-2017 + +MCE = MC(1).elements(element_idx); +if ~exist('oldSens','var') + oldSens = MCE.sensitivity_silent; +end + +hw_ch = MC(1).hardware_channel; +preamp_var = false; +switch(lower(MCE.type)) + case {'sensor'} + sInit.Reference = '94'; + sInit.Unit = 'dB re Pa'; + case {'preamp','ad','preamp_robo_fix'} + sInit.Reference = '1'; + sInit.Unit = 'V'; + case ('preamp_var') + preamp_var = true; + case ('none') + varargout(1) = {1}; + return + otherwise + error(['which element type is this - ' lower(MCE.type)]); +end + +if exist('oldReference','var') + sInit.Reference = oldReference; +end + +old_sens_str = [' {old: ' num2str(oldSens) '}']; + +%% GUI +pList = []; + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +ele = numel(pList)+1; +pList{ele}.datatype = 'text'; +pList{ele}.description = ['Calibrating: ' upper(MCE.type) '::' MCE.name '::' 'Hardware Channel: ' num2str(hw_ch) '...']; + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +calibrated_str = ''; +if MCE.calibrated == 0 + calibrated_str = '(UNCALIBRATED)'; +end + +if preamp_var + ele = numel(pList)+1; + pList{ele}.description = 'Ext. Preamp Gain [dB]'; %this text will be shown in the GUI + pList{ele}.helptext = 'External gain of preamp, e.g. BK Type 2610'; %this text should be shown when the mouse moves over the textfield for the description + pList{ele}.datatype = 'int'; %based on this type a different row of elements has to drawn in the GUI + pList{ele}.default = 20*log10(double(MCE.sensitivity_silent)); %default value, could also be empty, otherwise it has to be of the datatype specified above +else + ele = numel(pList)+1; + pList{ele}.datatype = 'text'; + pList{ele}.description = ['Current Sensitivity: ' num2str(MCE.sensitivity) ' ' calibrated_str old_sens_str]; + + ele = numel(pList)+1; + pList{ele}.datatype = 'line'; + + ele = numel(pList)+1; + pList{ele}.description = ['Reference [' num2str(sInit.Unit) ']']; %this text will be shown in the GUI + pList{ele}.helptext = 'Reference voltage or sound pressure level'; %this text should be shown when the mouse moves over the textfield for the description + pList{ele}.datatype = 'int'; %based on this type a different row of elements has to drawn in the GUI + pList{ele}.default = sInit.Reference; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pList)+1; + pList{ele}.description = 'Sampling Rate [Hz]'; %this text will be shown in the GUI + pList{ele}.helptext = 'Sampling Rate'; %this text should be shown when the mouse moves over the textfield for the description + pList{ele}.datatype = 'int'; %based on this type a different row of elements has to drawn in the GUI + pList{ele}.default = ita_preferences('samplingRate'); %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pList)+1; + pList{ele}.description = 'Length [s]'; %this text will be shown in the GUI + pList{ele}.helptext = 'Time to wait for the calibration to be finished'; %this text should be shown when the mouse moves over the textfield for the description + pList{ele}.datatype = 'int'; %based on this type a different row of elements has to drawn in the GUI + pList{ele}.default = 2; %default value, could also be empty, otherwise it has to be of the datatype specified above +end + +ele = numel(pList)+1; +pList{ele}.datatype = 'line'; + +%call GUI +if preamp_var + pList = ita_parametric_GUI(pList,['Ext. Preamp Gain: ' MCE.type '::' MCE.name ' - hwch: ' ... + num2str(hw_ch)]); + if isempty(pList) + varargout{1} = MCE.sensitivity; + else + varargout{1} = 10^(pList{1}/20); + end + return; +else + ele = numel(pList)+1; + pList{ele}.datatype = 'text'; + pList{ele}.description = 'Settings for Calibration Device'; + + pList = ita_parametric_GUI(pList,['Calibration: ' MCE.type '::' MCE.name ' - hwch: ' num2str(hw_ch)],'buttonnames',{'Accept','Calibrate'}); +end + +if isempty(pList) + ita_verbose_info(['Accepting sensitivity for ' MCE.type ' - ' MCE.name ' - hwch: ' num2str(hw_ch)],1) + varargout{1} = MCE.sensitivity; + return; +end + +%% Initialization +%default init struct +oldSens = MCE.sensitivity; +sInit.Channel = hw_ch; +if strcmpi(sInit.Unit,'V') + sInit.Reference = itaValue(pList{1},'V'); +else + sInit.Reference = itaValue(10^(pList{1}/20)*2e-5,'Pa'); +end +sInit.samplingRate = pList{2}; +sInit.length = pList{3}*sInit.samplingRate; + +%% Measurement +% record data +signalRecord = ita_NI_daq_run(sInit.length,niSession,'inputchannels',sInit.Channel,'samplingRate',sInit.samplingRate); +signalRecord = ita_filter_bandpass(signalRecord,'lower',20,'zerophase',false); +signalRecord = signalRecord / MC(1).sensitivity(MCE.type); % compensation of rest + +%% Evaluate +calibrationDomain = ita_preferences('calibrationDomain'); +if ~(strcmpi(calibrationDomain,'time') || strcmpi(calibrationDomain,'frequency')) + error('ita_measurement_chain_elements_calibration: Unknown calibration mode. Please set calibration mode in ita_preferences->ExpertSettings to either ''Time'' or ''Frequency''!'); +end + +validCalibFreqs = ita_preferences('calibrationFrequencies'); +if numel(validCalibFreqs)==0 || any(validCalibFreqs<0) + error('ita_measurement_chain_elements_calibration: No valid calibration Frequency defined. Please set list of valid calibration frequencies in ita_preferences->ExpertSettings!'); +end + +calibFreqTolerance = ita_preferences('calibrationFrequencyTolerance'); +if ~isnumeric(calibFreqTolerance) || ~isfinite(calibFreqTolerance) || calibFreqTolerance<0 + error('ita_measurement_chain_elements_calibration: Invalid calibration frequency tolerance. Please set calibration frequency tolerance in Cent (Default 100) in ita_preferences->ExpertSettings!'); +end + +if strcmpi(calibrationDomain,'time') + nBlocks = 5; + % Signal Segmentation + for idxBlock = 1:nBlocks + signalSegment(idxBlock) = ita_time_crop(signalRecord,[((idxBlock-1)*signalRecord.nSamples/nBlocks+1),(idxBlock*signalRecord.nSamples/nBlocks)],'samples'); %#ok + signalSegment(idxBlock).comment = ['Block number ' num2str(idxBlock)]; %#ok + end + signalSegment = signalSegment.merge; + + % Test single blocks for correct frequency + % in case of low battery or bad signal, the frequency will change + signalTimeWindow = ita_time_window(signalSegment, ... + [round(0.5*signalSegment.nSamples), 1, ... + round(0.5*signalSegment.nSamples)+1, signalSegment.nSamples], ... + 'samples',@hann); % multiply with window before FFT + freq_vec = signalTimeWindow.freqVector; + rmsVals = signalSegment.rms; + % look for maxMagnitude and its frequency -> should be somewhere near + % calibration frequency; if not, don't consider block for sensitivity calculation + [dummy,maxFrequencyIdx] = max(abs(signalTimeWindow.freq),[],1); %#ok + maxFrequency = freq_vec(maxFrequencyIdx); + idxRMSValid = 0; + validBlock = zeros(nBlocks,1); + for idxBlock = 1:nBlocks + ita_verbose_info(['maxFrequency of block number :' num2str(idxBlock) ': ' num2str(maxFrequency(idxBlock)) 'Hz'],1) + validBlock(idxBlock) = any((maxFrequency(idxBlock) > (validCalibFreqs*2.^(-calibFreqTolerance/1200))) & (maxFrequency(idxBlock) < (validCalibFreqs*2.^(calibFreqTolerance/1200)))); + if validBlock(idxBlock)==0 + ita_verbose_info('Invalid frequency detected',0); + else % Calculate rms value for each valid block + idxRMSValid = idxRMSValid+1; + rmsValid(idxRMSValid) = rmsVals(idxBlock); %#ok + end + end + ita_verbose_info(['Number of valid blocks: ' num2str(idxRMSValid)],2) + + % Take Median as most representative value + % check first if enough blocks have been verified as valid + if idxRMSValid == 0 + signalRecord.channelNames{1} = 'CALIBRATION FAILED::INPUT SIGNAL LOOKS LIKE THIS!'; + ita_plot_time(signalRecord); + ita_verbose_info('Signal does not have any of the allowed calibration frequencies.',0); + SensValid = MCE.sensitivity; + elseif idxRMSValid <= round(nBlocks/2) + ita_verbose_info('Measurement might not be accurate. Bad signal.',1); + SensValid = itaValue(median(rmsValid),signalRecord.channelUnits{1})/sInit.Reference; + else + SensValid = itaValue(median(rmsValid),signalRecord.channelUnits{1})/sInit.Reference; + end + +elseif strcmpi(calibrationDomain,'frequency') + str = num2cell(validCalibFreqs); + ok = 0; + while ~ok + [selection,ok] = listdlg('PromptString','Select a calibration frequency:',... + 'SelectionMode','single',... + 'ListString',str); + if ok + freqStart = validCalibFreqs(selection)*2^(-calibFreqTolerance/1200); + freqStop = validCalibFreqs(selection)*2^(+calibFreqTolerance/1200); + + % look for maxMagnitude and its frequency -> should be somewhere near + % the chosen calibration frequency; if not, give out warning + freq_vec = signalRecord.freqVector; + [dummy,maxIdx] = max(abs(signalRecord.freqData),[],1); %#ok + maxFrequency = freq_vec(maxIdx); + ita_verbose_info(['Frequency with highest Amplitude in recorded signal :' num2str(maxFrequency) 'Hz'],1) + if (maxFrequency < freqStart) || (maxFrequency > freqStop) + ita_verbose_info('Oh dear, the signal at you''re calibration frequency is not the strongest component in your signal. I hope you know what you''re doing!',0); + ita_plot_freq(signalRecord) + end + + % Calculate RMS at calibration frequency and consider leakage + calibValues = signalRecord.freq2value(freqStart,freqStop); % make row vector + calibRMS = sqrt(calibValues' * calibValues); + SensValid = itaValue(calibRMS,signalRecord.channelUnits{1})/sInit.Reference; + else + disp('ita_measurement_chain_elements_calibration: Please choose a calibration frequency!'); + end + end +end + +MC(1).elements(element_idx).sensitivity = SensValid; +SensValid = measurement_chain_elements_calibration_ni(niSession,MC,element_idx,oldSens,pList{1}); + +varargout(1) = {SensValid}; +end % function \ No newline at end of file diff --git a/applications/Measurement/NI_DAQ/measurement_chain_output_calibration_ni.m b/applications/Measurement/NI_DAQ/measurement_chain_output_calibration_ni.m new file mode 100644 index 0000000000000000000000000000000000000000..c69a310a8d9b990790def86565db74a6d04aa1af --- /dev/null +++ b/applications/Measurement/NI_DAQ/measurement_chain_output_calibration_ni.m @@ -0,0 +1,190 @@ +function MS = measurement_chain_output_calibration_ni(MS,input_chain_number,ele_idx,old_sens) +% calibration routine for output elements using the NI hardware (ita_NI_daq_run) +% otherwise identical to ita_measurement_chain_output_calibration (except ModulITA, Aurelio etc.) + +% Author: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca +% Created: 10-May-2017 + +if ~exist('old_sens','var') + old_sens_str = ''; +else + old_sens_str = [' {old: ' num2str(old_sens) '}']; +end + +sensFactor = 1; +MC = MS.outputMeasurementChain(1); %always get the latest measurement chain to be on the safe side + +if MC.elements(ele_idx).calibrated ~= -1 + %% GUI + pListExtra = {}; + + MCE = MC.elements(ele_idx); + if any(strfind(lower(MCE.name),'robo')) || any(strfind(lower(MCE.type),'robo')) + default_output2input = 'preamp'; + elseif ismember(MCE.type,{'actuator'}) + default_output2input = 'sensor'; + elseif ismember(MCE.type,{'loudspeaker'}) + default_output2input = 'sensor'; + + ele = numel(pListExtra)+1; + pListExtra{ele}.description = 'Distance to Loudspeaker [m]'; %this text will be shown in the GUI + pListExtra{ele}.helptext = 'distance in meters'; %this text should be shown when the mouse moves over the textfield for the description + pListExtra{ele}.datatype = 'double'; %based on this type a different row of elements has to drawn in the GUI + pListExtra{ele}.default = 1; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pListExtra)+1; + pListExtra{ele}.description = 'Microphone on the floor?'; %this text will be shown in the GUI + pListExtra{ele}.helptext = 'semi-anechoic chamber, microphone on the floor'; %this text should be shown when the mouse moves over the textfield for the description + pListExtra{ele}.datatype = 'bool'; %based on this type a different row of elements has to drawn in the GUI + pListExtra{ele}.default = false; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pListExtra)+1; + pListExtra{ele}.description = 'Window start time[s]'; %this text will be shown in the GUI + pListExtra{ele}.helptext = 'starting time of symmetrical window function'; %this text should be shown when the mouse moves over the textfield for the description + pListExtra{ele}.datatype = 'int'; %based on this type a different row of elements has to drawn in the GUI + pListExtra{ele}.default = 0.05; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pListExtra)+1; + pListExtra{ele}.description = 'Window end time[s]'; %this text will be shown in the GUI + pListExtra{ele}.helptext = 'end time of symmetrical window function'; %this text should be shown when the mouse moves over the textfield for the description + pListExtra{ele}.datatype = 'int'; %based on this type a different row of elements has to drawn in the GUI + pListExtra{ele}.default = 0.1; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pListExtra)+1; + pListExtra{ele}.datatype = 'line'; + else + default_output2input = 'ad'; + end + hw_ch = MC.hardware_channel; + + pList = []; + + ele = numel(pList)+1; + pList{ele}.datatype = 'line'; + + ele = numel(pList)+1; + pList{ele}.datatype = 'text'; + pList{ele}.description = ['Calibrating: ' upper(MCE.type) '::' MCE.name '::' 'Hardware Channel: ' num2str(hw_ch) '...']; + + ele = numel(pList)+1; + pList{ele}.datatype = 'line'; + + calibrated_str = ''; + if MCE.calibrated == 0 + calibrated_str = '(UNCALIBRATED)'; + end + + ele = numel(pList)+1; + pList{ele}.datatype = 'text'; + pList{ele}.description = ['Current Sensitivity: ' num2str(MCE.sensitivity) ' ' calibrated_str old_sens_str]; + + ele = numel(pList)+1; + pList{ele}.datatype = 'line'; + + ele = numel(pList)+1; + pList{ele}.description = 'output2input'; %this text will be shown in the GUI + pList{ele}.helptext = 'Output is connected to this element'; %this text should be shown when the mouse moves over the textfield for the description + pList{ele}.datatype = 'char_popup'; %based on this type a different row of elements has to drawn in the GUI + pList{ele}.list = 'ad|preamp|sensor'; + pList{ele}.default = default_output2input; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pList)+1; + pList{ele}.description = 'outputamplification'; %this text will be shown in the GUI + pList{ele}.helptext = 'in dBFS'; %this text should be shown when the mouse moves over the textfield for the description + pList{ele}.datatype = 'char'; %based on this type a different row of elements has to drawn in the GUI + pList{ele}.default = MS.outputamplification; %default value, could also be empty, otherwise it has to be of the datatype specified above + + ele = numel(pList)+1; + pList{ele}.datatype = 'line'; + + pList = [pList pListExtra]; + + %call GUI + pList = ita_parametric_GUI(pList,['Calibration: ' MCE.type '::' MCE.name ' - hwch: ' num2str(hw_ch)],'buttonnames',{'Accept','Calibrate'}); + + if isempty(pList) + ita_verbose_info(['Accepting sensitivity for ' MCE.type ' - ' MCE.name ' - hwch: ' num2str(hw_ch)],1) + MC.elements(ele_idx).sensitivity = MCE.sensitivity; %set sensitivity + MS.outputMeasurementChain = MC; + return; + else + output2input = pList{1}; + MS.outputamplification = pList{2}; + + %% measurement + %try to get the sensitivity of the chain. modulita and robo could + %be uninitialized + old_sens = MCE.sensitivity; + + MS.inputChannels = MS.inputMeasurementChain(input_chain_number).hardware_channel; + %TODO: Create latency vector + if MS.latencysamples == 0 + MS.run_latency; + end + % for electrical measurements, get best SNR with autoranging + if ~strcmpi(output2input,'sensor') + MS.run_autorange(0,pList{2}); + end + + inputChannels = MS.inputChannels; + outputChannels = MS.outputChannels; + samplingRate = MS.samplingRate; + final_excitation = MS.final_excitation; + latencysamples = MS.latencysamples; + + % measure TF + a = ita_NI_daq_run(final_excitation,MS.niSession,'InputChannels',inputChannels, ... + 'OutputChannels', outputChannels,'latencysamples',latencysamples,'samplingRate',samplingRate); + + a = a * MS.compensation / MS.outputamplification_lin; + a.signalType = 'energy'; + + if ~isempty(pListExtra) + %% compensation of distance + distance = itaValue(pList{3},'m'); + travel_time = distance / ita_constants('c'); + + a = ita_time_shift(a,-double(travel_time),'time'); + a = a * distance; + + %% floor compensation + if pList{4} + a = ita_amplify(a,'-6dB'); + end + + %% time windowing + a = ita_time_window(a,[pList{5} pList{6}],'time','symmetric'); + + end + + %% get FRF up to this point + frf_upto = MC.response(lower(MC.elements(ele_idx).type)); + if ~isempty(frf_upto) + a = a*ita_invert_spk_regularization(frf_upto,[1 MS.samplingRate/2],'filter'); + end + + %% get sensitivity of element around 1kHz + value = itaValue ( mean(abs(a.freq2value(950:1050))) , a.channelUnits{1}); + MC.elements(ele_idx).response = a / value; + switch lower(output2input) + case {'ad'} + value = value / sensFactor / MS.inputMeasurementChain(input_chain_number).sensitivity('preamp') ... + / MS.outputMeasurementChain.sensitivity(MCE.type); + + case {'preamp'} + value = value / sensFactor / MS.inputMeasurementChain(input_chain_number).sensitivity('sensor')... + / MS.outputMeasurementChain.sensitivity(MCE.type); + + case {'sensor'} + value = value / sensFactor / MS.inputMeasurementChain(input_chain_number).sensitivity()... + / MS.outputMeasurementChain.sensitivity(MCE.type); + otherwise + error('element type unknown') + end + + MC.elements(ele_idx).sensitivity = value; + MS.outputMeasurementChain = MC; + end + MS = measurement_chain_output_calibration_ni(MS,input_chain_number,ele_idx,old_sens); +end +end % function \ No newline at end of file diff --git a/applications/RoomAcoustics/BuildingAcoustics/ita_soundInsulationIndexAirborne.m b/applications/RoomAcoustics/BuildingAcoustics/ita_soundInsulationIndexAirborne.m index 548c9fd95d50cede827d16d5a90b407a905450a9..9327305586f80d254763fc3b9f5af8ed5d458fcd 100644 --- a/applications/RoomAcoustics/BuildingAcoustics/ita_soundInsulationIndexAirborne.m +++ b/applications/RoomAcoustics/BuildingAcoustics/ita_soundInsulationIndexAirborne.m @@ -8,91 +8,89 @@ function varargout = ita_soundInsulationIndexAirborne(varargin) %% input parsing -sArgs = struct('pos1_data','anything','bandsperoctave',3,'freqVector',[],'createPlot',false); +sArgs = struct('pos1_data','anything','bandsperoctave',3,'freqVector',[],'createPlot',false,'type','ISO'); [data,sArgs] = ita_parse_arguments(sArgs,varargin); -%% additional parameters -if sArgs.bandsperoctave == 3 - refSurf = 32; -elseif sArgs.bandsperoctave == 1 - refSurf = 10; -else - error([upper(mfilename) ':wrong input for badnsperoctave']); -end -refFreq = 500; - %% reference curves -if sArgs.bandsperoctave == 1 - refCurve = [36 45 52 55 56]; % Reference curve according to ISO 717-1 - freq = [125 250 500 1000 2000]; % Frequencies according to ISO 717-1 - lFreq = length(refCurve); +if strcmpi(sArgs.type,'iso') % Reference curve and frequencies according to ISO 717-1 + outputStr = 'R_W'; + roundingFactor = 0.1; + deficiencyLimit = Inf; + if sArgs.bandsperoctave == 1 + refCurve = [36 45 52 55 56]-52; + freq = [125 250 500 1000 2000]; + refSurf = 10; + elseif sArgs.bandsperoctave == 3 + freq = [100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150]; + refCurve = [33 36 39 42 45 48 51 52 53 54 55 56 56 56 56 56]-52; + refSurf = 32; + else + error([upper(mfilename) ':wrong input for bandsperoctave']); + end +elseif strcmpi(sArgs.type,'astm') % Reference curve and frequencies according to ASTM E413 + outputStr = 'STC'; + roundingFactor = 1; + deficiencyLimit = 8; + sArgs.bandsperoctave = 3; + freq = [125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000]; + refCurve = [-16 -13 -10 -7 -4 -1 0 1 2 3 4 4 4 4 4 4]; + refSurf = 32; else - freq = [100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150]; % Frequencies according to ISO 717-1 - refCurve = [33 36 39 42 45 48 51 52 53 54 55 56 56 56 56 56]; % Reference curve according to ISO 717-1 - lFreq = length(refCurve); + error([upper(mfilename) ':wrong input for type']); end +freq = freq(:); +refCurve = refCurve(:); +dbStep = 1; %% prepare data -msgExtrap = 'Sound insulation data will be extrapolated'; -if isa(data,'itaSuper') - freqVector = data.freqVector; - soundInsulation = 20.*log10(interp1(freqVector,data.freqData(:,1),freq,'spline','extrap')); -else +if ~isa(data,'itaSuper') freqVector = sArgs.freqVector; if ~isempty(freqVector) - soundInsulation = interp1(freqVector,data,freq,'spline','extrap'); % not in dB?!? + data = itaResult(10.^(data(:)./20),freqVector(:),'freq'); else error([upper(mfilename) ':not enough input data']); end end +freqVector = data.freqVector; +soundInsulation = 20.*log10(interp1(freqVector,data.freqData(:,1),freq,'spline','extrap')); -if isempty(find(freqVector <= freq(1),1,'first')) || isempty(find(freqVector >= freq(end),1,'first')) - warning(upper(mfilename),msgExtrap); +if min(freqVector) > freq(1) || max(freqVector) < freq(end) + warning([upper(mfilename) ': Sound insulation data will be extrapolated']); end %% sound insulation index -soundInsulation = round(soundInsulation*10)/10; -delta = refCurve-soundInsulation; -soundInsulationIndexTest = sum(delta(delta>0)); +soundInsulation = round(soundInsulation/roundingFactor)*roundingFactor; +soundInsulationIndex = min(floor(soundInsulation-refCurve)); +delta = max(0,refCurve + soundInsulationIndex - soundInsulation); counter = 0; % stopping criterion -% shift reference curve until 32dB is reached +% shift reference curve until limits are reached -Diff = 0; -while abs(soundInsulationIndexTest -refSurf)>1 && counter < 1e6 - if counter == 0 % Anpassung der refTerzkurve - Diff = round(mean(delta)*10)/10; - elseif sum(soundInsulationIndexTest) < refSurf-1 - Diff = Diff+0.1; - else - Diff = Diff-0.1; - end - - delta = (refCurve+Diff)- soundInsulation; - soundInsulationIndexTest = sum(delta(delta>0)); - +while sum(delta) < refSurf && all(delta <= deficiencyLimit) && counter < 1e3 + soundInsulationIndex = soundInsulationIndex + dbStep; + delta = max(0,refCurve + soundInsulationIndex - soundInsulation); counter = counter+1; end - -soundInsulationIndex = round((refCurve(freq == refFreq)+Diff)*100)/100; +soundInsulationIndex = soundInsulationIndex - dbStep; +delta = max(0,refCurve + soundInsulationIndex - soundInsulation); +deficiencies = itaResult(delta,freq,'freq')*itaValue(1,'dB'); +deficiencies.allowDBPlot = false; %% output if sArgs.createPlot - plotResult = itaResult; - plotResult.freqVector = freq; - plotResult.freqData(:,1) = 10.^((refCurve+Diff)./20); - plotResult.freqData(:,3) = 10.^(soundInsulation./20); - plotResult.freqData(:,2) = ones(lFreq,1)*10.^(soundInsulationIndex./20); - - plotResult.channelNames{3} = 'sound insulation'; - plotResult.channelNames{1} = 'shifted reference curve'; - plotResult.channelNames{2} = ['R_W = ' num2str(soundInsulationIndex) 'dB']; - - plotResult.plot_freq; - xlim([min(freq) max(freq)]); + fgh = ita_plot_freq(data); + plotResult = itaResult([10.^((refCurve+soundInsulationIndex)./20), [ones(sum(freq<=500),1)*10.^(soundInsulationIndex./20); nan(sum(freq>500),1)]],freq,'freq'); + ita_plot_freq(plotResult,'figure_handle',fgh,'axes_handle',gca,'hold'); + bar(gca,deficiencies.freqVector,deficiencies.freq,'hist'); + [maxDef,maxIdx] = max(deficiencies.freq); + legend({'Sound transmission loss','Shifted reference curve',[outputStr ' = ' num2str(soundInsulationIndex) 'dB'],['Deficiencies (sum: ' num2str(sum(deficiencies.freq)) 'dB, max: ' num2str(maxDef) 'dB at ' num2str(deficiencies.freqVector(maxIdx)) 'Hz)']}); + ylim([0 max(max(soundInsulation),max(refCurve)+soundInsulationIndex)+15]); end varargout{1} = soundInsulationIndex; % reference curve specified at the freqVector specified by the input itaResult: -if nargout == 2 - varargout{2} = interp1(freq, 10.^((refCurve+Diff)./20), freqVector, 'linear'); +if nargout >= 2 + varargout{2} = interp1(freq, 10.^((refCurve+soundInsulationIndex)./20), freqVector, 'linear'); + if nargout == 3 + varargout{3} = deficiencies; + end end diff --git a/applications/RoomAcoustics/ita_speech_transmission_index.m b/applications/RoomAcoustics/ita_speech_transmission_index.m index a9420531dec879f7a839d45e270a119deff4ef76..6de82efde1110a7ac71595a167e1ced93520f3ba 100644 --- a/applications/RoomAcoustics/ita_speech_transmission_index.m +++ b/applications/RoomAcoustics/ita_speech_transmission_index.m @@ -35,6 +35,16 @@ function varargout = ita_speech_transmission_index(varargin) sArgs = struct('pos1_ir','itaAudio', 'levels', [], 'SNR', [],'plot',false,'analytic', false, 'gender', 'male'); [ir,sArgs] = ita_parse_arguments(sArgs,varargin); +if ~strcmpi(ir.signalType,'energy') + ita_verbose_info('Your IR does not have the correct signalType, I will fix this, but be careful!',0); + ir.signalType = 'energy'; +end + +if ir.trackLength < 1.6 + ita_verbose_info('IR is shorter than ISO 60268-16 recommends, I hope you know what you are doing! I will extend the data',0); + ir = ita_extend_dat(ir,1.6); +end + if strcmpi(sArgs.gender,'male') genderIndex = 1; elseif strcmpi(sArgs.gender,'female') @@ -45,7 +55,7 @@ end %% constants fk = ita_ANSI_center_frequencies([125 8000],1); -fm = ita_ANSI_center_frequencies([62 1250],3)./100; +fm = [0.63 0.8 1 1.25 1.6 2 2.5 3.15 4 5 6.3 8 10 12.5]; fm(5) = 1.6; % rounding errors I_k_rt = 10.^([46 27 12 6.5 7.5 8 12].'./10); @@ -77,7 +87,15 @@ I_k = 10.^(L(:)./10) + 10.^((L(:)-SNR(:))./10); h_k = ita_filter_fractional_octavebands(ir,'bandsperoctave',1,'freqRange',[125 8000]); h_k_sq = h_k.^2; % modulation transfer function values -m_k_fm = bsxfun(@rdivide,abs(h_k_sq.freq2value(fm)),abs(h_k_sq.freq2value(0)).*(1+10.^(-SNR./10))).'; +m_k_fm = zeros(numel(fk),numel(fm)); +for iM = 1:numel(fm) + % to get an FFT bin exactly at fm + newLength = floor(h_k.trackLength*fm(iM))/fm(iM); + h_k_sq_tmp = ita_time_crop(h_k_sq,[0 newLength],'time'); + m_k_fm(:,iM) = (abs(h_k_sq_tmp.freq2value(fm(iM)))./(abs(h_k_sq_tmp.freq2value(0)).*(1+10.^(-SNR./10)))).'; +end +% old version: +% m_k_fm = bsxfun(@rdivide,abs(h_k_sq.freq2value(fm)),abs(h_k_sq.freq2value(0)).*(1+10.^(-SNR./10))).'; % correction terms for masking and reception threshold I_k_am = zeros(numel(fk),1); @@ -120,7 +138,7 @@ if sArgs.analytic % RT for analytic result RT = ita_roomacoustics(ir,'T30','freqRange',[125 8000],'bandsperoctave',1); RT = RT.T30.freq; - m_analytic = 1./sqrt(1 + (2*pi.*bsxfun(@times,fm,RT)./13.8).^2); + m_analytic = 1./(bsxfun(@times,sqrt(1 + (2*pi.*bsxfun(@times,fm,RT)./13.8).^2),(1+10.^(-SNR.'./10)))); m_analytic = min(bsxfun(@times,m_analytic,correctionTerm),1); SNR_eff_analytic = min(max(real(10.*log10(m_analytic./(1-m_analytic))),-15),15); TI_analytic = (SNR_eff_analytic + 15)./30; @@ -144,7 +162,7 @@ varargout(1) = {STI}; if nargout > 1 varargout(2) = {MTI}; if nargout > 2 && sArgs.analytic - varargout(3) = STI_analytic; + varargout(3) = {STI_analytic}; end end diff --git a/applications/SpatialAudio/SourceTranslation/m/calculateHandT_slideSource.m b/applications/SpatialAudio/SourceTranslation/m/calculateHandT_slideSource.m index 005779c9122c3ffa139dd223d0b928bf41489abe..6a3df27e401f4464f0e90c835cf91897669a1613 100644 --- a/applications/SpatialAudio/SourceTranslation/m/calculateHandT_slideSource.m +++ b/applications/SpatialAudio/SourceTranslation/m/calculateHandT_slideSource.m @@ -47,7 +47,7 @@ coordinate.cart = loc; grid = params.display.grid + coordinate; % Spherical Harmonics: num_angels x (Nc+1)^2. -Ynm = ita_sph_base(grid, N, 'Williams', true); +Ynm = ita_sph_base(grid, N); % prepare hankel functions for translated source for k_ind = 1 : K diff --git a/applications/SpatialAudio/SourceTranslation/m/initParams.m b/applications/SpatialAudio/SourceTranslation/m/initParams.m index a67129e7921e151188d9b568fc9a43f532b76226..ced435157acbb5c8fd9c4a4d8711c6382bc89ae3 100644 --- a/applications/SpatialAudio/SourceTranslation/m/initParams.m +++ b/applications/SpatialAudio/SourceTranslation/m/initParams.m @@ -168,8 +168,7 @@ params.display.grid.r = params.array.r; % Spherical Harmonics: num_angels x (Nc+1)^2. params.display.Ynm = ita_sph_base(params.display.grid,... - params.source.N,... - 'Williams', true); + params.source.N); % Build up a reverse grid. th = params.display.grid.theta; @@ -207,8 +206,7 @@ params.array.mics = params.array.grid.nPoints; % Number of microphones % Spherical Harmonics: num_angels x (Narray+1)^2. params.array.Ynm = ita_sph_base(params.array.grid,... - params.array.N,... - 'Williams', true); + params.array.N); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Parameters that are relative to the center diff --git a/applications/SpatialAudio/SourceTranslation/m/simulateHandT_firstTranslation.m b/applications/SpatialAudio/SourceTranslation/m/simulateHandT_firstTranslation.m index bf3f6788b42f274ea39dc16bb42847aaebea5fd0..99943b33f6e1512c879939e75a5aa9a308407b51 100644 --- a/applications/SpatialAudio/SourceTranslation/m/simulateHandT_firstTranslation.m +++ b/applications/SpatialAudio/SourceTranslation/m/simulateHandT_firstTranslation.m @@ -53,7 +53,7 @@ for k_ind = 1 : K grid = params.display.grid - coordinate; % Spherical Harmonics: num_angels x (Nc+1)^2. - Ynm = ita_sph_base(grid, N, 'Williams', true); + Ynm = ita_sph_base(grid, N); hn_disPoints_x_N = ita_sph_besselh([0:N], 1,... kr_disPoints_x_freqs(:,k_ind)); diff --git a/applications/SpatialAudio/ita_decodeAmbisonics.m b/applications/SpatialAudio/ita_decodeAmbisonics.m index 70b32e61a2cd9346b69320ef1c70e946a973a1ef..1834509106356a63362c28a2e56865805a6c549a 100644 --- a/applications/SpatialAudio/ita_decodeAmbisonics.m +++ b/applications/SpatialAudio/ita_decodeAmbisonics.m @@ -57,7 +57,7 @@ for k=1:numel(weights) Bformat(:,k)=weights(k).*Bformat(:,k); %Apply weighting factors end % SH and inversion -Y = ita_sph_base(LS,TheOrder,'williams',false); % generate basefunctions +Y = ita_sph_base(LS, TheOrder, 'real'); % generate basefunctions Yinv=pinv(Y); % calculate Pseudoinverse diff --git a/applications/SpatialAudio/ita_encodeAmbisonics.m b/applications/SpatialAudio/ita_encodeAmbisonics.m index 1751b3a203ce6c81a8b60cd7aa9b248ec93793b3..78fdf156a09dab45552d85d074d6be61e6bf87fc 100644 --- a/applications/SpatialAudio/ita_encodeAmbisonics.m +++ b/applications/SpatialAudio/ita_encodeAmbisonics.m @@ -17,7 +17,7 @@ opts.audioSource=''; %Signal of the source (itaAudio or filename) opts = ita_parse_arguments(opts, varargin); % Encode Source -Bformat = ita_sph_base(sourcePos, opts.order, 'Williams', false); +Bformat = ita_sph_base(sourcePos, opts.order, 'real'); if ~isempty(opts.audioSource) if isa(opts.audioSource,'char') diff --git a/applications/SpatialAudio/ita_panVBAP.m b/applications/SpatialAudio/ita_panVBAP.m new file mode 100644 index 0000000000000000000000000000000000000000..9f816483d0aeb15c70d6431d2ed619ae835a72e3 --- /dev/null +++ b/applications/SpatialAudio/ita_panVBAP.m @@ -0,0 +1,56 @@ +function weights = panVBAP(pos_LS,pos_VS) +%panVBAP - Calculate weights for VBAP +% +% This function receives the position of the loudspeakers and the position +% of the virtual source. Both input must be given as objects of the class +% itaCoordinates. +% +% The output is the set of frequency independent weights used to pan the +% virtual source on the given array. +% +% Call: weights = panVBAP(pos_LS,pos_VS) +% +% Author: Michael Kohnen -- Email: mko@akustik.rwth-aachen.de +% Former author: Bruno Masiero -- Email: bma@akustik.rwth-aachen.de +% Created: 13-Jun-2011 +% Last modified: 04-May-2017 +%$ENDHELP$ + +%% Preliminary tests +% if pos_LS.isPlane +% Number_of_active_loudspeakers = 2; %Extended Stereo +% else + Number_of_active_loudspeakers = 3; +% end + +weights = zeros(pos_VS.nPoints,pos_LS.nPoints); + +if pos_VS.r == 0 + error('No direction for Virtual Source was given!') +end + +%% Find the closest loudspeakers +% Calculate the distance of each loudspeaker to the virtual source with the +% help of the itaCoordinate overloaded function itaCoordinate.r. +% To sort the distance in ascending value, use the function sort. + aux = pos_LS - pos_VS; + dist = aux.r; + [junk,index] = sort(dist,'ascend'); + index = index(1:Number_of_active_loudspeakers); + + active_loudspeakers = pos_LS.n(index); + +%% Calculate the weights for the active loudspeakers +% Create a base matrix with the direction of the active loudspeakers and +% multiply the direction of the virtual source with the inverse of this +% matrix. +% Don't forget to normalize yor results with C = 1; +for idx = 1:pos_VS.nPoints + p = pos_VS.n(idx).cart; + L = active_loudspeakers.cart; + g = p*pinv(L); + % Re-normalize. + g = abs(g)/norm(g); + weights(idx,index) = g; +end + diff --git a/applications/SphericalHarmonics/itaSamplingSphReal.m b/applications/SphericalHarmonics/itaSamplingSphReal.m index f0c4b6cbd702c5e6b23503c78cd7c2a57cd165fb..2ffc9ddd7e37d55d4f6ff4d567c3014157edca0a 100644 --- a/applications/SphericalHarmonics/itaSamplingSphReal.m +++ b/applications/SphericalHarmonics/itaSamplingSphReal.m @@ -14,7 +14,7 @@ classdef itaSamplingSphReal < itaSamplingSph end methods(Access = protected) function this = set_base(this, nmax) - this.Y = ita_sph_base(this,nmax,'Williams',false); + this.Y = ita_sph_base(this, nmax, 'real'); ita_verbose_info(['Real valued spherical harmonics up to order ' num2str(nmax) ' calculated.']); end end diff --git a/applications/SphericalHarmonics/ita_sph_base.m b/applications/SphericalHarmonics/ita_sph_base.m index 5d5a14be74b19bd008016114d594c3876bda9561..ddaaa5b139a592819f5c8706c33bb11d940be513 100644 --- a/applications/SphericalHarmonics/ita_sph_base.m +++ b/applications/SphericalHarmonics/ita_sph_base.m @@ -1,83 +1,79 @@ -function Y = ita_sph_base(s, nmax, type, complex) +function Y = ita_sph_base(varargin) %ITA_SPH_BASE - creates spherical harmonics (SH) base functions -% function Y = ita_sph_base(s, nmax, type) -% -% calculates matrix with spherical harmonic base functions +% function Y = ita_sph_base(sampling, Nmax, options) +% Y is a matrix with dimensions [nr_points x nr_coefs] +% calculates matrix with spherical harmonic basis functions % for the grid given in theta and phi -% give type of normalization +% % % the definition was taken from: % E. G. Williams, "Fourier Acoustics", -% Academic Press, San Diego, 1999. p.??? +% Academic Press, San Diego, 1999. p. 190 % -% Martin Pollow (mpo@akustik.rwth-aachen.de) -% Institute of Technical Acoustics, RWTH Aachen, Germany -% 03.09.2008 +% This definition includes the Condon-Shotley phase term (-1)^m +% +% Syntax: +% Y = ita_sph_base(sampling, Nmax, options) +% +% Options (default): +% 'norm' ('orthonormal') : Normalization type +% 'real' (false) : Return real valued SH +% +% Example: +% Y = ita_sph_base(sampling, Nmax, 'norm', 'Williams') +% +% See also: +% ita_toolbox_gui, ita_read, ita_write, ita_generate +% +% Reference page in Help browser +% doc ita_sph_base % -% Deleted and new implementation -% (careful, history seems to be lost, but checkout of old version is possible): -% Bruno Masiero (bma@akustik.rwth-aachen.de) -% Institute of Technical Acoustics, RWTH Aachen, Germany -% 12.04.2011 - % -% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved. -% You can find the license for this m-file in the application 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. % +% +% +% Martin Pollow (mpo@akustik.rwth-aachen.de) +% Institute of Technical Acoustics, RWTH Aachen, Germany +% 03.09.2008 -% check input -if nargin < 4 - % Use complex version of the spherical harmonics - complex = true; -end -if nargin < 3 - type = 'Williams'; - % disp(['using default normalization: ' type]); -end +sArgs = struct('pos1_sampling', 'itaCoordinates', ... + 'pos2_Nmax', 'int', ... + 'norm', 'orthonormal', ... + 'real', false); +[sampling, Nmax, sArgs] = ita_parse_arguments(sArgs, varargin); %check for invalid sampling -if s.nPoints < 1 +if sampling.nPoints < 1 Y = []; + ita_verbose_info('The sampling needs to consist of at least one point.', 0) return end % vectorize grid angles -theta = s.theta; -phi = s.phi; - -nr_points = numel(theta); -if nr_points ~= numel(phi) - error('theta and phi must have same number of points'); -end - -nr_coefs = (nmax+1)^2; -nm = 1:nr_coefs; +theta = sampling.theta; +phi = sampling.phi; -% avoid "~" for backward compatibility -% [~,m] = ita_sph_linear2degreeorder(nm); -[n,m] = ita_sph_linear2degreeorder(nm); %#ok +nm = 1:(Nmax+1)^2; - -%% Check if either real or complex bases wanted - -%% Complex SH base -% [nr_points nr_coefs] +[~,m] = ita_sph_linear2degreeorder(nm); exp_term = exp(1i*phi*m); % calculate a matrix containing the associated legendre functions % function Pnm = ass_legendre_func % calculate the matrix of associated Legrendre functions -Pnm = zeros(nr_points, nr_coefs); +Pnm = zeros(sampling.nPoints, (Nmax+1)^2); -for ind = 0:nmax +for ind = 0:Nmax % define the linear indices for the i'th degree index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order index_m_pos = ita_sph_degreeorder2linear(ind,1:ind); index_m_pos0 = ita_sph_degreeorder2linear(ind,0:ind); % define positive orders with correct normalization - switch lower(type) + switch lower(sArgs.norm) case {'orthonormal','williams'} % the Pnm's used here are the Pnm's from Williams multiplied with the % orthonormality factor sqrt((2n+1)./(4*pi).*(n-m)! ./ (n+m)!) @@ -89,13 +85,20 @@ for ind = 0:nmax bsxfun(@times,(-1).^(0:ind)*sqrt(2),legendre(ind,cos(theta.'),'norm').'); case {'schmidt','sch','semi-normalized'} - Pnm(:,index_m_pos0) = legendre(ind,cos(theta.'),'sch').'; + Pnm(:,index_m_pos0) = ... + bsxfun(@times,(-1).^(0:ind),legendre(ind,cos(theta.'),'sch').'); + case {'ambix'} + Pnm(:,index_m_pos0) = ... + bsxfun(@times,(-1).^(0:ind)/sqrt(8*pi),legendre(ind,cos(theta.'),'sch').'); otherwise error('Wow! I do not know this normalization!') end % copy the Pnm data to the left side of the Toblerone spectrum + % the phase term over theta is not included up to this point and Pnm is + % already normalized. We do not need to use the complex conjugate of + % Pnm or renormalize. if ind > 0 Pnm(:,index_m_neg) = bsxfun(@times,(-1).^(1:ind),Pnm(:,index_m_pos)); end @@ -104,22 +107,12 @@ end % compose the spherical harmonic base functions Y = Pnm .* exp_term; -if ~complex - % now the conversion is done outside of this function - Y = ita_sph_complex2real(Y')'; - - % below the old code: -% for ind = 1:nmax -% % define the linear indices for the i'th degree -% index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order -% index_m_pos = ita_sph_degreeorder2linear(ind,1:ind); -% -% for m = 1:length(index_m_neg) -% C = ((-1)^m*Y(:,index_m_pos(m)) + Y(:,index_m_neg(m)))/sqrt(2); -% S = ((-1)^m*Y(:,index_m_neg(m)) - Y(:,index_m_pos(m)))/sqrt(2)/1i; -% -% Y(:,ita_sph_degreeorder2linear(ind,m)) = C; -% Y(:,ita_sph_degreeorder2linear(ind,-m)) = S; -% end -% end +if strcmp(sArgs.norm, 'ambix') + % multiply by sqrt(2-delta_m) + mask = ~(m | zeros(size(m))); + Y(:,mask) = Y(:,mask) * sqrt(2); +end + +if sArgs.real + Y = ita_sph_complex2real(Y.').'; end \ No newline at end of file diff --git a/applications/SphericalHarmonics/ita_sph_complex2real.m b/applications/SphericalHarmonics/ita_sph_complex2real.m index 54261416dc4b526c36c12b900fa4d0996d936616..e4ee1d1942f79ae4c9f912c803fe88b2d0ec62cd 100644 --- a/applications/SphericalHarmonics/ita_sph_complex2real.m +++ b/applications/SphericalHarmonics/ita_sph_complex2real.m @@ -1,98 +1,74 @@ -function SH = ita_sph_complex2real(SH) -% converts a complex valued base or SH vector into its real base / SH vectors +function SH = ita_sph_complex2real(varargin) +%ITA_SPH_COMPLEX2REAL - Transforms from complex to valued basis functions +% This function transfroms spherical harmonic coefficients with complex +% valued basis functions to their corresponding real valued +% representation. The sign convention of the real valued basis functions +% can be chosen. The definition in 'Williams - Fourier Acoustics' is +% assumed for the complex basis functions (including the Condon-Shotley +% phase). The phase conventions for the real valued spherical harmonics +% can be 'raven' (which equals the ambix phase convention as defined in +% "Nachbar - AMBIX: A Suggested Ambisonics Format (revised by Zotter)") +% or 'zotter' as defined in "Zotter - Analysis and synthesis of +% sound-radiation with spherical arrays" % +% Syntax: +% SH_real = ita_sph_complex2real(SH_cplx, options) +% +% Options (default): +% 'phase' ('ambix') : Phase definitions as in the AmbiX format +% +% Example: +% SH_real = ita_sph_complex2real(SH_cplx) +% +% See also: +% ita_toolbox_gui, ita_read, ita_write, ita_generate +% +% Reference page in Help browser +% doc ita_sph_complex2real % -% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved. -% You can find the license for this m-file in the application 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. % -% SH can be the vector/matrix of coefficients or the maximum order to -% calculate a matrix. - -if isnumeric(SH) - if numel(SH) == 1 - % return a matrix - SH = ita_sph_complex2real_mat(SH); - else - % assume this as SH coefs - SH = ita_sph_complex2real_coefs(SH); - end -else - error('please check syntax') -end - -end - -function SH = ita_sph_complex2real_coefs(SH) - -% assume the SH dimension as first dimension -sizeY = size(SH); +% Rewrite, original author unknown +% +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: 12-Jun-2017 -% check for max 2 dimension -if numel(sizeY) > 2 - reshape(SH, sizeY(1), []); -end +sArgs = struct('pos1_SH', 'double', ... + 'phase', 'ambix'); +[SH, sArgs] = ita_parse_arguments(sArgs, varargin); -% check for SH in 1st dimension -if sizeY(1) == 1 && sizeY(2) > 1 - SH = SH.'; - isTransposed = true; -else - isTransposed = false; +if ~isnatural(sqrt(size(SH, 1)) - 1) + ita_verbose_info('The dimensions of the input data do not match.', 0) + SH = []; + return end -nmax = sqrt(size(SH,1)) - 1; -if nmax ~= round(nmax) - error('ita_sph_complex2real something wrong with the number of coefficients (they do not fill up to a full order)'); -end +Nmax = floor(sqrt(size(SH,1))-1); -for ind = 1:nmax +for ind = 1:Nmax % define the linear indices for the i'th degree index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order index_m_pos = ita_sph_degreeorder2linear(ind,1:ind); - + for m = 1:length(index_m_neg) cPos = SH(index_m_pos(m),:); cNeg = SH(index_m_neg(m),:); - + rPos = ((-1)^m * cPos + cNeg) / sqrt(2); - rNeg = ((-1)^m * cNeg - cPos) / sqrt(2) .* 1i; - + switch sArgs.phase + case {'ambix', 'raven'} + rNeg = (cNeg - (-1)^m*cPos) / sqrt(2) .* 1i; + case 'zotter' + rNeg = (-cNeg + (-1)^m * cPos) / sqrt(2) .* 1i; + otherwise + ita_verbose_info('I do not know this phase convention. Aborting...', 0); + SH = []; + return + end SH(ita_sph_degreeorder2linear(ind,m),:) = rPos; SH(ita_sph_degreeorder2linear(ind,-m),:) = rNeg; end -end - -if isTransposed - SH = SH.'; -end - -% now bring to the old shape -SH = reshape(SH,sizeY); -end - -function T = ita_sph_complex2real_mat(nmax) - -if numel(nmax) > 1 - error; -end - -nSH = (nmax+1).^2; - -lin = (1:nSH).'; -multiplicator = 2*mod(lin,2) - 1; - -[deg,order] = ita_sph_linear2degreeorder(lin); -linNeg = ita_sph_degreeorder2linear(deg,-order); - -isPos = order > 0; -isNeg = order < 0; - -T = sparse(eye(nSH)); -T_orig = T; -T_swap = T .* diag(multiplicator); - -T(lin(isPos),:) = (T_swap(lin(isPos),:) + T_orig(linNeg(isPos),:)) / sqrt(2); -T(lin(isNeg),:) = (T_swap(lin(isNeg),:) - T_orig(linNeg(isNeg),:)) / sqrt(2) .* 1i; -end +end \ No newline at end of file diff --git a/applications/SphericalHarmonics/ita_sph_mimo_error_simulation.m b/applications/SphericalHarmonics/ita_sph_mimo_error_simulation.m index 0af3103f5f82f33a06893a384973bd698ba32b03..cb67b9949680e51e35312337b5c262b42a8054ec 100644 --- a/applications/SphericalHarmonics/ita_sph_mimo_error_simulation.m +++ b/applications/SphericalHarmonics/ita_sph_mimo_error_simulation.m @@ -426,13 +426,12 @@ sArgs = struct('pos1_receiverCoords','itaCoordinates',... 'pos5_k','double',... 'r',[],... 'r_eq',1,... - 'norm',false,... - 'shType','complex'); + 'norm',false); [receiverCoords,receiverNmax,sourceCoords,sourceNmax,k,sArgs] = ita_parse_arguments(sArgs,varargin); [receiverLookDirection, sourceLookDirection,r] = array_orientation(receiverCoords,sourceCoords); -yReceiver = ita_sph_base(receiverLookDirection,receiverNmax,'Williams',strcmp(sArgs.shType,'complex'))'; -ySource = ita_sph_base(sourceLookDirection,sourceNmax,'Williams',strcmp(sArgs.shType,'complex'))'; +yReceiver = ita_sph_base(receiverLookDirection,receiverNmax)'; +ySource = ita_sph_base(sourceLookDirection,sourceNmax)'; % avoid sArgs in parfor since it is a broadcast variable if ~isempty(sArgs.r) diff --git a/applications/SphericalHarmonics/ita_sph_modal_strength.m b/applications/SphericalHarmonics/ita_sph_modal_strength.m index 0b5502e083b08c134acbd9b47fa6b9ca7df645f1..3b891b0702882949b2347c82709929b5b820e1f7 100644 --- a/applications/SphericalHarmonics/ita_sph_modal_strength.m +++ b/applications/SphericalHarmonics/ita_sph_modal_strength.m @@ -98,7 +98,7 @@ for idxRad = 1:numel(uniqueRad) end case {'loudspeaker','ls'} if strcmp(type,'rigid') && isempty(sArgs.scatterer) - bn = diag(-1i.^n*Z_air)*bn.*ita_sph_besselh(n,sArgs.hankelKind,sArgs.dist*kVec); + bn = diag((-1i) * (-1)^sArgs.hankelKind * Z_air)*bn.*ita_sph_besselh(n,sArgs.hankelKind,sArgs.dist*kVec); else ita_verbose_info('This design type makes no sense for a loudspeaker array.',0); varargout{1} = []; diff --git a/applications/SphericalHarmonics/ita_sph_real2complex.m b/applications/SphericalHarmonics/ita_sph_real2complex.m index 5a457a2bd86a5cb35c33e76e72edebf2c6ef97ee..684c3756693ecb9280a087ff10c6de74202187e4 100644 --- a/applications/SphericalHarmonics/ita_sph_real2complex.m +++ b/applications/SphericalHarmonics/ita_sph_real2complex.m @@ -1,5 +1,41 @@ -function SH = ita_sph_real2complex(SH) -% converts a real base or SH vector into its complex base / SH vectors +function SH = ita_sph_real2complex(varargin) +%ITA_SPH_REAL2COMPLEX - Transforms from real to complex valued basis functions +% This function transfroms spherical harmonic coefficients with real +% valued basis functions to their corresponding complex valued +% representation. The sign convention of the real valued basis functions +% can be chosen. The definition in 'Williams - Fourier Acoustics' is +% specified for the complex basis functions (including the Condon-Shotley +% phase). The phase conventions for the real valued spherical harmonics +% can be 'raven' (which equals the ambix phase convention as defined in +% "Nachbar - AMBIX: A Suggested Ambisonics Format (revised by Zotter)") +% or 'zotter' as defined in "Zotter - Analysis and synthesis of +% sound-radiation with spherical arrays" +% +% Syntax: +% SH_cplx = ita_sph_real2complex(SH_real, options) +% +% Options (default): +% 'phase' ('ambix') : Phase definitions as in the AmbiX format +% +% Example: +% SH_cplx = ita_sph_real2complex(SH_real) +% +% See also: +% ita_toolbox_gui, ita_read, ita_write, ita_generate +% +% Reference page in Help browser +% doc ita_sph_real2complex + +% +% 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. +% + +% Rewrite, original author unknown +% +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: 12-Jun-2017 + % % @@ -10,44 +46,20 @@ function SH = ita_sph_real2complex(SH) % SH can be the vector/matrix of coefficients or the maximum order to % calculate a matrix. -if isnumeric(SH) - if numel(SH) == 1 - % return a matrix - SH = ita_sph_real2complex_mat(SH); - else - % assume this as SH coefs - SH = ita_sph_real2complex_coefs(SH); - end -else - error('please check syntax') -end +sArgs = struct('pos1_SH', 'double' ,... + 'phase', 'ambix'); +[SH, sArgs] = ita_parse_arguments(sArgs, varargin); +if ~isnatural(sqrt(size(SH, 1)) - 1) + ita_verbose_info('The dimensions of the input data do not match.', 0) + SH = []; + return end -function SH = ita_sph_real2complex_coefs(SH) +Nmax = floor(sqrt(size(SH,1))-1); -% assume the SH dimension as first dimension -sizeY = size(SH); -% check for max 2 dimension -if numel(sizeY) > 2 - reshape(SH, sizeY(1), []); -end - -% check for SH in 1st dimension -if sizeY(1) == 1 && sizeY(2) > 1 - SH = SH.'; - isTransposed = true; -else - isTransposed = false; -end - -nmax = sqrt(size(SH,1)) - 1; -if nmax ~= round(nmax) - error('ita_sph_complex2real something wrong with the number of coefficients (they do not fill up to a full order)'); -end - -for ind = 1:nmax +for ind = 1:Nmax % define the linear indices for the i'th degree index_m_neg = ita_sph_degreeorder2linear(ind,-1:-1:-ind); % count in reverse order index_m_pos = ita_sph_degreeorder2linear(ind,1:ind); @@ -56,43 +68,19 @@ for ind = 1:nmax rPos = SH(index_m_pos(m),:); rNeg = SH(index_m_neg(m),:); - cPos = ((-1)^m .* rPos + 1i .* rNeg) ./ sqrt(2); - cNeg = (rPos - 1i .* (-1)^m .* rNeg) ./ sqrt(2); - + switch lower(sArgs.phase) + case {'ambix','raven'} + cPos = (rPos + 1i.* rNeg) * (-1)^m / sqrt(2); + cNeg = (rPos - 1i * rNeg) ./ sqrt(2); + case {'zotter'} + cPos = (rPos - 1i.* rNeg) * (-1)^m / sqrt(2); + cNeg = (rPos + 1i * rNeg) ./ sqrt(2); + otherwise + ita_verbose_info('I do not know this phase convention. Aborting...', 0); + SH = []; + return + end SH(ita_sph_degreeorder2linear(ind,m),:) = cPos; SH(ita_sph_degreeorder2linear(ind,-m),:) = cNeg; end -end - - -if isTransposed - SH = SH.'; -end - -% now bring to the old shape -SH = reshape(SH,sizeY); -end -function T = ita_sph_real2complex_mat(nmax) - -if numel(nmax) > 1 - error; -end - -nSH = (nmax+1).^2; - -lin = (1:nSH).'; -multiplicator = 2*mod(lin,2) - 1; - -[deg,order] = ita_sph_linear2degreeorder(lin); -linNeg = ita_sph_degreeorder2linear(deg,-order); - -isPos = order > 0; -isNeg = order < 0; - -T = sparse(eye(nSH)); -T_orig = T; -T_swap = T .* diag(multiplicator); - -T(lin(isPos),:) = (T_swap(lin(isPos),:) + 1i .* T_orig(linNeg(isPos),:)) / sqrt(2); -T(lin(isNeg),:) = (T_orig(linNeg(isNeg),:) - 1i .* T_swap(lin(isNeg),:)) / sqrt(2); -end +end \ No newline at end of file diff --git a/applications/SphericalHarmonics/ita_sph_real_valued_basefunctions.m b/applications/SphericalHarmonics/ita_sph_real_valued_basefunctions.m deleted file mode 100644 index 80375c45a0e407cb75a1f734982581f81d415171..0000000000000000000000000000000000000000 --- a/applications/SphericalHarmonics/ita_sph_real_valued_basefunctions.m +++ /dev/null @@ -1,53 +0,0 @@ -function SH = ita_sph_real_valued_basefunctions(sampling, nmax) - -% -% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved. -% You can find the license for this m-file in the application folder. -% - - -% ita_sph_realvalued_basefunctions(sampling, nmax) -% Returnes the amplitudes of the real valued basefunctions at the sampling's -% points up to order nmax -% -% according to Dissertation Zotter, eq. 26 (page 18) -% Martin Kunkemoeller 16.11.2010 -%% - -ita_verbose_obsolete('Marked as obsolete. Please report to mpo, if you still use this function.'); - -SH = zeros(sampling.nPoints, (nmax+1)^2); - -for n = 0:nmax - idxLeftSide = n^2+n : -1 : n^2+1; - idxRightSide = n^2+n+1: 1 : (n+1)^2; - - - - N = normalize_const(n); - P = legendre(n, cos(sampling.theta)).'; - if size(P,1) ~= sampling.nPoints - P = P.'; - end - - SH(:, idxRightSide) = repmat(N ,[sampling.nPoints 1]) .* P .* cos(sampling.phi * (0:n)); - SH(:, idxLeftSide) = repmat(N(2:end) ,[sampling.nPoints 1]) .* P(:,2:end) .* sin(-sampling.phi * (1:n)); -end -end - -function N = normalize_const(n) -% function for orthogonality (Dis Zotter, eq.31 (page 19)) -N = zeros(1,n+1); -for m = 0:n - N(m+1) = ((-1)^m) * sqrt((2*n+1)*(2-d(m,0))*factorial(n-m)/(4*pi*factorial(n+m))); -end -end - -function out = d(x,m) -% kronecker delta -out = zeros(size(x)); -% if sum(x==m), disp('hello'); end -out(x~=m)=0; -out(x==m)=1; - -end diff --git a/applications/SphericalHarmonics/ita_sph_real_valued_basefunctions_raven.m b/applications/SphericalHarmonics/ita_sph_real_valued_basefunctions_raven.m deleted file mode 100644 index 0f1696ef3da8b59754a944e56956336a9848165f..0000000000000000000000000000000000000000 --- a/applications/SphericalHarmonics/ita_sph_real_valued_basefunctions_raven.m +++ /dev/null @@ -1,55 +0,0 @@ -function SH = ita_sph_real_valued_basefunctions_raven(sampling, nmax) - -% -% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved. -% You can find the license for this m-file in the application folder. -% - - -% ita_sph_realvalued_basefunctions(sampling, nmax) -% Returnes the amplitudes of the real valued basefunctions at the sampling's -% points up to order nmax -% -% according to Dissertation Zotter, eq. 26 (page 18) -% Martin Kunkemoeller 16.11.2010 -%% - -ita_verbose_obsolete('Marked as obsolete. Please report to mpo, if you still use this function.'); - -SH = zeros(sampling.nPoints, (nmax+1)^2); - -for n = 0:nmax - idxLeftSide = n^2+n : -1 : n^2+1; - idxRightSide = n^2+n+1: 1 : (n+1)^2; - - - - N = normalize_const(n); - P = legendre(n, cos(sampling.theta)).'; - if size(P,1) ~= sampling.nPoints - P = P.'; - end - - %% ATTENTION!!!!! - SH(:, idxRightSide) = repmat(N ,[sampling.nPoints 1]) .* P .* cos(sampling.phi * (0:n)); - SH(:, idxLeftSide) = repmat(N(2:end) ,[sampling.nPoints 1]) .* P(:,2:end) .* sin(sampling.phi * (1:n)); %% !!!!!!! P(:,2:end) .* sin(-sampling.phi * (1:n)) in ITA-Toolbox!!!! No minus for RAVEN, just as in Zotters Diss - %% ATTENTION!!!!! -end -end - -function N = normalize_const(n) -% function for orthogonality (Dis Zotter, eq.31 (page 19)) -N = zeros(1,n+1); -for m = 0:n - N(m+1) = ((-1)^m) * sqrt((2*n+1)*(2-d(m,0))*factorial(n-m)/(4*pi*factorial(n+m))); -end -end - -function out = d(x,m) -% kronecker delta -out = zeros(size(x)); -% if sum(x==m), disp('hello'); end -out(x~=m)=0; -out(x==m)=1; - -end diff --git a/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_dodecahedron.m b/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_dodecahedron.m index 0ef3dbc67f5d50f2c913a47ad323a51c26cc4c0b..b42390006ed48c37112da006a63998bb81d3b1e6 100644 --- a/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_dodecahedron.m +++ b/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_dodecahedron.m @@ -44,6 +44,9 @@ phi3 = 4*pi/3; theta = [repmat(theta1,3,1); repmat(theta2,3,1); repmat(theta3,3,1); repmat(theta4,3,1)].'; phi = repmat([[phi1; phi2; phi3]; pi/3+[phi1; phi2; phi3]],2,1).'; -r = 0.15 .* ones(size(theta)); +r = ones(size(theta)); -s = itaSamplingSph([r(:) theta(:) phi(:)],'sph'); \ No newline at end of file +s = itaSamplingSph([r(:) theta(:) phi(:)],'sph'); +s.nmax = floor(sqrt(s.nPoints)-1); + +end \ No newline at end of file diff --git a/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_equalarea.m b/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_equalarea.m index d637eaded541402bf952241ec51dcc63ce6d32a2..b038c94f5f705e7084c1b279e8c9107263856c9c 100644 --- a/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_equalarea.m +++ b/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_equalarea.m @@ -52,8 +52,12 @@ while true coordsCart = eq_point_set(2,sArgs.nPoints).'; sampling = itaSamplingSph(coordsCart,'cart'); Y = ita_sph_base(sampling,Nmax); - condNum = cond(Y); - if condNum < sArgs.condSHT + if sArgs.condSHT ~= inf + condNum = cond(Y); + if condNum < sArgs.condSHT + break; + end + else break; end sArgs.nPoints = sArgs.nPoints+1; diff --git a/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_spiral_points.m b/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_spiral_points.m index a6ce5d3f22f3ae845fd3f24c4bfad00bd843c5b2..7229b9f32dea181c8b60d8e4446809487c1f3754 100644 --- a/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_spiral_points.m +++ b/applications/SphericalHarmonics/ita_sph_sampling/ita_sph_sampling_spiral_points.m @@ -53,8 +53,12 @@ while true coordsCart = calculate_spiral_points(sArgs.nPoints); sampling = itaSamplingSph(coordsCart,'sph'); Y = ita_sph_base(sampling,Nmax); - condNum = cond(Y); - if condNum < sArgs.condSHT + if sArgs.condSHT ~= inf + condNum = cond(Y); + if condNum < sArgs.condSHT + break; + end + else break; end sArgs.nPoints = sArgs.nPoints+1; diff --git a/applications/VirtualAcoustics/Raven/itaRavenProject.m b/applications/VirtualAcoustics/Raven/itaRavenProject.m index 70134cebee3ac12941cfe989d3733bac6233f4a5..22c791f43c4d19cb86132d7b3ec55b9cce85bac3 100644 --- a/applications/VirtualAcoustics/Raven/itaRavenProject.m +++ b/applications/VirtualAcoustics/Raven/itaRavenProject.m @@ -274,17 +274,26 @@ classdef itaRavenProject < handle %------------------------------------------------------------------ function setRavenExe(obj, newRavenExe) - obj.ravenExe = newRavenExe; - - if (exist(obj.ravenIniFile,'file')) - obj.raven_ini.SetValues('Global', {'PathRavenExe'}, {obj.ravenExe}); - obj.raven_ini.WriteFile(obj.ravenIniFile); + + + 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 - obj.raven_ini = IniConfig(); - obj.raven_ini.AddSections({'Global'}); - obj.raven_ini.AddKeys('Global', {'PathRavenExe'}, {obj.ravenExe}); + error('[itaRaven]: Error: Path to new Raven binary not found!'); + end - + end %------------------------------------------------------------------ @@ -444,29 +453,30 @@ classdef itaRavenProject < handle % 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]); + + if (~exist(obj.ravenExe,'file')) + error('[itaRaven]: Error: Cannot find Raven binary file!'); + end + 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('Getting results...'); + disp('[R] Simulation seems to be finished. Getting results...'); obj.gatherResults(); - disp('Done.'); + disp('[R] Done.'); obj.simulationDone = true; @@ -505,8 +515,12 @@ classdef itaRavenProject < handle delete(obj.ravenLogFile); end % system([obj.ravenExe ' "' obj.ravenProjectFile '" >> ' obj.ravenLogFile]); - dos([obj.ravenExe ' "' obj.ravenProjectFile '"'], '-echo'); + prevPath = pwd; + cd(fileparts(obj.ravenExe)); + dos(['"' obj.ravenExe '"' ' "' obj.ravenProjectFile '"'],'-echo'); disp('Done.'); + cd(prevPath); + % restore the initial project name obj.setProjectName(savedProjectName); @@ -749,12 +763,12 @@ classdef itaRavenProject < handle 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)),'%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; - currentMaterial.pf; + ita_plot_freq(currentMaterial,'LineWidth',2); % change format of plot title(''); @@ -764,7 +778,7 @@ classdef itaRavenProject < handle set(gca,'YLim',[0 1]); leg = findobj(gcf,'Tag','legend'); set(leg,'Location','NorthWest'); - set(leg,'FontSize',9); + set(leg,'FontSize',12); % remove [1] in legend entry for iMat=1:numberMaterials @@ -5418,7 +5432,7 @@ classdef itaRavenProject < handle 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 diff --git a/applications/VirtualAcoustics/Raven/ita_raven_demo.m b/applications/VirtualAcoustics/Raven/ita_raven_demo.m index 5e2a04e7600efed40161fbb87b64502d62eaac82..965aff3fa5a682b59be371096eb30c220f23b1c5 100644 --- a/applications/VirtualAcoustics/Raven/ita_raven_demo.m +++ b/applications/VirtualAcoustics/Raven/ita_raven_demo.m @@ -13,8 +13,9 @@ ravenProjectPath = '..\RavenInput\Classroom\Classroom.rpf'; if (~exist(ravenProjectPath,'file')) [filename, pathname] = uigetfile('Classroom.rpf', 'Please select raven project file!'); ravenProjectPath = [pathname filename]; - ravenBasePath = ravenProjectPath(1:end-34); end +ravenBasePath = ravenProjectPath(1:end-34); + rpf = itaRavenProject(ravenProjectPath); %% Simulationsparameter einstellen % Image sources up to second order diff --git a/applications/VirtualAcoustics/Raven/readRavenMaterial.m b/applications/VirtualAcoustics/Raven/readRavenMaterial.m index e5323d2f8b2b8e8c48855c2575a2c6b0f9d97ac2..14cf2d49ce9548c36277649718d49f004572fdc9 100644 --- a/applications/VirtualAcoustics/Raven/readRavenMaterial.m +++ b/applications/VirtualAcoustics/Raven/readRavenMaterial.m @@ -1,24 +1,24 @@ function [absorp, scatter] = readRavenMaterial(materialName, pathToMaterials) - - -% -% 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. -% + % + % 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. + % if nargin < 2 pathToMaterials = '..\RavenDatabase\MaterialDatabase'; - end + end - if isequal(materialName(end-4:end), '.mat') - materialName = materialName(1:end-4); + if (length(materialName) > 4) + if isequal(materialName(end-4:end), '.mat') + materialName = materialName(1:end-4); + end end - ini = IniConfig(); + ini = IniConfig(); ini.ReadFile(fullfile(pathToMaterials, [materialName '.mat'])); absorp = ini.GetValues('Material', 'absorp'); - + if nargout > 1 scatter = ini.GetValues('Material', 'scatter'); end diff --git a/applications/VirtualAcoustics/VA/itaVA.m b/applications/VirtualAcoustics/VA/itaVA.m index 738420da5cf78862ee9e3155aa707c5e85feb721..c58575dfeb96fad2d8253a59e9a5699bd5f09b4b 100644 --- a/applications/VirtualAcoustics/VA/itaVA.m +++ b/applications/VirtualAcoustics/VA/itaVA.m @@ -235,8 +235,15 @@ classdef itaVA < handle connected = VAMatlab( 'IsTrackerConnected', this.handle ); end + function disconnectTracker( this ) + % Disconnects from the NatNet tracking server + VAMatlab( 'DisconnectTracker', this.handle ) + end + + % -- Tracked listener -- % + function setTrackedListener( this, listener_id ) - % Connects a VA listener with the tracked rigid body + % Connects a VA listener with the tracked listener rigid body % % Parameters: % @@ -245,38 +252,101 @@ classdef itaVA < handle VAMatlab( 'SetTrackedListener', this.handle, listener_id ); end - function setTrackedSource( this, source_id ) - % Connects a VA source with the tracked rigid body + function setTrackedListenerRigidBodyIndex( this, index ) + % Sets the index of the rigid body to be tracked for listener (default is 1) + VAMatlab( 'SetTrackedListenerRigidBodyIndex', this.handle, index ) + end + + function setTrackedListenerRigidBodyTranslation( this, translation ) + % Sets the pivot point translation for the tracked listener rigid body + % + % Parameters: + % + % translation [double-3x1] Translation in local coordinate system of rigid body [m] + % + VAMatlab( 'SetTrackedListenerRigidBodyTranslation', this.handle, translation ) + end + + function setTrackedListenerRigidBodyRotation( this, rotation ) + % Sets the rotation of orientation for the tracked listener rigid body + % + % Given rotation has to be a Matlab quaternion type (order: w(real), i, j, k) + % + % Parameters: + % + % rotation [quaternion] Rotation of rigid body + % + VAMatlab( 'SetTrackedListenerRigidBodyRotation', this.handle, rotation ) + end + + % -- Tracked real-world listener -- % + + function setTrackedRealWorldListener( this, listener_id ) + % Connects a VA real-world listener with the tracked real-world rigid body % % Parameters: % - % source_id [integer-1x1] VA listener id + % listener_id [integer-1x1] VA real-world listener id % - VAMatlab( 'SetTrackedSource', this.handle, source_id ); + VAMatlab( 'SetTrackedRealWorldListener', this.handle, listener_id ); + end + + function setTrackedRealWorldListenerRigidBodyIndex( this, index ) + % Sets the index of the rigid body to be tracked for real-world listener (default is 1) + VAMatlab( 'SetTrackedRealWorldListenerRigidBodyIndex', this.handle, index ) end - function disconnectTracker( this ) - % Disconnects from the NatNet tracking server - VAMatlab( 'DisconnectTracker', this.handle ) + function setTrackedRealWorldListenerRigidBodyTranslation( this, translation ) + % Sets the pivot point translation for the tracked real-world listener rigid body + % + % Parameters: + % + % translation [double-3x1] Translation in local coordinate system of rigid body [m] + % + VAMatlab( 'SetTrackedRealWorldListenerRigidBodyTranslation', this.handle, translation ) end - function setRigidBodyIndex( this, index ) - % Sets the index of the rigid body to be tracked (default is 1) - VAMatlab( 'SetRigidBodyIndex', this.handle, index ) + function setTrackedRealWorldListenerRigidBodyRotation( this, rotation ) + % Sets the rotation of orientation for the tracked real-world listener rigid body + % + % Given rotation has to be a Matlab quaternion type (order: w(real), i, j, k) + % + % Parameters: + % + % rotation [quaternion] Rotation of rigid body + % + VAMatlab( 'SetTrackedRealWorldListenerRigidBodyRotation', this.handle, rotation ) + end + + % -- Tracked source -- % + + function setTrackedSource( this, source_id ) + % Connects a VA source with the tracked source rigid body + % + % Parameters: + % + % source_id [integer-1x1] VA source id + % + VAMatlab( 'SetTrackedSource', this.handle, source_id ); + end + + function setTrackedSourceRigidBodyIndex( this, index ) + % Sets the index of the rigid body to be tracked for source (default is 1) + VAMatlab( 'SetTrackedSourceRigidBodyIndex', this.handle, index ) end - function setRigidBodyTranslation( this, translation ) - % Sets the pivot point translation for the tracked rigid body + function setTrackedSourceRigidBodyTranslation( this, translation ) + % Sets the pivot point translation for the tracked source rigid body % % Parameters: % % translation [double-3x1] Translation in local coordinate system of rigid body [m] % - VAMatlab( 'SetRigidBodyTranslation', this.handle, translation ) + VAMatlab( 'SetTrackedSourceRigidBodyTranslation', this.handle, translation ) end - function setRigidBodyRotation( this, rotation ) - % Sets the rotation of orientation for the tracked rigid body + function setTrackedSourceRigidBodyRotation( this, rotation ) + % Sets the rotation of orientation for the tracked source rigid body % % Given rotation has to be a Matlab quaternion type (order: w(real), i, j, k) % @@ -284,8 +354,9 @@ classdef itaVA < handle % % rotation [quaternion] Rotation of rigid body % - VAMatlab( 'SetRigidBodyRotation', this.handle, rotation ) + VAMatlab( 'SetTrackedSourceRigidBodyRotation', this.handle, rotation ) end + %% --= Functions =-- @@ -481,6 +552,24 @@ classdef itaVA < handle [id] = VAMatlab('createSoundSourceExplicitRenderer', this.handle, name,renderer); end + function [signalSourceID] = createTextToSpeechSignalSource(this, name) + % Creates a text to speech signal + % + % Parameters: + % + % name [string] Displayed name (optional, default: '') + % + % Return values: + % + % signalSourceID [string] Signal source ID + % + + if this.handle==0, error('Not connected.'); end; + + if ~exist('name','var'), name = ''; end + [signalSourceID] = VAMatlab('createTextToSpeechSignalSource', this.handle, name); + end + function [] = deleteListener(this, listenerID) % Deletes a listener from the scene % @@ -1209,6 +1298,24 @@ classdef itaVA < handle [info] = VAMatlab('getSignalSourceInfos', this.handle); end + function [params] = getSignalSourceParameters(this, ID,args) + % Returns the current signal source parameters + % + % Parameters: + % + % ID [string] Signal source identifier + % args [mstruct] Requested parameters + % + % Return values: + % + % params [mstruct] Parameters + % + + if this.handle==0, error('Not connected.'); end; + + [params] = VAMatlab('getSignalSourceParameters', this.handle, ID,args); + end + function [info] = getSoundInfo(this, soundID) % Returns information on a loaded sound % @@ -2279,6 +2386,24 @@ classdef itaVA < handle VAMatlab('setReproductionModuleMuted', this.handle, sModuleID,bMuted); end + function [] = setSignalSourceParameters(this, ID,params) + % Sets signal source parameters + % + % Parameters: + % + % ID [string] Signal source identifier + % params [mstruct] Parameters + % + % Return values: + % + % None + % + + if this.handle==0, error('Not connected.'); end; + + VAMatlab('setSignalSourceParameters', this.handle, ID,params); + end + function [] = setSoundSourceAuralizationMode(this, soundSourceID,auralizationMode) % Returns the auralization mode of a sound source % diff --git a/applications/VirtualAcoustics/VA/itaVA_NCTC_gui.fig b/applications/VirtualAcoustics/VA/itaVA_NCTC_gui.fig new file mode 100644 index 0000000000000000000000000000000000000000..86f70e9df9a998fee333770de2077cfe98670b21 Binary files /dev/null and b/applications/VirtualAcoustics/VA/itaVA_NCTC_gui.fig differ diff --git a/applications/VirtualAcoustics/VA/itaVA_NCTC_gui.m b/applications/VirtualAcoustics/VA/itaVA_NCTC_gui.m new file mode 100644 index 0000000000000000000000000000000000000000..f59685bb3e82642489ef1111d24d1e40a9d23c66 --- /dev/null +++ b/applications/VirtualAcoustics/VA/itaVA_NCTC_gui.m @@ -0,0 +1,176 @@ +function varargout = itaVA_NCTC_gui(varargin) +% ITAVA_NCTC_GUI MATLAB code for itaVA_NCTC_gui.fig +% ITAVA_NCTC_GUI, by itself, creates a new ITAVA_NCTC_GUI or raises the existing +% singleton*. +% +% H = ITAVA_NCTC_GUI returns the handle to a new ITAVA_NCTC_GUI or the handle to +% the existing singleton*. +% +% ITAVA_NCTC_GUI('CALLBACK',hObject,eventData,handles,...) calls the local +% function named CALLBACK in ITAVA_NCTC_GUI.M with the given input arguments. +% +% ITAVA_NCTC_GUI('Property','Value',...) creates a new ITAVA_NCTC_GUI or raises the +% existing singleton*. Starting from the left, property value pairs are +% applied to the GUI before itaVA_NCTC_gui_OpeningFcn gets called. An +% unrecognized property name or invalid value makes property application +% stop. All inputs are passed to itaVA_NCTC_gui_OpeningFcn via varargin. +% +% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one +% instance to run (singleton)". +% +% See also: GUIDE, GUIDATA, GUIHANDLES + +% Edit the above text to modify the response to help itaVA_NCTC_gui + +% Last Modified by GUIDE v2.5 20-Jun-2017 15:50:42 + +% Begin initialization code - DO NOT EDIT +gui_Singleton = 1; +gui_State = struct('gui_Name', mfilename, ... + 'gui_Singleton', gui_Singleton, ... + 'gui_OpeningFcn', @itaVA_NCTC_gui_OpeningFcn, ... + 'gui_OutputFcn', @itaVA_NCTC_gui_OutputFcn, ... + 'gui_LayoutFcn', [] , ... + 'gui_Callback', []); +if nargin && ischar(varargin{1}) + gui_State.gui_Callback = str2func(varargin{1}); +end + +if nargout + [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); +else + gui_mainfcn(gui_State, varargin{:}); +end +% End initialization code - DO NOT EDIT + + +% --- Executes just before itaVA_NCTC_gui is made visible. +function itaVA_NCTC_gui_OpeningFcn(hObject, eventdata, handles, varargin) +% This function has no output args, see OutputFcn. +% hObject handle to figure +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +% varargin command line arguments to itaVA_NCTC_gui (see VARARGIN) + +% Choose default command line output for itaVA_NCTC_gui +handles.output = hObject; + +handles.va = itaVA; +handles.repros = []; +handles.current_repro_params = struct(); + +% Update handles structure +guidata(hObject, handles); + +% UIWAIT makes itaVA_NCTC_gui wait for user response (see UIRESUME) +% uiwait(handles.figure1); + + +% --- Outputs from this function are returned to the command line. +function varargout = itaVA_NCTC_gui_OutputFcn(hObject, eventdata, handles) +% varargout cell array for returning output args (see VARARGOUT); +% hObject handle to figure +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Get default command line output from handles structure +varargout{1} = handles.output; + + +% --- Executes on button press in pushbutton_connect. +function pushbutton_connect_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton_connect (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +handles.va.connect( handles.edit_va_server_ip.String ) + +if ~handles.va.isConnected + error( 'Not connected to VA' ) +end + +repros = handles.va.getReproductionModules(); +for n=1:numel( repros ) + repro = repros( n ); + if strcmp( repro.class, 'NCTC' ) + handles.repros = [ handles.repros repro ]; + end +end + +if numel( handles.repros ) == 0 + error( 'No NCTC reproduction module found, please activate in VA configuration' ) +end + +in_args.info = true; +handles.current_repro_params = handles.va.callModule( [ 'NCTC:' handles.repros( 1 ).id ], in_args ); + +handles.slider_cross_talk_cancellation_factor; + +% Update handles structure +guidata(hObject, handles); + + +function edit_va_server_ip_Callback(hObject, eventdata, handles) +% hObject handle to edit_va_server_ip (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit_va_server_ip as text +% str2double(get(hObject,'String')) returns contents of edit_va_server_ip as a double + + +% --- Executes during object creation, after setting all properties. +function edit_va_server_ip_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit_va_server_ip (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +% --- Executes on selection change in listbox_repros. +function listbox_repros_Callback(hObject, eventdata, handles) +% hObject handle to listbox_repros (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: contents = cellstr(get(hObject,'String')) returns listbox_repros contents as cell array +% contents{get(hObject,'Value')} returns selected item from listbox_repros + + +% --- Executes during object creation, after setting all properties. +function listbox_repros_CreateFcn(hObject, eventdata, handles) +% hObject handle to listbox_repros (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: listbox controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +% --- Executes on slider movement. +function slider_cross_talk_cancellation_factor_Callback(hObject, eventdata, handles) +% hObject handle to slider_cross_talk_cancellation_factor (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'Value') returns position of slider +% get(hObject,'Min') and get(hObject,'Max') to determine range of slider + + +% --- Executes during object creation, after setting all properties. +function slider_cross_talk_cancellation_factor_CreateFcn(hObject, eventdata, handles) +% hObject handle to slider_cross_talk_cancellation_factor (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: slider controls usually have a light gray background. +if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor',[.9 .9 .9]); +end diff --git a/applications/VirtualAcoustics/VA/itaVA_example_text_to_speech.m b/applications/VirtualAcoustics/VA/itaVA_example_text_to_speech.m new file mode 100644 index 0000000000000000000000000000000000000000..70f6dfee236e85a709e61cbf4c8edb8941227e5f --- /dev/null +++ b/applications/VirtualAcoustics/VA/itaVA_example_text_to_speech.m @@ -0,0 +1,79 @@ +%% itaVA simple example code for a text-to-speech sound source + +% Quick setup +itaVAq +va.reset(); +va.setOutputGain( .25 ); + +% Create a text-to-speech audio signal source +X = va.createTextToSpeechSignalSource( 'tts_demo' ); + +% Create a virtual sound source and set a position +S = va.createSoundSource( 'itaVA_Source' ); +va.setSoundSourcePosition( S, [ -2 1.7 -1 ] ) + +% Connect the signal source to the virtual sound source +va.setSoundSourceSignalSource( S, X ) + +% Create a listener with an HRTF and position him +H = va.loadHRIRDataset( '$(DefaultHRIR)' ); +L = va.createListener( 'itaVA_Listener', 'default', H ); +va.setListenerPosition( L, [ 0 1.7 0 ] ) + +% Set the listener as the active one +va.setActiveListener( L ) + +%% Control TTS +tts_in = struct(); +tts_in.list_voices = true; +tts_out = va.getSignalSourceParameters( X, tts_in ) +assert( tts_out.number > 0 ) + +tts_in = struct(); +tts_in.voice = 'Heather'; +tts_in.id = 'id_vr'; +tts_in.prepare_text = 'virtual acoustics is a real-time auralization framework for scientific research in Virtual Reality created by the institute of technical acoustics, RWTH aachen university'; +tts_in.direct_playback = true; +va.setSignalSourceParameters( X, tts_in ) +pause( 5 ) + +tts_in = struct(); +tts_in.id = 'id_ok'; +tts_in.prepare_text = 'OK'; +va.setSignalSourceParameters( X, tts_in ) + +tts_in = struct(); +tts_in.id = 'id_alright'; +tts_in.prepare_text = 'Alright'; +va.setSignalSourceParameters( X, tts_in ) + +tts_in_ok.play_speech = 'id_ok'; +tts_in_alright.play_speech = 'id_alright'; + +va.setSignalSourceParameters( X, tts_in_alright ) +pause( 1 ) +va.setSignalSourceParameters( X, tts_in_ok ) +pause( 1.5 ) +va.setSignalSourceParameters( X, tts_in_ok ) +pause( 1 ) +va.setSignalSourceParameters( X, tts_in_alright ) + +tts_in = struct(); +tts_in.id = 'id_excidea'; +tts_in.prepare_text = 'That is an excellent idea!'; +tts_in.direct_playback = true; +va.setSignalSourceParameters( X, tts_in ) + +while( true ) + txt = inputdlg; + if isempty( txt ) + break; + end + tts_in = struct(); + tts_in.direct_playback = true; + tts_in.prepare_text = txt{1}; + if isempty( tts_in.prepare_text ) + break; + end + va.setSignalSourceParameters( X, tts_in ) +end diff --git a/applications/VirtualAcoustics/VA/itaVA_example_tracked_listener.m b/applications/VirtualAcoustics/VA/itaVA_example_tracked_listener.m index 4b5e25b64bb6ff73bb8d50e66651d2637fdac9b4..0d578efaa7c8f13a9aec87e08c09eb7cc3bccebc 100644 --- a/applications/VirtualAcoustics/VA/itaVA_example_tracked_listener.m +++ b/applications/VirtualAcoustics/VA/itaVA_example_tracked_listener.m @@ -8,18 +8,9 @@ va = itaVA( 'localhost' ) L = va.createListener( 'itaVA_Tracked_Listener' ); % OptiTrack tracker conneection and listener updates -va.setTrackedListener( L ) +va.setTrackedListener( L ) % For virtual scene / rendering +va.setTrackedRealWorldListener( L ) % For CTC reproductions va.connectTracker - -% apply pivot point offset to a rigid body -% Hint: the method .calibrate of itaOptirack() calculates the individual -% offset between a head-mounted rigid body and the center of the interaural -% axis of a listener -va.setRigidBodyIndex( 1 ) % set index of rigid body that should be manipulated (cf. Motive) -va.setRigidBodyTranslation( [0 -0.08 0] ) % translation in local coordinate system of rigid body [m] - % move rigid body by 8 cm in - % negative y direction - pause( 12 ) % Observe how you can move the virtual listener in VAGUI va.disconnectTracker diff --git a/applications/VirtualAcoustics/VA/itaVA_experimental_gui.m b/applications/VirtualAcoustics/VA/itaVA_experimental_gui.m index 01a9db67de87d883975fd5a4a6849d338aaaf614..569d18b7b905f67cde7f52ae30a17a10cbea56de 100644 --- a/applications/VirtualAcoustics/VA/itaVA_experimental_gui.m +++ b/applications/VirtualAcoustics/VA/itaVA_experimental_gui.m @@ -59,6 +59,9 @@ handles.module_id = 'PrototypeGenericPath:MyGenericRenderer'; handles.va_source_id = -1; handles.va_signal_id = ''; handles.va_listener_id = -1; +handles.module_channels = -1; +handles.module_filter_length_samples = -1; +handles.listbox_sourcesignals_last_index = -1; refresh_workspace_vars( hObject, handles ); refresh_sourcesignals( hObject, handles ); @@ -116,10 +119,12 @@ disp( [ 'Using channel prototype generic path renderer with identifier: ' gpg_re handles.module_id = [ gpg_renderer.class ':' gpg_renderer.id ]; in_args.info = true; out_args = handles.va.callModule( handles.module_id, in_args ); -disp( [ 'Your experimental renderer "' handles.module_id '" has ' num2str( out_args.numchannels ) ' channels and an FIR filter length of ' num2str( out_args.irfilterlengthsamples ) ' samples' ] ) +handles.module_channels = out_args.numchannels; +handles.module_filter_length_samples = out_args.irfilterlengthsamples; +disp( [ 'Your experimental renderer "' handles.module_id '" has ' num2str( handles.module_channels ) ' channels and an FIR filter length of ' num2str( out_args.irfilterlengthsamples ) ' samples' ] ) -handles.edit_va_channels.String = out_args.numchannels; -handles.edit_va_fir_taps.String = out_args.irfilterlengthsamples; +handles.edit_va_channels.String = handles.module_channels; +handles.edit_va_fir_taps.String = handles.module_filter_length_samples; handles.edit_va_fs.String = '44.100'; % @todo get from VA audio streaming settings % Useful FIRs @@ -141,6 +146,8 @@ handles.va_signal_id = ''; handles.va.setOutputMuted( handles.checkbox_global_muted.Value ); handles.va.setOutputGain( handles.slider_volume.Value ); +refresh_sourcesignals( hObject, handles ); +refresh_workspace_vars( hObject, handles ); % Update handles (store values) guidata( hObject, handles ); @@ -219,31 +226,35 @@ filename_selected = filename_list{ index_selected }; if handles.va.isConnected - va_signal_id_old = handles.va_signal_id; + last_index = handles.listbox_sourcesignals_last_index; - if ~isempty( va_signal_id_old ) && strcmp( handles.va.getSoundSourceSignalSource( handles.va_source_id ), va_signal_id_old ) - playstate = handles.va.getAudiofileSignalSourcePlaybackState( va_signal_id_old ); + if last_index == index_selected + % play/pause current + va_signal_id = handles.va.getSoundSourceSignalSource( handles.va_source_id ); + assert( ~isempty( va_signal_id ) ) + playstate = handles.va.getAudiofileSignalSourcePlaybackState( va_signal_id ); if strcmpi( playstate, 'playing' ) - handles.va.setAudiofileSignalSourcePlaybackAction( handles.va_signal_id, 'pause' ); + handles.va.setAudiofileSignalSourcePlaybackAction( va_signal_id, 'pause' ); else - handles.va.setAudiofileSignalSourcePlaybackAction( handles.va_signal_id, 'play' ); + handles.va.setAudiofileSignalSourcePlaybackAction( va_signal_id, 'play' ); end else - if ~isempty( va_signal_id_old ) + % new selected, remove old? + if last_index ~= -1 + va_signal_id = handles.va.getSoundSourceSignalSource( handles.va_source_id ); handles.va.setSoundSourceSignalSource( handles.va_source_id, '' ); - handles.va.deleteSignalSource( va_signal_id_old ); + handles.va.deleteSignalSource( va_signal_id ); end + % create new handles.va_signal_id = handles.va.createAudiofileSignalSource( filepath_selected, filename_selected ); handles.va.setSoundSourceSignalSource( handles.va_source_id, handles.va_signal_id ); handles.va.setAudiofileSignalSourcePlaybackAction( handles.va_signal_id, 'play' ); is_looping = handles.checkbox_loop.Value; handles.va.setAudiofileSignalSourceIsLooping( handles.va_signal_id, is_looping ); + handles.listbox_sourcesignals_last_index = index_selected; - if ~isempty( va_signal_id_old ) - handles.va.deleteSignalSource( va_signal_id_old ); - end end end @@ -271,21 +282,35 @@ function pushbutton_start_va_Callback(hObject, eventdata, handles) % handles structure with handles and user data (see GUIDATA) itaVA_experimental_start_server + function refresh_workspace_vars( hObject, handles ) +% Updates workspace variables in listbox base_ws_vars = evalin( 'base', 'whos' ); stringlist = ''; for i=1:numel( base_ws_vars ) if( strcmp( base_ws_vars( i ).class, 'itaAudio' ) ) - stringlist = [ stringlist; { base_ws_vars( i ).name } ]; + audio_var = evalin( 'base', base_ws_vars( i ).name ); + if handles.module_channels == audio_var.nChannels + stringlist = [ stringlist; { base_ws_vars( i ).name } ]; + end end end -handles.listbox_filters.String = stringlist; +if ~isempty( stringlist ) + handles.listbox_filters.String = stringlist; +else + if handles.va.isConnected + warning( [ 'No itaAudio objects with matching ' handles.module_channels ' channels found in current workspace' ] ) + else + %warning( [ 'No itaAudio objects found in current workspace' ] ) + end +end function refresh_sourcesignals( hObject, handles ) filelist = dir( pwd ); +handles.listbox_sourcesignals_last_index = -1; stringlist = ''; fullfile_stringlist = ''; @@ -298,8 +323,18 @@ for i=1:numel( filelist ) end end -handles.listbox_sourcesignals.String = stringlist; -handles.listbox_sourcesignals.UserData = fullfile_stringlist; +if ~isempty( stringlist ) + handles.listbox_sourcesignals.String = stringlist; +else + %warning( [ 'No WAV files found in current workfolder (' pwd ')' ] ) +end +if ~isempty( fullfile_stringlist ) + handles.listbox_sourcesignals.UserData = fullfile_stringlist; +end + +% Update handles structure +guidata( hObject, handles ); + % --- Executes on button press in pushbutton_refresh_workspace_vars. diff --git a/applications/VirtualAcoustics/openDAFF/OpenDAFFv1.5/daffv15_write.m b/applications/VirtualAcoustics/openDAFF/OpenDAFFv1.5/daffv15_write.m index 09278c2031402cea89af37fda4778ff27ad632b0..a154e59a58a74f90ad9be404bc61e36adc1ef9f3 100644 --- a/applications/VirtualAcoustics/openDAFF/OpenDAFFv1.5/daffv15_write.m +++ b/applications/VirtualAcoustics/openDAFF/OpenDAFFv1.5/daffv15_write.m @@ -902,7 +902,8 @@ for i=1:props.dataset.numrecords % Back to Matlab indices (+1) i1 = recordDesc{i}.effOffset(c) + 1; - i2 = i1 + recordDesc{i}.effLength(c) - 1; +% i2 = i1 + recordDesc{i}.effLength(c) - 1; + i2 = recordDesc{i}.effLength(c) - 1; % Write down the filter coefficients switch props.quantization diff --git a/external_packages/sofa/itaFunctions/ita_sofa_getCoordinates.m b/external_packages/sofa/itaFunctions/ita_sofa_getCoordinates.m index cfac1179843b097319967b5abb23e21c872897ce..0483ae6a825f9de55052ecb59d900c161f87e404 100644 --- a/external_packages/sofa/itaFunctions/ita_sofa_getCoordinates.m +++ b/external_packages/sofa/itaFunctions/ita_sofa_getCoordinates.m @@ -50,7 +50,7 @@ switch dataType case 'spherical' % if strcmpi(dataUnit,'degree') coordinates.phi_deg = data(:,1); - coordinates.elevation = data(:,2); + coordinates.theta_deg = data(:,2)+90; coordinates.r = data(:,3); % else % warning('ITA_READ_SOFA: unit is not checked in spherical type'); diff --git a/ita_toolbox_license.m b/ita_toolbox_license.m index eeec95f3502c6e87510abac022b48770ffe50e69..a2aecb97514b1e24e8fa501b4ad9ecdb2c56446e 100644 --- a/ita_toolbox_license.m +++ b/ita_toolbox_license.m @@ -17,7 +17,7 @@ fprintf(2,'********************************************************************* %% ask if usejava('desktop') %Only if jvm available (non_cluster) commandwindow(); - choice = questdlg('Do you agree to the terms of the license agreement? The complete licnese ("license.txt") can be found in the root directory of the ITA-Toolbox.','License Agreement - ITA-Toolbox:', ... + choice = questdlg('Do you agree to the terms of the license agreement? The complete license ("license.txt") can be found in the root directory of the ITA-Toolbox.','License Agreement - ITA-Toolbox:', ... 'Yes','No','No'); switch lower(choice) case 'yes' diff --git a/kernel/ClassStuff/@itaCoordinates/plot_map_projection.m b/kernel/ClassStuff/@itaCoordinates/plot_map_projection.m new file mode 100644 index 0000000000000000000000000000000000000000..44342506345d4e18e900927c8d14f7e514ef7b1a --- /dev/null +++ b/kernel/ClassStuff/@itaCoordinates/plot_map_projection.m @@ -0,0 +1,135 @@ +function varargout = plot_map_projection(this, varargin) +%PLOT_MAP_PROJECTION - Plot a map projection of data on a sphere +% This function plots a map projection for a dataset on a full or +% partial spherical surface +% +% Syntax: +% plot_map_projection(this, data, options) +% +% Options (default): +% 'proj' ('wagner4') : Specifies the projection type (See mapping toolbox documentation for supported projections) +% 'phi0' (0) : Longitute angle phi0 = 0 -> x-axis at origin of the map +% 'shading' ('interp') : Specifies the shading algorithm +% 'db' (false) : db plot +% 'limits' ([]) : Colorbar limits +% +% Example: +% plot_map_projection(this, data, 'db', 'limits', [-50, 0]) +% +% See also: +% itaCoordinates, ita_sph_sampling, surf, scatter +% +% Reference page in Help browser +% doc plot_map_projection + +% +% 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. +% + + +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: 06-Jun-2017 + +sArgs = struct('pos1_data','double',... + 'phi0',0,... + 'proj','wagner4',... + 'db',false,... + 'shading','interp',... + 'limits',[],... + 'fgh',[],... + 'colormap',[]); +[data,sArgs] = ita_parse_arguments(sArgs,varargin); + +% convert to longitude, latitude coordinates required by the mapping toolbox +[azi,ele,~] = cart2sph(this.x,this.y,this.z); +lat = ele*180/pi; +% if phi0 = 0, x = 1 is the origin of the map at lat = 0, lon = 0 +lon = (azi-sArgs.phi0)*180/pi; +% longitude, latitude limits for the plot +latLim = [min(lat),max(lat)]; +lonLim = [min(lon),max(lon)]; + +if isempty(sArgs.fgh) + sArgs.fgh = figure(); +end +% project sampling coordinates +mstruct = defaultm(sArgs.proj); +mstruct.origin = [0,0]; +mstruct = defaultm(mstruct); +mstruct.maplonlimit = lonLim; +mstruct.maplatlimit = latLim; +mstruct.flinewidth = 1; +mstruct.flatlimit = latLim; +mstruct.flonlimit = lonLim; +[x,y] = mfwdtran(mstruct,lat,lon); + +LatTicks = 0:30:latLim(2); +LatTicks = [-fliplr(30:30:abs(latLim(1))),LatTicks(1:end)]; +LonTicks = 0:50:lonLim(2); +LonTicks = [-fliplr(50:50:abs(lonLim(1))),LonTicks(1:end)]; + +[latTicks,lonTicks] = meshgrid(latLim(1),LonTicks); +[XTicks,~] = mfwdtran(mstruct,latTicks,lonTicks); +[latTicks,lonTicks] = meshgrid(LatTicks,lonLim(2)); +[~,YTicks] = mfwdtran(mstruct,latTicks,lonTicks); + + +% calculate delaunay triangulation for the meshgrid +tri = delaunay(x,y); +if sArgs.db + data = 20*log10(abs(data)); +else + data = abs(data); +end +trisurf(tri,x,y,data); + +cb = colorbar; +cb.Label.String = 'Magnitude'; +if ~isempty(sArgs.limits) + caxis(sArgs.limits); +end + + +% use custom colormap if given +if ~isempty(sArgs.colormap) + colormap(sArgs.fgh,sArgs.colormap); +end +% set shading, will be set to interpolate by default +shading(sArgs.shading); + +% change view to top view +view([0,90]) + +% create axes corresponding to the chosen map projection +axh = axesm(sArgs.proj, ... + 'MapLatLim', latLim, ... + 'MapLonLim', lonLim, ... + 'Grid','on', ... + 'Frame','on', ... + 'LabelFormat','none', ... + 'PLineLocation',30, ... + 'MLineLocation',50, ... + 'Origin',[0,0]); +axh.YLim = latLim * pi/180; +axh.XLim = lonLim * pi/180; + +% set ticks +axh.YTick = YTicks; +axh.YTickLabels = LatTicks; +axh.XTick = XTicks; +axh.XTickLabels = LonTicks; + +axh.YTickLabel = strcat(axh.YTickLabel, '^\circ'); +axh.XTickLabel = strcat(axh.XTickLabel, '^\circ'); + +% set labels for x- and y-axes +axh.YLabel.String = 'Latitude'; +axh.XLabel.String = 'Longitude'; + +if nargout + varargout{1} = sArgs.fgh; +end + +%end function +end \ No newline at end of file diff --git a/kernel/ClassStuff/@itaCoordinates/surf.m b/kernel/ClassStuff/@itaCoordinates/surf.m index 734ccc82f12ddb5d46ca02aeb7b8e85846306c80..b210030979352fc41c5ce852ffd62b88f3aeacb8 100644 --- a/kernel/ClassStuff/@itaCoordinates/surf.m +++ b/kernel/ClassStuff/@itaCoordinates/surf.m @@ -38,16 +38,19 @@ function varargout = surf(this, varargin) % Author: % Martin Pollow, mpo@akustik.rwth-aachen.de, 4.1.2010 -narginchk(1,inf); -defaultProperties = {'EdgeAlpha',0, 'FaceColor', 'interp'}; -sArgs = struct('pos1_data','double', 'radius', [],'magnitude',0,defaultProperties{:},'parent',[]); -if ~isempty(varargin) - [data,sArgs] = ita_parse_arguments(sArgs,varargin); -end +defaultProperties = {'EdgeAlpha', 0, ... + 'FaceColor', 'interp'}; -numArgIn = nargin; +sArgs = struct('pos1_data','double', ... + 'radius', [], ... + 'magnitude', 0, ... + 'parent', [], ... + 'complex', false, ... + 'colorbar', true, ... + defaultProperties{:}); +[data,sArgs] = ita_parse_arguments(sArgs,varargin); if sArgs.magnitude data = abs(data); @@ -64,7 +67,7 @@ end % parent = varargin{2}; % varargin = varargin(3:end); % end -if (numArgIn == 1) +if (nargin == 1) % only coordinates given r = this.r; color = zeros(size(r)); @@ -76,7 +79,7 @@ else else r = data.'; end - isComplex = ~all(isreal(r));% & min(r(:)) < 0; + isComplex = ~all(isreal(r)) || sArgs.complex; if isComplex % if the radius is complex, set the color color = mod(angle(r),2*pi); @@ -193,15 +196,19 @@ switch colorbar_settings case 'geometry' % do nothing case 'complex' - colormap hsv - colorbar; + colormap(hsv) + if sArgs.colorbar + colorbar; + end caxis([0 2*pi]); case 'magnitude' colormap jet % caxis([0 max(color)]); if min(color) < max(color) caxis([min(color) max(color)]); - colorbar; + if sArgs.colorbar + colorbar; + end elseif all(isfinite(color)) caxis([0 max(color)]); else diff --git a/kernel/DSP/Edit/ita_extend_dat.m b/kernel/DSP/Edit/ita_extend_dat.m index 4b89bb3970ace67079189c78801dd1ff52f4ba75..abec6d1435c03e82f033de66c869e549616af327 100644 --- a/kernel/DSP/Edit/ita_extend_dat.m +++ b/kernel/DSP/Edit/ita_extend_dat.m @@ -7,8 +7,8 @@ function [ varargout ] = ita_extend_dat( varargin ) % Syntax: ita_extend_dat( dat, nSamples, Options) % Syntax: [dat1, dat2] = ita_extend_dat( dat1, dat2 ) % -% FFTdegree is up to a value of 30 -% nSamples is a value greater than 30 +% FFTdegree is up to a value of (including) 35 +% nSamples is a value greater than 35 % % Options (default): % 'forcesamples' (false): Interpret second arguments as samples, even if lower than 35 @@ -59,7 +59,7 @@ if make_same_length asData2 = ita_metainfo_rm_historyline(asData2); else if sArgs.symmetric %pdi added - if new_number_samples < 35 + if new_number_samples <= 35 new_number_samples = round(2.^new_number_samples/2)*2; end old_number_samples = asData.nSamples; diff --git a/kernel/DSP/Edit/ita_normalize_dat.m b/kernel/DSP/Edit/ita_normalize_dat.m index 0e9e0f4ca249c9f235495ea9016a3879af870f19..31a10214facde8c2b00b36cde35ccd06005d0564 100644 --- a/kernel/DSP/Edit/ita_normalize_dat.m +++ b/kernel/DSP/Edit/ita_normalize_dat.m @@ -74,7 +74,7 @@ if sArgs.allchannels gainApplied = nan(1, result.nChannels); for ch_idx = 1:result.nChannels gainApplied(1, ch_idx) = max(abs(result.dat(ch_idx,:))); - result.dat(ch_idx,:) = result.dat(ch_idx,:) ./ gainApplied; + result.dat(ch_idx,:) = result.dat(ch_idx,:) ./ gainApplied(1, ch_idx); end else gainApplied = max(max(abs(result.timeData))); diff --git a/kernel/DSP/Edit/ita_normalize_spk.m b/kernel/DSP/Edit/ita_normalize_spk.m index 1b84ee7f2db2b1284319be6adb0741ae296ca951..240119ff4e2ec0043dfede7936b0170ec51fdad0 100644 --- a/kernel/DSP/Edit/ita_normalize_spk.m +++ b/kernel/DSP/Edit/ita_normalize_spk.m @@ -66,9 +66,10 @@ end %% Normalize if sArgs.allchannels + gainApplied = nan(1, result.nChannels); for ch_idx = 1:result.nChannels - gainApplied = max(abs(result.freqData(:,ch_idx))); - result.freqData(:,ch_idx) = result.freqData(:,ch_idx) ./ gainApplied; + gainApplied(1, ch_idx) = max(abs(result.freqData(:,ch_idx))); + result.freqData(:,ch_idx) = result.freqData(:,ch_idx) ./ gainApplied(1, ch_idx); end else gainApplied = max(max(abs(result.freqData))); diff --git a/kernel/DSP/Filter/ita_filter_fractional_octavebands.m b/kernel/DSP/Filter/ita_filter_fractional_octavebands.m index cf565f261c464944cff326998de4d0f3c406d532..d4014e2c02cbadc1dcf334f0fce1b9702104d50f 100644 --- a/kernel/DSP/Filter/ita_filter_fractional_octavebands.m +++ b/kernel/DSP/Filter/ita_filter_fractional_octavebands.m @@ -41,6 +41,8 @@ input = sArgs.data; %% Call mpb filter input = ita_mpb_filter(input,'oct',sArgs.bandsperoctave,'octavefreqrange',sArgs.freqRange,'zerophase',sArgs.zerophase,'order',sArgs.order); +% units are not copied in mpf_filter routine, do it here (?) +input.channelUnits(:) = input.channelUnits(1); if nargin == 0 % setin base diff --git a/kernel/DataAudio_IO/ita_read/ita_read_flac.m b/kernel/DataAudio_IO/ita_read/ita_read_flac.m new file mode 100644 index 0000000000000000000000000000000000000000..be6a05ec900b739f793e09698d40b134e91817f0 --- /dev/null +++ b/kernel/DataAudio_IO/ita_read/ita_read_flac.m @@ -0,0 +1,37 @@ +function result = ita_read_flac(filename,varargin) +%ITA_READ_FLAC - Read Free Lossless Audio Codec Files +% This function is completely based on the MATLAB audioread. +% +% It returns a itaAudio object containing the files data and metadata. +% +% See also ita_read, ita_write, audioread. + +% +% 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. +% + + +%% Return type of data this function can read +if nargin == 0 + result{1}.extension = '.flac'; + result{1}.comment = 'FLAC Files (*.flac)'; + return +end + + +if ~exist(filename,'file') + error('ITA_READ_FLAC: File does not exist'); +end + +try + [data,fs] = audioread(filename); + result = itaAudio; + result.samplingRate = fs; + result.timeData = data; +catch + error('ITA_READ_FLAC: Something went wrong'); +end + +end + diff --git a/kernel/DataAudio_IO/private/getDocNode.m b/kernel/DataAudio_IO/private/getDocNode.m index a5619a77abb01d45a39ef02f0c678e50065098f6..39cf0527ca1ed6bbc03ec184a5dd5d13930f839b 100644 --- a/kernel/DataAudio_IO/private/getDocNode.m +++ b/kernel/DataAudio_IO/private/getDocNode.m @@ -3,7 +3,11 @@ function docNode = getDocNode(transFunc, inputImp, bsName) xml = getXMLStrings(); docNode = com.mathworks.xml.XMLUtils.createDocument(xml.DocType); +domImpl = docNode.getImplementation(); +doctype = domImpl.createDocumentType(xml.DocType, [], ''); +docNode.appendChild(doctype); docRootNode = docNode.getDocumentElement; + % docRootNode.setAttribute('attr_name','attr_value'); versionNode = docNode.createElement(xml.Version); @@ -11,6 +15,9 @@ versionNode.appendChild... (docNode.createTextNode(xml.ActualVersion)); docRootNode.appendChild(versionNode); +optionsNode = getOptionsNode(docNode, bsName); +docRootNode.appendChild(optionsNode); + mainNode = docNode.createElement(xml.SaveTypeQuadripole); docRootNode.appendChild(mainNode); @@ -26,10 +33,9 @@ xml = getXMLStrings(); quadMSNode = docNode.createElement(xml.QuadripoleMeasured); -optionsNode = getOptionsNode(docNode, bsName); [freqVecNode, transferFunctionNode, inputImpedanceNode] = lsData2DomNodes(docNode, transFunc, inputImp); -quadMSNode.appendChild(optionsNode); + quadMSNode.appendChild(freqVecNode); quadMSNode.appendChild(transferFunctionNode); quadMSNode.appendChild(inputImpedanceNode); diff --git a/kernel/PlayrecMidi/ita_playrec_init.m b/kernel/PlayrecMidi/ita_playrec_init.m new file mode 100644 index 0000000000000000000000000000000000000000..360044dd009608737d497b930c2cb508093f6a27 --- /dev/null +++ b/kernel/PlayrecMidi/ita_playrec_init.m @@ -0,0 +1,54 @@ +function ita_playrec_init(varargin) +%ITA_PLAYREC_INIT - Initialize Playrec +% This function initializes creates a handle for playrec and initialzes +% playrec and PortAudio. +% +% Syntax: +% ita_playrec_init(options) +% +% Options: +% 'handle' : handle to playrec +% 'samplingRate' : sampling rate of the audio interface +% 'recDeviceID' : portaudio device id of the recording device +% 'playDeviceID' : portaudio device id of the playback device +% +% +% Example: +% ita_playrec_init() +% +% See also: +% ita_portaudio, ita_portaudio_run, ita_playrec +% +% Reference page in Help browser +% doc ita_playrec_init + +% +% 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. +% + + +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: 08-May-2017 + +sArgs = struct('handle', [], ... + 'samplingRate', ita_preferences('samplingRate'), ... + 'playDeviceID', ita_preferences('playDeviceID'), ... + 'recDeviceID', ita_preferences('recDeviceID')); +sArgs = ita_parse_arguments(sArgs, varargin); + +if isempty(sArgs.handle) + hPlayRec = ita_playrec; +else + hPlayRec = sArgs.handle; +end + +if ~hPlayRec('isInitialised') + hPlayRec('init', sArgs.samplingRate, sArgs.playDeviceID, sArgs.recDeviceID); + ita_verbose_info('Initializing Playrec... waiting 1 second...', 0); + pause(1); +else + ita_verbose_info('Playrec is already initialized.',0); +end +%end function +end \ No newline at end of file diff --git a/kernel/PlotRoutines/Plottools/ita_plottools_cursors.m b/kernel/PlotRoutines/Plottools/ita_plottools_cursors.m index 6e0b6305c6108174274a0ef298a8d7f28c9cb962..ada855fb5d2e99a7e533d2dbc400074c26ed7e88 100644 --- a/kernel/PlotRoutines/Plottools/ita_plottools_cursors.m +++ b/kernel/PlotRoutines/Plottools/ita_plottools_cursors.m @@ -138,7 +138,7 @@ switch state 'Parent',current_axh,... 'UserData',lnh,... 'Visible','off'); %do not display cursor initially - %[yl(1) yl(2)] + ph2 = line([xv2 xv2],[-1e7 1e7], ... 'Color',PassiveColor, ... 'Marker',marker, ... @@ -150,6 +150,11 @@ switch state 'UserData',lnh,... 'Visible','off'); %do not display cursor initially + % hide the cursors from the plot legend + set(get(get(ph1,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); + set(get(get(ph2,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); + + setappdata(ph1,'sample_CurrentPos',x_init(1)); setappdata(ph2,'sample_CurrentPos',x_init(2)); diff --git a/kernel/ita_path_handling.m b/kernel/ita_path_handling.m index 6e91e365b6832772d028aa07ddcfb7dd4f9b1e65..b3cd89896bf0eb1017bdd50e4856a8b361ffb046 100644 --- a/kernel/ita_path_handling.m +++ b/kernel/ita_path_handling.m @@ -18,7 +18,7 @@ function varargout = ita_path_handling(varargin) if nargin > 0 error('There should not be any input arguments to ita_path_handling'); end -ignoreList = {'.git','.svn','private','tmp','prop-base','props','text-base','template','doc'}; +ignoreList = {'.git','.svn','private','tmp','prop-base','props','text-base','template','doc','helpers'}; %% toolbox prefix string fullpath = ita_toolbox_path(); diff --git a/tutorials/helpers/plot_basis_functions.m b/tutorials/helpers/plot_basis_functions.m new file mode 100644 index 0000000000000000000000000000000000000000..6a2cdcfeb9eb8d374b66205bb1501684b71f1a95 --- /dev/null +++ b/tutorials/helpers/plot_basis_functions.m @@ -0,0 +1,38 @@ +function plot_basis_functions(sampling, Nmax, Y) +% Helper function for plotting the SH basis functions +% +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: Jun-2017 + +% +% This file is part of the ITA-Toolbox. All rights reserved. +% You can find the license for this m-file in the application folder. +% + + +spatVec = zeros(sampling.nPoints, (Nmax + 1)^2); +for idx = 1:(Nmax + 1)^2 + shVec = zeros((Nmax + 1)^2, 1); + shVec(idx) = 1; + spatVec(:, idx) = Y(:,1:(Nmax+1)^2) * shVec; +end + + +mPlot = numel(-Nmax:Nmax); +nPlot = Nmax+1; +figure; + +for n = 0:Nmax + idxmPlot = ceil(mPlot / 2); + for m = -n:n + nm = ita_sph_degreeorder2linear(n, m); + subplot(nPlot,mPlot, (n * mPlot) + idxmPlot + m) + surf(sampling,spatVec(:,nm),'complex',~isreal(Y)); + view([140,21]) + if (n == 0) + title('Spherical harmonic basis functions') + end + end +end + +end \ No newline at end of file diff --git a/tutorials/ita_tutorial_itaCoordinates.m b/tutorials/ita_tutorial_itaCoordinates.m index d34393b34533a4196cfefc8844c2557bec1128ca..b881997ea87b9e618a9687ba3f7c90fd679aa0f7 100644 --- a/tutorials/ita_tutorial_itaCoordinates.m +++ b/tutorials/ita_tutorial_itaCoordinates.m @@ -128,11 +128,11 @@ surf(coord, data); % chosen. unitSphere = ones(coord.nPoints,1); -surf(coord, unitSphere, data); +surf(coord, data,'radius',unitSphere); % note that now the color resembles the magnitude of the data % by the way, the short form does the same: -surf(coord, 1, data) +surf(coord, data,'radius',1) %% Plotting phase details of the directivity @@ -141,7 +141,7 @@ surf(coord, 1, data) % convert the matlab angle (-pi..pi) to positive phase (0..2pi) phase = angle(data) + pi; % phase on unit sphere -surf(coord, 1, phase); +surf(coord, phase,'radius',1); caxis([0 2*pi]); colormap hsv % here you have to adjust the color axis manually, as the surf plot @@ -160,7 +160,7 @@ surf(coord_half, data_half) % plop factor gives the relative size of the triangles that we destroy % together with the larger ones. -surf(coord_half, 1, data_half, 'plop', 4) +surf(coord_half, data_half, 'radius',1,'plop', 4) % The higher the plop factor, the more patches are preserved. %% Sample plot with lightning and additional parameters @@ -178,7 +178,7 @@ camlight % For better quality, the absolute maximum should be determined and the % same axes length be chosen. -return; %this is only for the publish() function +% return; %this is only for the publish() function % take a lower spatial resolution of 22.5°/22.5° (theta/phi) coordSmall = ita_generateSampling_equiangular(22.5,22.5); % create an audio object @@ -191,6 +191,6 @@ ao.freq = bsxfun(@times, 1e4./ao.freqVector, cardioid.'); ao.freq = ao.freq + rand(size(ao.freq)); % and plot the directivities of this audio object for some frequencies for freqs = 1000:1000:16000 - surf(coordSmall, ao, freqs) + surf(coordSmall, ao.freq2value(freqs)) pause(0.3); end diff --git a/tutorials/ita_tutorial_spherical_harmonics.m b/tutorials/ita_tutorial_spherical_harmonics.m new file mode 100644 index 0000000000000000000000000000000000000000..6c17a2d98ffebce09b18a25116476006fd2b08f0 --- /dev/null +++ b/tutorials/ita_tutorial_spherical_harmonics.m @@ -0,0 +1,136 @@ +%% Sound field synthesis and analysis in spherical harmonics +% +% <<../../pics/ita_toolbox_logo_wbg.png>> +% +% This tutorial demonstrates some of the spherical harmonic signal +% processing techniques implemented in the ITA-Toolbox +% +% For information about the presented methods, refer to +% Williams - Fourier Acoustics, +% Rafaely - Fundamentals of Spherical Array Processing + +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: Jun-2017 + +% +% This file is part of the ITA-Toolbox. All rights reserved. +% You can find the license for this m-file in the application folder. +% + +%% First off, clear the workspace +ccx +cd(fullfile(ita_toolbox_path, 'tutorials')) +addpath('helpers') + +%% Let's start by creating a spherical sampling of order N_max +NmaxSampling = 20; +sampling = ita_sph_sampling_gaussian(NmaxSampling); + +% take a look at the sampling +sampling.scatter; + +% we can calculate the complex valued spherical harmonic basis functions +% for the sampling with +Y = ita_sph_base(sampling,NmaxSampling); +Y_real = ita_sph_base(sampling, NmaxSampling, 'real'); + +% plot the basis functions +plot_basis_functions(sampling, 2, Y) +plot_basis_functions(sampling, 2, Y_real) +%% Indexing +% The linear index nm for spherical harmonic coefficients is defined as +% nm = n.^2 + n + m + 1 +% Functions exist to convert from one indexing to the other and vice versa +n = 2; +m = 1; +nm = ita_sph_degreeorder2linear(n, m) +[n, m] = ita_sph_linear2degreeorder(nm) + + + +%% Plane wave incident +Nmax = 7; +real = true; +k = linspace(0.5, 5, 128); + +DOA = itaCoordinates([1,0,0]); +y_vec = ita_sph_base(DOA, Nmax, 'real', real)'; + +% calculate the modal strength on a rigid sphere for a plane wave incidence +% from the DOA +B = ita_sph_modal_strength(sampling, Nmax, k, 'rigid'); + +% calculate the resulting spherical harmonic coefficients +pnm = zeros((Nmax + 1)^2, numel(k)); +for idx = 1:numel(k) + pnm(:,idx) = B(:,:,idx) * y_vec; +end +%% Let's look at the pressure on the sphere +% The wave number index for plotting can be chosen by changing idxPlot +idxPlot = 50; +figure +sampling.surf(pnm(:,idxPlot)) +title(['kr = ', num2str(k(idxPlot) * uniquetol(sampling.r))]); + +%% Beamforming + +% define the look directions of the beamformer +lookDirections = ita_sph_sampling_equiangular(30); + +% calculate the sh basis functions for all look directions +Y_lookDirs = ita_sph_base(lookDirections, Nmax)'; + +% calculate the radial filter to compensate for the modal strength of the +% sphere +radialFilter = zeros((Nmax+1)^2, (Nmax+1)^2, numel(k)); +for idx = 1:numel(k) + radialFilter(:,:,idx) = pinv(B(:,:,idx)); +end + +% perform the beamforming +beamformerOutput = zeros(lookDirections.nPoints, numel(k)); +for idx = 1:numel(k) + beamformerOutput(:,idx) = 4*pi / (Nmax + 1)^2 * Y_lookDirs' * radialFilter(:,:,idx) * pnm(:,idx); +end + +% plot a map projection for all look directions used in the beamforming +figure +lookDirections.plot_map_projection(beamformerOutput(:,1)) +%% Sound pressure from a vibrating spherical cap at the north pole +% coordinates of the north pole +northPole = itaCoordinates([0,0,1]); + +NmaxCap = 20; +% SH basis functions for the north pole +Y_cap = ita_sph_base(northPole, NmaxCap); + + +% apterture angle of 30 degree for the spherical cap +alpha = 30 * pi / 180; +rCap = sin(alpha); +G = ita_sph_aperture_function_sla(northPole, NmaxCap, rCap, 'diag'); + +% radial distance from the sphere +distance = 3; +% radiation impedance and propagation term from the sphere to a point in +% free space with a distance of 3 meters from the origin +H = ita_sph_modal_strength(northPole, NmaxCap, k, 'rigid', 'transducer', 'ls', 'dist', distance); + +% radial velocity of the cap +u = 0.01; + +% grid around the sphere for plotting purposes +grid = ita_sph_sampling_equiangular(20); +grid_Y = ita_sph_base(grid, NmaxCap); + +% resulting sound pressure on the specified grid in 3 meter distance +p = zeros(grid.nPoints, numel(k)); +for idx = 1:numel(k) + p(:,idx) = conj(grid_Y) * H(:,:,idx) * G * Y_cap' * u; +end + +% plot the resulting sound pressure at a spherical surface in the specified +% distance +idxPlot = 50; +surf(grid, p(:,idxPlot)); +title(['k = ', num2str(k(idxPlot) * uniquetol(northPole.r))]);