Commit abee4652 authored by jme's avatar jme

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

parents 44eb987d 8cd5cb36
......@@ -19,6 +19,11 @@ function cThis = interp(this,varargin)
% set to 1 if no range extrapolation is required
% order ... order of spherical harmonics matrix (default: 50)
% epsilon ... regularization coefficient (default: 1e-8)
% shiftToEar to a shift to approximate ear position to
% improve sh transformation (see [2])
% shiftAxis shift along this axis ('x','y' (default),'z')
% shiftOffset shift ears (L - R) by these values
% (default: [-0.0725 0.0725])
%
% OUTPUT:
% itaHRTF object
......@@ -26,19 +31,24 @@ function cThis = interp(this,varargin)
% .timeData: interpolated / range-extrapolated HRIRs for defined field points
% .dirCoord: itaCoordinates object
%
% Required: SphericalHarmonics functions of ITA Toolbox
%
%
% [1] Pollow, Martin et al., "Calculation of Head-Related Transfer Functions
% for Arbitrary Field Points Using Spherical Harmonics Decomposition",
% Acta Acustica united with Acustica, Volume 98, Number 1, January/February 2012,
% pp. 72-82(11)
%
% [2] Richter, Jan-Gerrit et al. "Spherical harmonics based hrtf datasets:
% Implementation and evaluation for real-time auralization",
% Acta Acustica united with Acustica, Volume 100, Number 4, July/August 2014,
% pp. 667-675(9)
%
%
%
% Author: Florian Pausch <fpa@akustik.rwth-aachen.de>
% Version: 2016-02-05
% TODO: check why this is still not working (coordinate assignment???)
sArgs = struct('order',50,'eps',1e-5);
sArgs = struct('order',50,'eps',1e-5,'shiftToEar',false,'shiftAxis','y','shiftOffset',[-0.0725 0.0725]);
sArgs = ita_parse_arguments(sArgs,varargin,2);
if ~isa(varargin{1},'itaCoordinates'),error('itaHRTF:interp', ' An itaCoordinate object is needed!')
end
......@@ -47,7 +57,7 @@ field_in = varargin{1};
% only take unique direction coordinates (round to 0.01deg resolution)
tempfield = unique(round([field_in.phi_deg*100 field_in.theta_deg*100]),'rows'); % may cause problems with older Matlab versions (<=R2013)!
tempfield = tempfield./100;
temp_r = this.dirCoord.r(1);
temp_r = field_in.r(1);
field = itaCoordinates(size(tempfield,1));
field.r = repmat(temp_r,size(tempfield,1),1);
field.phi_deg = tempfield(:,1);
......@@ -67,46 +77,34 @@ if ~isequal(this.dirCoord.r(1),field.r(1))
kr0 = k*this.dirCoord.r(1); % measurement radius
kr1 = k*field.r(1); % extrapolation radius
hankel_r0 = ita_sph_besselh(1:Nmax,2,kr0);
hankel_r1 = ita_sph_besselh(1:Nmax,2,kr1);
hankel_r0 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr0);
hankel_r1 = ita_sph_besselh(ita_sph_linear2degreeorder(1:Nmax),2,kr1);
hankel_div = hankel_r1 ./ hankel_r0;
hankel_rep = hankel_div(:,1);
end
dweights = 1 + (0:Nmax).*((0:Nmax)+1); % calculate regularization weights
nSH = (Nmax+1).^2;
I = sparse(eye(nSH));
n = ita_sph_linear2degreeorder(1:round(nSH)).';
dweights_rep = zeros(sum(2*(0:Nmax)'+1),1);
dweights_rep(1)=dweights(1);
counter = 2;
for n=1:Nmax
nTimes = 2*n+1;
dweights_rep(counter:counter+nTimes-1)=dweights(n+1)*ones(nTimes,1);
if ~isequal(this.dirCoord.r(1),field.r(1))
hankel_rep=[hankel_rep, repmat(hankel_div(:,n),1,2*n+1)];
%% move data from earcenter
copyData = this;
if sArgs.shiftToEar
ear_d = sArgs.shiftOffset;
for ear=1:2
movedData = moveFullDataSet(this.ch(ear:2:this.nChannels),sArgs,ear_d(ear));
copyData.freqData(:,ear:2:copyData.nChannels) = movedData;
end
counter = counter + nTimes;
end
%% move data from earcenter
ear_d = [-0.07 0.07];
%% Weights
[~,w]= this.dirCoord.spherical_voronoi; % calculate weighting coefficients (Voronoi surfaces <-> measurement points)
% sG_full = ita_sph_sampling_equiangular(73,144,'theta_type','[]','phi_type','[)');
% % sG_full = ita_sph_sampling_gaussian(71);
% sG = sG_full;
% isOnArc = sG.theta < 2.75;
% sG.cart = sG.cart(isOnArc,:);
% sG.weights = sG.weights(isOnArc);
% % rescale the weights for the sum to be 4*pi again
% w = sG.weights .* 4.*pi ./ sum(sG.weights);
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,'real'); % calculate real-valued SHs using the measurement grid
D = I .* diag(1 + n.*(n+1)); % decomposition order-dependent Tikhonov regularization
Y = ita_sph_base(this.dirCoord,Nmax,'orthonormal',true); % calculate real-valued SHs using the measurement grid
%% Calculate HRTF data for field points
if Nmax > 25
......@@ -117,73 +115,58 @@ end
hrtf_arbi = zeros(this.nBins,2*field.nPoints); % columns: LRLRLR...
for ear=1:2
freqData_temp = this.freqData(:,ear:2:end);
% sh transformation
freqData_temp = copyData.freqData(:,ear:2:end);
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
% newCoords = this.dirCoord;
% newCoords.y = newCoords.y + ear_d(ear);
%
% data.channelCoordinates = newCoords;
% data.freqData = freqData_temp;
% data.freqVector = this.freqVector;
% data.c_meas = 344;
% data = process_result_delay_correction(data,this.dirCoord.r(1));
% % calculate weighted SH coefficients using a decomposition order-dependent Tikhonov regularization
%
% freqData_temp = data.freqData;
% Y = ita_sph_base(data.channelCoordinates,Nmax,'real');
% %% test the sh transformation results by plotting both spatial and sh data with surf
% s = itaSamplingSph(field);
% s.nmax = Nmax;
% figure
% surf(s,a0(:,10))
% figure
% surf(s,freqData_temp(10,:))
a0 = (Y.'*W*Y + epsilon*D) \ Y.'*W * freqData_temp.';
%a0 = (Y.'*W*Y + epsilon*D) \ Y.' *
%this.freqData(:,ear:2:end).'; % fpa version
%a0 = pinv(Y)* this.freqData(:,ear:2:end).'; % jck version
% range extrapolation
if ~isequal(this.dirCoord.r(1),field.r(1))
% calculate range-extrapolated HRTFs
a1 = a0 .* hankel_rep.';
Yest = ita_sph_base(field,N,'real'); % use real-valued SH's
%%% test here to see extrapolation results in spatial domain
% surf(s,a1(:,10))
% reconstruction to spatial data
Yest = ita_sph_base(field,N,'orthonormal',true); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a1).'; % interpolated + range-extrapolated HRTFs
else
Yest = ita_sph_base(field,Nmax,'real'); % use real-valued SH's
% reconstruction to spatial data
Yest = ita_sph_base(field,Nmax,'orthonormal',true); % use real-valued SH's
hrtf_arbi(:,ear:2:end) = (Yest*a0).'; % interpolated HRTFs
end
end
%% move back to head center
% todo
% for ear=1:2
%
% newCoords = field;
% newCoords.y = newCoords.y - ear_d(ear);
%
% freqData = hrtf_arbi(:,ear:2:end);
%
%
% data.channelCoordinates = newCoords;
% data.freqData = freqData;
% data.freqVector = this.freqVector;
% data.c_meas = 344;
% data = process_result_delay_correction(data,this.dirCoord.r(1));
%
% hrtf_arbi(:,ear:2:end) = data.freqData;
%
% end
% set new direction coordinates
sph = zeros(field.nPoints*2 ,3);
sph(1:2:end,:) = field.sph;
sph(2:2:end,:) = field.sph;
% write new HRTF data set
cAudio = itaAudio(hrtf_arbi, 44100, 'freq');
cAudio.channelCoordinates.sph= sph;
cThis = this;
cThis.freqData = hrtf_arbi;
cThis.channelCoordinates.sph= sph;
cThis = itaHRTF(cAudio);
cThis.freqData = hrtf_arbi;
%% move back to head center
if sArgs.shiftToEar
ear_d_back = -ear_d;
% movedData = zeros(size(hrtf_arbi));
for ear=1:2
movedData = moveFullDataSet(cThis.ch(ear:2:cThis.nChannels),sArgs,ear_d_back(ear));
cThis.freqData(:,ear:2:cThis.nChannels) = movedData;
end
end
if ~isequal(cThis.dirCoord.r(1),field.r(1))%???
cThis.dirCoord.r = field.r;
......@@ -195,15 +178,44 @@ end
end
function result = process_result_delay_correction(result, target_d)
% shift every measurement point to target_d by applying a
% phase shift to the channel: (simplified!)
freq_L = result.freqVector;
add_phase_L = (result.channelCoordinates.r - target_d)* (freq_L./result.c_meas)' .* 2.*pi;
function [ data ] = moveFullDataSet(data,options,offsetShift)
fullCoords = data.channelCoordinates;
freqVector = data.freqVector;
shiftedData = zeros(size(data.freqData));
axis = options.shiftAxis;
for index = 1:length(freqVector)
shiftedData(index,:) = moveHRTF(fullCoords,data.freqData(index,:),freqVector(index),axis,offsetShift);
end
data = shiftedData;
end
function [data] = moveHRTF(s, data, frequency, axis, offset)
% the offset is given in m
result.freqData = result.freqData .* exp(1i.*add_phase_L');
origAxis = s.r;
if (size(data,2) > size(data,1))
data = data.';
end
offset = real(offset); % ??
switch axis
case 'x'
s.x = s.x + offset;
case 'y'
s.y = s.y + offset;
case 'z'
s.z = s.z + offset;
end
result.channelCoordinates.r = target_d;
newAxis = s.r;
k = 2*pi*frequency/340;
% the phase is moved by the difference of the axis points
data = data .* exp(1i*k*(newAxis - origAxis));
% amplitude manipulation did not yield better results
% data = data .* newAxis ./ origAxis;
end
......@@ -71,12 +71,35 @@ assert( betarange( 2 ) >= 0.0 )
assert( betarange( 1 ) <= 180.0 )
%% Assemble metadata
%% Assemble metadata (if not already present)
keyname = 'Generation script';
if isempty(metadata) || ~any( strcmpi( { metadata(:).name }, keyname ) )
metadata = daffv17_add_metadata( metadata, keyname, 'String', 'writeDAFFFile.m' );
end
keyname = 'Generation toolkit';
if ~any( strcmpi( { metadata(:).name }, keyname ) )
metadata = daffv17_add_metadata( metadata, keyname, 'String', 'ITA-Toolkit' );
end
keyname = 'Generation date';
if ~any( strcmpi( { metadata(:).name }, keyname ) )
metadata = daffv17_add_metadata( metadata, keyname, 'String', date );
end
keyname = 'Git Version';
if ~any( strcmpi( { metadata(:).name }, keyname ) )
versionHash = ita_git_getMasterCommitHash;
metadata = daffv17_add_metadata( metadata, keyname, 'String', versionHash );
end
keyname = 'Web resource';
if ~any( strcmpi( { metadata(:).name }, keyname ) )
metadata = daffv17_add_metadata( metadata, keyname, 'String', 'http://www.ita-toolkit.org' );
end
metadata = daffv17_add_metadata( metadata, 'Generation script', 'String', 'writeDAFFFile.m' );
metadata = daffv17_add_metadata( metadata, 'Generation toolkit', 'String', 'ITA-Toolkit' );
metadata = daffv17_add_metadata( metadata, 'Generation date', 'String', date );
metadata = daffv17_add_metadata( metadata, 'Web resource', 'String', 'http://www.ita-toolkit.org' );
channels=this.nChannels/this.nDirections;
if(channels<1)
warning('Number of channels per record was not detected correctly, assuming 2 channel records');
......
function varargout = ita_HRTF_postProcessing_smoothLowFreq(varargin)
%ITA_HRTF_POSTPROCESSING_SMOOTHLOWFREQ - Fills low frequency of HRTF measurement
% This function will take a HRTF measurement and fills in low frequency content.
% Data is interpolated to 1 at 0 Hz
% Phase information is preserved.
%
% Syntax:
% audioObjOut = ita_HRTF_postProcessing_smoothLowFreq(audioObjIn, options)
%
% Options (default):
% 'cutOffFrequency' (100) : The lowest point of measured data
% 'upperFrequency' (300) : Frequency information up to this frequency is used during interpolation but not changed
% 'timeShift' (0) : Time shift the data to zero to help
% with unwrap phase
%
% Example:
% audioObjOut = ita_HRTF_postProcessing_smoothLowFreq(audioObjIn,'cutOffFrequency',100)
%
%
% Reference page in Help browser
% <a href="matlab:doc ita_HRTF_postProcessing_smoothLowFreq">doc ita_HRTF_postProcessing_smoothLowFreq</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: 21-Aug-2017
%% 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('pos1_data','itaAudio', 'cutOffFrequency', 100,'upperFrequency',300,'timeShift',0);
[input,sArgs] = ita_parse_arguments(sArgs,varargin);
%% body
% for multi instances, do a for loop
for instIndex = 1:length(input)
% first, shift the HRTF to 0. This will allow better phase interpolation
if sArgs.timeShift
shiftSamples = ita_start_IR(input(instIndex));
[tmp] = ita_time_shift(input(instIndex),-shiftSamples,'samples');
else
tmp = input(instIndex);
end
% interpolate to 0
binIdxAtLower = tmp.freq2index(sArgs.cutOffFrequency);
binIdxAtUpper = tmp.freq2index(sArgs.upperFrequency);
dataFromLowerToUpper = tmp.freqData(binIdxAtLower:binIdxAtUpper,:);
interpValues1 = interp1([0 tmp.freqVector(binIdxAtLower:binIdxAtUpper)], [ones(1,tmp.nChannels); abs(dataFromLowerToUpper)], tmp.freqVector(1:binIdxAtUpper));
interpPhase = interp1(tmp.freqVector(binIdxAtLower:binIdxAtUpper),unwrap(angle(dataFromLowerToUpper)),tmp.freqVector(1:binIdxAtUpper),'linear','extrap');
% set the interpolated values into the shifted audio
tmp.freqData(1:binIdxAtUpper,:) = interpValues1.*exp(1i.*interpPhase);
%% Add history line
tmp = ita_metainfo_add_historyline(tmp,mfilename,varargin);
if sArgs.timeShift
% shift the audio back to its original position
data_full(instIndex) = ita_time_shift(tmp,shiftSamples,'samples');
else
data_full(instIndex) = tmp;
end
end
%% Set Output
varargout(1) = {data_full};
%end function
end
......@@ -122,6 +122,7 @@ classdef itaMotorControlNanotec < itaMotorControl
pause(this.waitForSerialPort);
end
end
function moveTo(this,varargin)
% if itaCoordinates are given, the position is passed to all
......@@ -178,6 +179,13 @@ classdef itaMotorControlNanotec < itaMotorControl
end
function freeFromStopButton(this)
for index = 1:length(this.motorList)
this.motorList{index}.freeFromStopButton();
end
end
function reference(this)
this.clear_commandlist;
for index = 1:length(this.motorList)
......@@ -190,6 +198,10 @@ classdef itaMotorControlNanotec < itaMotorControl
end
this.wait4everything;
for index = 1:length(this.motorList)
this.motorList{index}.setReferenced(true);
end
end
function prepareForContinuousMeasurement(this,varargin)
......@@ -207,17 +219,21 @@ classdef itaMotorControlNanotec < itaMotorControl
% get the preangle and the speed
sArgs.preAngle = 0;
sArgs.speed = 2;
sArgs.postAngle = 10;
[sArgs notFound] = ita_parse_arguments(sArgs, varargin);
% first, do a reference move
ita_verbose_info('Moving to reference',1)
this.reference
ita_verbose_info('Disable reference position',1)
this.motorList{1}.disableReference(1);
ita_verbose_info('Moving to pre angle',1)
this.moveTo(motorName,-sArgs.preAngle,'absolut',false,'speed',1)
this.moveTo(motorName,-sArgs.preAngle,'absolut',false,'speed',2)
% now prepare the big move but don't start it
this.moveTo(motorName,360+sArgs.preAngle+15,'speed',sArgs.speed,'absolut',false,'start',0,'limit',0);
moveAngle = 360 + sArgs.preAngle + sArgs.postAngle;
this.moveTo(motorName,moveAngle,'speed',sArgs.speed,'absolut',false,'start',0,'limit',0,'direct',1);
this.preparedList = motorName;
ita_verbose_info('Finished preparing',2)
end
......@@ -246,6 +262,8 @@ classdef itaMotorControlNanotec < itaMotorControl
end
this.wait4everything
this.preparedList = [];
ita_verbose_info('Enable reference position',1)
% this.sendControlSequenceAndPrintResults(':port_in_a=7');
end
function success = add_to_commandlist(this, string_to_send)
......@@ -456,7 +474,11 @@ classdef itaMotorControlNanotec < itaMotorControl
ita_verbose_info('Position reached',1);
else
ita_verbose_info('Position NOT reached! - Check for errors!', 0);
this.send_commandlist(this.failed_command_repititions); % mpo: bugfix: send_commandlist needs argument
%assuming stop button clicked
%call freefrombutton function
% this.freeFromStopButton
% this.send_commandlist(this.failed_command_repititions); % mpo: bugfix: send_commandlist needs argument
% this.isReferenced = false;
end
this.clear_receivedlist;
......
......@@ -22,6 +22,7 @@ classdef itaMotorNanotec < itaHandle
stop(this);
active = isActive(this);
setActive(this,value);
% isInitialized(this);
getStatus(this);
getMotorID(this);
......@@ -29,7 +30,14 @@ classdef itaMotorNanotec < itaHandle
% basic moves: requires execution to halt while something is moving
this = moveToReferencePosition(this);
this = startMoveToPosition(this);
disableReference(this,value);
% this = freeFromStopButton(this);
ret = prepareMove(this,position,varargin);
end
end
methods
function setReferenced(this,value)
this.mIsReferenced = value;
end
end
end
\ No newline at end of file
......@@ -81,6 +81,10 @@ classdef itaMotorNanotec_Arm < itaMotorNanotec
name = this.motorName;
end
function disableReference(this,value)
end
function sendConfiguration(this)
motorControl = this.mMotorControl;
motorControl.add_to_commandlist(sprintf('#%d:port_in_a7\r' , this.motorID));
......
......@@ -99,6 +99,10 @@ classdef itaMotorNanotec_Slayer < itaMotorNanotec
end
function disableReference(this,value)
end
function this = moveToReferencePosition(this)
% Prepare reference move (arm)
......
......@@ -16,7 +16,7 @@ classdef itaMotorNanotec_Turntable < itaMotorNanotec
properties(Constant, Hidden = true)
sArgs_default_motor = struct( ...
'wait', true, ...
'speed', 0.5, ...
'speed', 2, ...
'VST', 'adaptiv', ...
'limit', false, ...
'continuous', false, ...
......@@ -117,8 +117,8 @@ classdef itaMotorNanotec_Turntable < itaMotorNanotec
% Set direction:
motorControl.add_to_commandlist(sprintf('#%dd=0\r' , this.motorID));
% Calculate and set lower speed:
%stepspersecond = (this.sArgs_default_motor.speed/0.9*this.sArgs_default_motor.gear_ratio);
motorControl.add_to_commandlist(sprintf('#%du=%.2f\r' , this.motorID, 25));
stepspersecond = (this.sArgs_default_motor.speed/0.9*this.sArgs_default_motor.gear_ratio);
motorControl.add_to_commandlist(sprintf('#%du=%.2f\r' , this.motorID, stepspersecond));
% Calculate and set upper speed:
stepspersecond = (this.sArgs_default_motor.speed/0.9*this.sArgs_default_motor.gear_ratio);
motorControl.add_to_commandlist(sprintf('#%do=%.2f\r' , this.motorID, stepspersecond));
......@@ -136,6 +136,9 @@ classdef itaMotorNanotec_Turntable < itaMotorNanotec
this.mMotorControl.add_to_commandlist(sprintf('#%dA\r' , this.motorID));
end
function disableReference(this,value)
end
function started = prepareMove(this,position,varargin)
......@@ -152,6 +155,7 @@ classdef itaMotorNanotec_Turntable < itaMotorNanotec
end
sArgs.continuous = false;
sArgs.direct = true;
sArgs.speed = this.sArgs_default_motor.speed;
if ~isempty(varargin)
sArgs = ita_parse_arguments(sArgs,varargin);
......@@ -170,8 +174,11 @@ classdef itaMotorNanotec_Turntable < itaMotorNanotec
end
if this.old_position.phi ~= position.phi
angle = mod(position.phi(1)/2/pi*360+360, 720)-360;
if ~sArgs.direct
angle = mod(position.phi(1)/2/pi*360+360, 720)-360;
else
angle = position.phi_deg(1);
end
ret = this.prepare_move(angle, 'absolut', true, 'wait', true, 'speed', sArgs.speed);
this.old_position = position;
started = ret;
......@@ -394,7 +401,7 @@ classdef itaMotorNanotec_Turntable < itaMotorNanotec
% Divide by 0.9 because each (half)-step is equal to 0.9 degree
% and multiply by the gear_ratio because the given angle value
% is for the turntable and not for the motor:
steps = (angle/0.9*this.sArgs_motor.gear_ratio)
steps = (angle/0.9*this.sArgs_motor.gear_ratio);
% Check if absolut or relative position mode:
if this.sArgs_motor.absolut == true
% Absolut position mode:
......
......@@ -180,6 +180,7 @@ classdef test_itaEimarMotorControl < itaMeasurementTasksScan
function moveTo(this,varargin)
% Check if it is initialized:
% varargin{1}.phi_deg
if ~this.isInitialized
ita_verbose_info('Not initialized - I will do that for you...!',0);
this.initialize
......@@ -207,18 +208,35 @@ classdef test_itaEimarMotorControl < itaMeasurementTasksScan
numRepetitions = this.measurementSetup.repetitions;
timePerRepetition = this.measurementSetup.twait*length(this.measurementSetup.outputChannels);
speed = 360/(numRepetitions*timePerRepetition);
motorControl = this.getMotorControl;
motorName = motorControl.motorList{1}.getMotorName;
if strcmp(motorName,'HRTFArc')
% preangletime
preAngle = 45;
preAngleTime = preAngle/speed; % it takes 2 seconds to start recording
% preangletime
preAngleTime = 2/64*numRepetitions; % it takes 2 seconds to start recording
preAngle = preAngleTime*speed;
postAngle = -35;
postAngleTime = postAngle/speed;
preAngle = min(preAngle,15);
preAngle = max(preAngle,8);
numTotalRepetitions = numRepetitions+ceil(preAngleTime/(timePerRepetition))+9;
this.measurementSetup.repetitions = numTotalRepetitions;
additionalReps = ceil((postAngleTime+preAngleTime + 1)/timePerRepetition);
numTotalRepetitions = numRepetitions+additionalReps;
this.measurementSetup.repetitions = numTotalRepetitions;
else
% preangletime
preAngle = 15;
preAngleTime = preAngle/speed; % it takes 2 seconds to start recording