Skip to content

GitLab

toolbox
toolbox
Commit f95c0917 authored by Jan-Gerrit Richter's avatar Jan-Gerrit Richter
Browse files
Options

added new arc postprocess function for both step wise and continuous measurement 

delete old arc functions (all by zillekens (2013)
parent 8adabeec
Hide whitespace changes
Inline Side-by-side
Showing with 709 additions and 746 deletions
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_coord_flute_mics.m deleted 100644 → 0
View file @ 8adabeec
function [ coord_mics ] = ita_HRTFarc_coord_flute_mics( phi, r_z, d, a, h, gamma, beta, alpha,xOff,yOff)
%TEST_ZILLEKENS_COORD_FLUTE_MICS returns an itaCoordinates containing the position of the
%flutes mics
%
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d)
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d, a)
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d, a, h)
% coord_mics = TEST_ZILLEKENS_COORD_FLUTE_MICS(phi, r_z, d, a, h, gamma, beta, alpha)
%
% phi - vector of the angles phi
% r_z - distance of the mics to center of flute
% d - length of boom arm
%
% optional:
% a - distance between joint and boom arm
% h - height of stand to joint in center of turntable (default 0)
%
% optional uncertainties:
% gamma - angle between mic holder and flute (default 0)
% beta - angle between mic stand and mic boom arm (default 0)
% alpha - angle offset of reference point of turntable (default 0)
%
% See also TEST_ZILLEKENS_HT_DHTRAFO TEST_ZILLEKENS_ERROR_DISTANCE_LS
% Author: Stefan Zillekens
% Created: 2013-09-02
if(nargin < 6)
gamma = 0; % angle beetween mic holder and flute
beta = 0; % angle beetween mic stand and mic boom
alpha = 0; % angle offset of reference point of motor
end
if(nargin < 5)
h = 0;
end
if(nargin < 4)
a = 0;
end
% init
npos = numel(phi);
nmic = numel(r_z);
coord_mics = itaCoordinates(npos*nmic);
XYZ = zeros(npos*nmic,3);
for idp = 1:npos;
for idm = 1:nmic;
% homogeneous transformation matrices:
T1 = ita_HRTFarc_ht_DHtrafo(-pi/2+phi(idp)+alpha, h, 0, beta);
T2 = ita_HRTFarc_ht_DHtrafo(+pi/2,a,d,gamma);
T3 = ita_HRTFarc_ht_DHtrafo(0, r_z(idm), 0, 0);
T = T1 * T2 * T3;
XYZ((idp-1)*nmic+idm, :) = T(1:3,4)';
end
end
coord_mics.cart = XYZ;
coord_mics.x = coord_mics.x+xOff;
coord_mics.y = coord_mics.y+yOff;
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_cropMeasurementData.m deleted 100644 → 0
View file @ 8adabeec
function [ output_args ] = ita_HRTFarc_cropMeasurementData( iMS,dataPath )
%ITA_HRTFARC_CROPMEASUREMENTDATA crop of raw measurement data
% This function is intended to be used after a measurement with the HRTF
% arc. During the measurement, the cropping is not done to save time.
% The function assumes a data structure from itaEimar
% It will read all files from dataPath/data, crop them
% It will rename the folder to data_raw. Cropped data will be in
% dataPath/data to maintain backwards compability.
savePath = sprintf('%s/cropped/',dataPath);
dataPathInternal = sprintf('%s/data/',dataPath);
mkdir(savePath)
% load all data, crop and save again
dataDir = dir(dataPathInternal);
wb = itaWaitbar(length(dataDir)-2);
for index = 3:length(dataDir)
data = ita_read(sprintf('%s/%s',dataPathInternal,dataDir(index).name));
data_cropped = iMS.crop(data);
ita_write(data_cropped,sprintf('%s/%s',savePath,dataDir(index).name));
wb.inc;
end
wb.close;
rawPath = sprintf('%s/data_raw/',dataPath);
movefile(dataPathInternal,rawPath);
movefile(savePath,dataPathInternal);
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_error_distance_LS.m deleted 100644 → 0
View file @ 8adabeec
function [ err ] = ita_HRTFarc_error_distance_LS( x, toa, phi, r_z, SR )
%TEST_ZILLEKENS_ERROR_DISTANCE_LS returns the result of a least-mean square
%error
%
% err = TEST_ZILLEKENS_ERROR_DISTANCE_LS( x, toa, phi, r_z )
% err = TEST_ZILLEKENS_ERROR_DISTANCE_LS( x, toa, phi, r_z, SR )
%
% x - containg a vector of LS coordinates, latency sample between DA
% converters, alpha, beta, gamma, d, height, speed of sound
% toa - vector of time of arrivals
% phi - vector of phi position of flute mics
% r_z - vector with distances of mics to center of flute
% SR - sampling rate
%
% See also TEST_ZILLEKENS_COORD_FLUTE_MICS TEST_ZILLEKENS_OPTIMIZE_COORD_ARCLS
% Author: Stefan Zillekens
% Created: 2013-09-09
if nargin < 5
SR = 44100;
end
nopt = 8; % number of values to optimize
npos = numel(phi);
nmic = numel(r_z);
nLS = (numel(x)-nopt)/3;
if numel(toa) ~= nLS*npos*nmic;
error('Dimensions do not fit')
end
%% uncertainties:
a = x(end-7); % distance from joint to boom arm ~4cm
d = x(end-6); % length of boom arm
% height = x(end-4);
c_meas = x(end-5); % Speed of Sound
alpha = x(end-4); % angle offset of reference point of TT ~0
% alpha = 0;
beta = x(end-3); % angle between mic stand and mic boom ~pi/2
gamma = x(end-2); % angle between mic holder and flute ~0
xOff = x(end-1);
yOff = x(end);
% latency = x(end-4:end); % relative changes in latency of D/A converter
latency = zeros(1,8);
%% calculate mic coordinations via homogeneous transforms
height = 0;
coord_mics = ita_HRTFarc_coord_flute_mics(phi,r_z,d,a,height,gamma,beta,alpha,xOff,yOff);
c_meas = 344;
%% error calculation
micsCart = coord_mics.cart;
% tao MxLSxP
toa_reshape = reshape(toa,nmic,nLS,npos);
% only for the measurements from 28/29.05
% toa_reshape(:,8:8:32,:) = 0;
% toa_reshape(:,19:20,:) = 0;
toa_corrected = reshape(toa_reshape,1,nmic*nLS*npos);
%p1 3xMxLSxP
micPos = reshape(micsCart.',3,nmic,1,npos);
micPos = repmat(micPos,[1,1,nLS,1]);
micPos = reshape(micPos,3,nLS*npos*nmic);
%p2 3xMxLSxP
lsPos = reshape(x(1:end-nopt),3,nLS);
lsPos = repmat(lsPos,[1,1,nmic,npos]);
lsPos_reshape = permute(lsPos,[1 3 2 4]);
lsPos = reshape(lsPos_reshape,3,nLS*npos*nmic);
diffMatrix = lsPos-micPos;
dist = sqrt(diffMatrix(1,:).^2+diffMatrix(2,:).^2+diffMatrix(3,:).^2);
dist_reshape = reshape(dist,nmic,nLS,npos);
% dist_reshape(:,8:8:32,:) = 0;
% dist_reshape(:,19:20,:) = 0;
dist_corrected = reshape(dist_reshape,1,nmic*nLS*npos);
err = (dist_corrected-toa_corrected/SR*c_meas).';
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_ht_DHtrafo.m deleted 100644 → 0
View file @ 8adabeec
function [ T ] = ita_HRTFarc_ht_DHtrafo( theta, d, a, alpha )
%TEST_ZILLEKENS_HT_DHTRAFO Homogeneous transformation matrix via
%Denavit-Hartenberg (DH) parameter
%
% T = TEST_ZILLEKENS_HT_DHTRAFO( theta, d, a, alpha )
%
% Rot('z_-1',theta) * Trans([0,0,d]) * Trans([a,0,0]) * Rot('x',aplpha)
%
% Rot(axis,angle) rotation around axis of coordinate system with angle
% Trans([x,y,z]) translation in x, y, z direction.
%
% angles in radiant
%
% EXAMPLE:
% T = test_zillekens_ht_DHtrafo( -pi, 1, 0, pi/8 )
%
% See also TEST_ZILLEKENS_HT_TRANS TEST_ZILLEKENS_HT_ROT
% Author: Stefan Zillekens
% Created: 2013-07-30
% References:
% [1] Roßmann, Jürgen: Mensch-Maschine-Interaktion und Robotik I,
% Institut für Mensch-Maschine-Interaktion, Aachen, 2007
T = ita_HRTFarc_ht_rot('z',theta) ...
* ita_HRTFarc_ht_trans([0,0,d]) ...
* ita_HRTFarc_ht_trans([a,0,0]) ...
* ita_HRTFarc_ht_rot('x',alpha);
%%
% T = zeros(4);
% T(1:2,1) = [ cos(theta) sin(theta) ];
% T(1:3,2) = [ -cos(alpha)*sin(theta) cos(alpha)*cos(theta) sin(alpha) ];
% T(1:3,3) = [ sin(alpha)*sin(theta) -sin(alpha)*cos(theta) cos(alpha) ];
% T(:,4) = [ a*cos(theta) a*sin(theta) d 1 ];
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_ht_rot.m deleted 100644 → 0
View file @ 8adabeec
function [ T_rot ] = ita_HRTFarc_ht_rot( rot_axis, phi )
%TEST_ZILLEKENS_HT_ROT Homogenous transformation matrix for a rotation
%
% T_rot = TEST_ZILLEKENS_HT_TRANS( 'rot_axis', angle )
%
% rot_axis - rotatation axis { 'x', 'y', 'z' }
% phi - rotation angle in radiant
%
% EXAMPLE:
% T_rot = test_zillekens_ht_rot( 'z', pi )
% % Rotation around z-axis with the angle pi
%
% See also TEST_ZILLEKENS_HT_TRANS TEST_ZILLEKENS_HT_DHTRAFO
% Author: Stefan Zillekens
% Created: 2013-07-30
% References:
% [1] Roßmann, Jürgen: Mensch-Maschine-Interaktion und Robotik I,
% Institut für Mensch-Maschine-Interaktion, Aachen, 2007
if ~isscalar(phi)
error('Expecting a scalar')
end
T_rot = eye(4);
if rot_axis == 'x'
T_rot(2:3,2:3) = [ cos(phi), -sin(phi); ...
sin(phi), cos(phi) ];
elseif rot_axis == 'y'
T_rot(1:2:3,1:2:3) = [ cos(phi), sin(phi); ...
-sin(phi), cos(phi) ];
elseif rot_axis == 'z'
T_rot(1:2,1:2) = [ cos(phi), -sin(phi); ...
sin(phi), cos(phi) ];
else
error('Rotation axis should be x, y or z.')
end
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_ht_trans.m deleted 100644 → 0
View file @ 8adabeec
function [ T_trans ] = test_rbo_zillekens_ht_trans( coord )
%TEST_ZILLEKENS_HT_TRANS Homogenous transformation matrix for a translation
%
% Ttrans = TEST_ZILLEKENS_HT_TRANS( coord )
% Ttrans = TEST_ZILLEKENS_HT_TRANS( itaCoordinates )
%
% EXAMPLE:
% Ttrans = test_zillekens_ht_trans( [1, 0, 0] )
% % Translation of 1 in x direction
% Ttrans = test_zillekens_ht_trans( itaCoordinates )
%
% See also TEST_ZILLEKENS_HT_ROT TEST_ZILLEKENS_HT_DHTRAFO
% Author: Stefan Zillekens
% Created: 2013-08-22
% References:
% [1] Roßmann, Jürgen: Mensch-Maschine-Interaktion und Robotik I,
% Institut für Mensch-Maschine-Interaktion, Aachen, 2007
% convert from itaCoordinates
% if isa(coord, 'itaCoordinates')
% coord = coord.cart;
% end
%
% coord = coord(:);
if ~isvector(coord) || length(coord)~= 3;
error('Invalid size of input data. Expecting [x,y,z] or an itaCoordinate')
end
T_trans = eye(4);
T_trans(1:3,4) = coord;
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_minimum_phase_IR.m deleted 100644 → 0
View file @ 8adabeec
function [ ao_mp ] = ita_HRTFarc_minimum_phase_IR( ai )
%MINIMUM_PHASE_IR returns minimum-phase IR using Hilbert transformation
%
% ao_mp = MINIMUM_PHASE_IR( itaAudio )
% Author: Stefan Zillekens
% Created: 2013-06-19
% check the input
if ~isa(ai, 'itaAudio')
error('Expecting an itaAudio.')
end
% hilbert transform
X = hilbert(-log(abs(ai.freqData)));
% minimal phase
mp = imag(X);
% minimal-phase version of ai
ao_mp = ai;
ao_mp.freqData = abs(ai.freqData) .* exp(1i*mp);
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_optimize_coord_arcLS.m deleted 100644 → 0
View file @ 8adabeec
function [ opt_coord_LS, opt_val ] = ita_HRTFarc_optimize_coord_arcLS( toa, phi, r_z, SR,coordLS_ideal)
%TEST_ZILLEKENS_OPTIMIZE_COORD_ARCLS returns the optimized arc LS
%coordinates of the arc and some optimized values
%
% [opt_coord_LS, opt_val] = TEST_ZILLEKENS_OPTIMIZE_COORD_ARCLS(toa, phi, r_z)
%
% INPUT:
% toa - time of arrival samples in a vector
% phi - vector of measured phi angles in radiant
% r_z - distances between the mics on the flute
% SR - sampling rate
%
% OUTPUT:
% opt_coord_LS - itaCoordinates with LS positions
% opt_val - struct with following values:
% a - distance to boom arm
% d - length of boom arm
% c - Speed of Sound
% alpha - offset angle of turntable
% beta - angle beetween mic stand and mic boom arm
% gamma - angle between mic holder and mic boom arm
% latency - latency samples between the 5 DA converters
%
% Origin of the coordinate system is considered in the center of the
% turntable (rotary axis) on top of the mic stands joint.
%
% See also TEST_ZILLEKENS_ERROR_DISTANCE_LS TEST_ZILLEKENS_COORD_HRTFARC
% TEST_ZILLEKENS_COORD2VEC
% Author: Stefan Zillekens
% Created: 2013-09-09
% Jan Richter - 27.05.2015
% renamed and changed to work with new arc setup
if nargin < 4
SR = 44100;
end
%% init
npos = numel(phi);
nmic = numel(r_z);
%% optimize those parameters [ start value, lower limit, upper limit ]
a_opt = [ 0.04 0.03 0.05 ];% distance from joint to boom arm ~4cm
d_opt = [ 0.4 0.3 0.5 ];% length of boom arm
% h_opt = [-0.005 0 0.005]; % height offset from the arm
c_opt = [ 345 320 350 ];% Speed of Sound
al_opt = [ 0 -pi/2 pi/2 ];% angle offset of reference point of TT ~0
be_opt = [ 0 -pi/2 pi/2 ]; % angle between mic stand and mic boom ~pi/2
ga_opt = [ 0 -pi/6 pi/6 ];% angle between mic holder and flute ~0
x_opt = [ 0.01 -0.05 0.05 ];% x offset of the drehteller
y_opt = [ 0.01 -0.05 0.05 ];% y offset of the drehteller
opt_val = struct( 'a', a_opt(1), ...
'd', d_opt(1), ...
'c', c_opt(1), ...
'alpha', al_opt(1), ...
'beta', be_opt(1), ...
'gamma', ga_opt(1), ...
'x', x_opt(1), ...
'y', y_opt(1));
% create ideal arc coordinates
nLS = coordLS_ideal.nPoints;
v_coordLS_init = reshape(coordLS_ideal.cart',1,[]); % to row vector
v_coordLS_min = min(coordLS_ideal.cart)-[0.3, 0.3, 0.8]; % Lower limit
v_coordLS_max = max(coordLS_ideal.cart)+[0.3, 0.3, 0.8]; % Upper limit
opt = optimset( 'TolFun', 1e-8, ...
'TolX', 1e-4, ...
'Display', 'iter', ...
'MaxIter', 40 ...
... 'FinDiffRelStep', [repmat(1e-3, 1, nLS*3), 1e-3, 1e-3, 1, 1e-4, 1e-4, 1e-4, ones(1,5)] ...
);
%% little warning
nvar = numel(fieldnames(opt_val)); % +4 caused by latency
if ( nvar+3*nLS > nLS*npos*nmic)
warning('There are more unknowns than equations.');
end
%% optimize
disp(['Starting optimizer at ',datestr(now,'HH:MM:SS')]);
% tic;
[x, ~] = lsqnonlin(@(x) ita_HRTFarc_error_distance_LS(x, toa, phi, r_z, SR ), ...
[v_coordLS_init a_opt(1) d_opt(1) c_opt(1) al_opt(1) be_opt(1) ga_opt(1) x_opt(1) y_opt(1)], ... % Startvalue
[repmat(v_coordLS_min, 1, nLS) a_opt(2) d_opt(2) c_opt(2) al_opt(2) be_opt(2) ga_opt(2) x_opt(2) y_opt(2)], ... % Lower limit
[repmat(v_coordLS_max, 1, nLS) a_opt(3) d_opt(3) c_opt(3) al_opt(3) be_opt(3) ga_opt(3) x_opt(3) y_opt(3)], ... % Upper limit
opt);
% toc;
disp(['Ended optimizer at ',datestr(now,'HH:MM:SS')]);
%% multi start optimization
% opts = optimoptions(@lsqnonlin,'Display','iter');
% problem = createOptimProblem('lsqnonlin','objective',...
% @(x) test_zillekens_error_distance_LS(x, toa, phi, r_z, SR ) ...
% ,'x0',[v_coordLS_init a_opt(1) d_opt(1) c_opt(1) al_opt(1) be_opt(1) ga_opt(1) x_opt(1) y_opt(1)], ...
% 'lb',[repmat(v_coordLS_min, 1, nLS) a_opt(2) d_opt(2) c_opt(2) al_opt(2) be_opt(2) ga_opt(2) x_opt(2) y_opt(2)], ...
% 'ub',[repmat(v_coordLS_max, 1, nLS) a_opt(3) d_opt(3) c_opt(3) al_opt(3) be_opt(3) ga_opt(3) x_opt(3) y_opt(3)], ...
% 'options',opts);
% ms = MultiStart('UseParallel',true);
% [x,f] = run(ms,problem,20)
opt_val.err = ita_HRTFarc_error_distance_LS(x, toa, phi, r_z, SR );
%% Solved
opt_val.a = x(end-7);
opt_val.d = x(end-6);
% opt_val.h = x(end-4);
opt_val.c = x(end-5);
opt_val.alpha = x(end-4);
opt_val.beta = x(end-3);
opt_val.gamma = x(end-2);
opt_val.x = x(end-1);
opt_val.y = x(end);
% Assign LS positions to itaCoordinates
v_posLS_opt = x(1:3*nLS);
m_posLS_opt = reshape(v_posLS_opt, 3,[])'; % row vector to matrix
opt_coord_LS = itaCoordinates(m_posLS_opt);
end
\ No newline at end of file
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_optimize_coords.m deleted 100644 → 0
View file @ 8adabeec
function coordArc = ita_HRTFarc_optimize_coords(varargin)
% based on test_zillikens_example_optimize_coords
% Jan Richter 27.05.2015
% This script process the data measured with the flute setup to obtain the
% coordinates of the HRTF arc. This script expects the data in itaAudio
% files, that were measured with the itaEimar in the HRTF arc. So there is
% a folder with one file for each azimuth angle. This azimuth angle is
% save in the itaAudio's channelCoordinates.
% The optimized arc's coordinates and some values are saved in one file
% called 'coordArc', which contains a struct with the optimized data.
%% properties
sArgs = struct('dataPath',[],'micPos',[0.2, 0.1, -0.1, -0.2],'radius',1.2,'deltaTheta',5.04);
sArgs = ita_parse_arguments(sArgs,varargin);
pathFlute = sArgs.dataPath; % flute path
rArc = sArgs.radius; % radius of the arc
deltaTheta = sArgs.deltaTheta; % distance of the LS in degree (elevation)
r_z = sArgs.micPos; % distance of the mircrophones to flute's middle point
nmic = numel(r_z); % number of used microphones
% pathFlute = 'D:\IRKO_HRTF\Fltenmessung\2015July21_1436_flute_5';
% shift_origin = false; % shift origin to arc's center point
%% check folder of data
dirlist = dir(pathFlute);
% kick out the non-directories and . and ..
for ind = numel(dirlist):-1:1
if ~dirlist(ind).isdir || dirlist(ind).name(1) == '.'
dirlist(ind) = []; % take out the non dirs
end
end
%% process the data
disp('__')
folder = dirlist.name;
filelist = dir([pathFlute filesep folder filesep '*.ita']); % check for .ita files
nFiles = numel(filelist);% it is a folder with .ita files in it: continue
npos = nFiles;
%% read files, merge and reshape
ita_disp(['Reading data: ', folder]);
ao = itaAudio(nFiles,1);
phi = zeros(nFiles,1);
for idf = 1 : nFiles
filename = [pathFlute filesep folder, filesep, filelist(idf).name];
dataFile = ita_read(filename);
ao(idf,1) = merge(dataFile);
phi(idf) = unique(ao(idf).channelCoordinates.phi);
end
nLS = numel(dataFile); % number of loudspeakers
ao = ao.merge;
aoReshape = ao;
% from order: nmic, nLS, npos
aoReshape.time = reshape(ao.time, [ ao.nSamples, nmic, nLS, npos]);
%% TOA
disp('Calculating TOA');
toa = ita_HRTFarc_toa_mcm(aoReshape);
%% calculations LS positions
SR = aoReshape.samplingRate;
% jri: this was really dangerous if more than 9 positions where taken
% phi = unique(aoReshape.channelCoordinates.phi);
thetaValues = 0:deltaTheta:(nLS-1)*deltaTheta;
coordTarget = itaCoordinates([rArc*ones(nLS,1), deg2rad(flipud(thetaValues')) zeros(nLS,1)],'sph');
disp('Optimize loudspeaker positions');
[coordArc, values] = ita_HRTFarc_optimize_coord_arcLS(toa, phi, r_z, SR, coordTarget);
end
applications/Binaural-HRTF/HRTFarc/ita_HRTFarc_postProcess.m 0 → 100755
View file @ f95c0917
function result = ita_HRTFarc_postProcess(varargin)
%ITA_HRTF_ARC_POSTPROCESS - +++ Short Description here +++
% This function ++++ FILL IN INFO HERE +++
%
% Syntax:
% audioObjOut = ita_HRTF_arc_postProcess(audioObjIn, options)
%
% Options (default):
% 'opt1' (defaultopt1) : description
% 'opt2' (defaultopt1) : description
% 'opt3' (defaultopt1) : description
%
% Example:
% audioObjOut = ita_HRTF_arc_postProcess(audioObjIn)
%
% See also:
% ita_toolbox_gui, ita_read, ita_write, ita_generate
%
% Reference page in Help browser
% <a href="matlab:doc ita_HRTF_arc_postProcess">doc ita_HRTF_arc_postProcess</a>
% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved.
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder.
% </ITA-Toolbox>
% Author: Jan-Gerrit Richter -- Email: jri@akustik.rwth-aachen.de
% Created: 05-Oct-2018
%% Initialization and Input Parsing
% all fixed inputs get fieldnames with posX_* and dataTyp
% optonal inputs get a default value ('comment','test', 'opt1', true)
% please see the documentation for more details
sArgs = struct( 'dataFolder','', ...
'refFile','',...
'ms',[], ...
'tw',[0.006 0.008], ...
'rotationDirection',-1, ...
'normalize',[], ...
'itdMethod','xcorr', ...
'dataChannel',1:2);
[options] = ita_parse_arguments(sArgs,varargin);
%% first, load and prepare the reference
ref = ita_read(options.refFile);
ref = merge(ref);
ref_tw = ita_time_window(ref,options.tw,'time');
% make end sample be div by 4 for daff export
endSample = round(options.tw(2) .* ref(1).samplingRate)+1;
endSample = endSample + mod(endSample,4);
ref_tw = ita_time_crop(ref_tw,[1 endSample],'samples');
ref_finished = ita_smooth_notches(ref_tw,'bandwidth',1/2,...
'threshold', 3);
clear tmp;
% create a multi instance again
for index = 1:ref_finished.nChannels
tmp(index) = ref_finished.ch(index);