Commit 9b96001a authored by Jan-Gerrit Richter's avatar Jan-Gerrit Richter
Browse files

Merge branch 'master' into itaMotorControl

parents 89ecc3e2 5f83ef33
function varargout = PlotComet_3D(x_data, y_data, z_data, varargin)
%function creates a moving comet plot in 3 dimensions using the specificed
%x, y, and z line data
%
%PlotComet_3D(x_data, y_data, z_data)
%
%Optional Input Arguements:
%cFigure, handle to the current figure or axes to create the plot in (will
% add the line data to the plot regardless of whether hold is on
%
%Frequency, playback frequency, if not provided assumed to be 10 Hz
%
%blockSize, number of points in the comet tail, assumed to be 1/20 of
% length(x_data) if not provided
%
%PlotComet_3D(x_data, y_data, z_data, 'cFigure',FigureHandle,...
% 'Frequency',freq, 'blockSize',num_points)
%
%tailFormat, a structured variable containing the fomrating requirements of
% the tail (if different than the default red solid line)
% ex: Tail = struct('LineWidth',2,'Color','g','LineStyle','*');
% PlotComet_3D(x_data,y_data,z_data,'tailFormat',Tail);
%
%headFormat, a structued vairable containing the formating requirements of
% the head (if different than the default blue circle)
% ex: Head = struct('LineStyle','none','MarkerSize',8,'Marker','^','Color','k');
% PlotComet_3D(x_data,y_data,z_data,'headFormat',head);
% Note: 'LineSytle','none' is required since the head is plotted as a
% single point!
%
%Optional Output Arguements
%hFigure = figure handle to the comet plot
%
%Example: Create a surface plot and construct a moving line plot through it
%with connected cyan stars and dashed line on a 50 points long tail
%and large green square for a head.
%
%%create some data for the surface
%[X, Y] = meshgrid(linspace(-1, 1, 25), linspace(-1, 1, 25));
%Z = X .* exp(-X.^2 - Y.^2);
%%open a new figure and plot the surface
%surf(X, Y, Z);
%fig = gcf;
%%create the line data
%t = -pi:pi/500:pi;
%x = sin(5*t);
%y = cos(3*t);
%z = t;
%%define the tail formated structure
%Tail = struct('LineStyle','--','Marker','*','Color','c',...
% 'LineWidth',2,'MarkerSize',4);
%%define the head formated structure
%Head = struct('LineStyle','none','Marker','s','Color','g','MarkerSize',10);
%
%%Invoke the comet plot with 50 Hz playback, and a tail length of 50 pts
%PlotComet_3D(x,y,z,'cFigure',fig,'blockSize',50,'Frequency',50,...
% 'headFormat',Head,'tailFormat',Tail);
%
%Nick Hansford
%09/26/2013
%initialize the inputs
freq = 10;
blockSize = floor(1/20*length(x_data));
tailFormat = struct('LineWidth',1,'Color','r','LineStyle','-');
headFormat = struct('LineStyle','none','Marker','o','MarkerSize',6,...
'Color','b');
pointFormat = struct('LineStyle','none','Marker','o','MarkerSize',6,...
'Color','g');
plotPoints = 0;
returnFrames = 0;
colorSwitch = 0;
moveTail = 0;
%parse out the inputs
argCount = nargin - 3;
for index = 1:2:argCount
% caseVar = varargin{i}
switch varargin{index}
case 'cFigure'
cFigure = varargin{index+1};
%get the original size:
cAxes = get(cFigure,'CurrentAxes');
xLim = get(cAxes,'XLim');
yLim = get(cAxes,'YLim');
zLim = get(cAxes,'ZLim');
%resize if the new plot will exceed them
if xLim(1) > min(x_data)
xLim(1) = min(x_data);
end
if xLim(2) < max(x_data)
xLim(2) = max(x_data);
end
if yLim(1) > min(y_data)
yLim(1) = min(y_data);
end
if yLim(2) < max(y_data)
yLim(2) = max(y_data);
end
if zLim(1) > min(z_data)
zLim(1) = min(z_data);
end
if zLim(2) < max(z_data)
zLim(2) = max(z_data);
end
axis([xLim, yLim, zLim]);
case 'blockSize'
blockSize = varargin{index+1};
case 'Frequency'
freq = varargin{index+1};
case 'tailFormat'
tailFormat = varargin{index+1};
case 'headFormat'
headFormat = varargin{index+1};
case 'plotPoints'
plotPoints = varargin{index+1};
case 'returnFrames'
returnFrames = varargin{index+1};
case 'colorSwitch'
colorSwitch = varargin{index+1};
case 'moveTail'
moveTail = varargin{index+1};
end
end
%make sure the figure exists, if not, create it
if ~exist('cFigure')
cFigure = figure;
view(3);
axis([min(x_data), max(x_data),...
min(y_data), max(y_data),...
min(z_data), max(z_data)]);
end
%user can pass in current axes, if so, get the parent for the figure window
%and use that instead...should make this compatible with GUIs?
if strcmp(get(cFigure,'Type'),'axes')
cFigure = get(cFigure,'Parent');
end
%activate the figure window
figure(cFigure);
cAxes = get(cFigure,'CurrentAxes');
oldNextPlot = get(cAxes,'NextPlot');
set(cAxes,'NextPlot','add');
pauseTime = 1./freq;
n_start = 1;
n_stop = 1;
frames = [];
deleteList = [];
%put on the starting point
deleteList{end+1} = plot3(x_data(n_start:n_stop),...
y_data(n_start:n_stop),...
z_data(n_start:n_stop),tailFormat);
deleteList{end+1} = plot3(x_data(n_stop), y_data(n_stop), z_data(n_stop),headFormat);
if returnFrames
frames{end+1} = getframe(gcf);
end
% prepare the colormap
if colorSwitch
colormap('jet');
colorValues = colormap;
colorValues = [colorValues;[0 0 0]];
indexFactor = 3.882/1.335;
else
indexFactor = 1;
end
if length(moveTail) > 1
tailX_data = squeeze(moveTail(:,1,:));
tailY_data = squeeze(moveTail(:,2,:));
tailZ_data = squeeze(moveTail(:,3,:));
end
n_blocks = length(x_data)/blockSize;
i_blocks = 0;
%playback!
for n = 1:1:ceil(length(x_data)+length(x_data)/(n_blocks*indexFactor))
a = tic;
if n <= blockSize*(i_blocks + 1)
n_start = blockSize*i_blocks + 1;
n_stop = n;
else
i_blocks = i_blocks + 1;
n_start = blockSize*i_blocks + 1;
n_stop = n;
deleteList = deleteList(2:end);
end
%delete the previous plot
for index = 1:length(deleteList)
delete(deleteList{index});
end
deleteList = [];
if n > length(x_data)
% n_start = 1;
n_stop = length(x_data);
end
%new plot
deleteList{end+1} = plot3(x_data(n_start:n_stop),...
y_data(n_start:n_stop),...
z_data(n_start:n_stop),tailFormat);
if plotPoints
for index = 1:n_stop-1
if colorSwitch
colorIndex = round((n-(index-1))*indexFactor);
if colorIndex >= length(colorValues)
colorIndex = length(colorValues);
end
pointColor = colorValues(colorIndex,:);
pointFormat = struct('LineStyle','none','Marker','o','MarkerSize',6,'Color',pointColor);
end
if moveTail == 0
deleteList{end+1} = plot3(x_data(index), y_data(index), z_data(index),pointFormat);
else
if colorIndex ~= length(colorValues)
% point
deleteList{end+1} = plot3(tailX_data(index,colorIndex), tailY_data(index,colorIndex), tailZ_data(index,colorIndex),pointFormat);
end
% tail
deleteList{end+1} = plot3(tailX_data(index,1:colorIndex),...
tailY_data(index,1:colorIndex),...
tailZ_data(index,1:colorIndex),tailFormat);
end
end
deleteList{end+1} = plot3(x_data(n_stop), y_data(n_stop), z_data(n_stop),headFormat);
else
deleteList{end+1} = plot3(x_data(n_stop), y_data(n_stop), z_data(n_stop),headFormat);
end
drawnow;
if returnFrames
frames{end+1} = getframe(gcf);
end
%update playback refresh rate
b = toc(a);
pause(pauseTime-b);
end
if plotPoints
plot3(x_data(n_stop), y_data(n_stop), z_data(n_stop),pointFormat);
drawnow;
if returnFrames
frames{end+1} = getframe(gcf);
end
end
set(cAxes,'NextPlot',oldNextPlot);
% fprintf('Finished\n');
if nargout == 1
varargout{1} = cFigure;
end
if returnFrames
varargout{2} = frames;
end
Copyright (c) 2013, Nick Hansford
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
......@@ -9,7 +9,7 @@ function writeDAFFFile( this, file_path, metadata_user )
%
% Output: none
metadata = [];
metadata = this.mMetadata;
if nargin >= 3
metadata = metadata_user;
end
......@@ -70,8 +70,11 @@ metadata = daffv17_add_metadata( metadata, 'Generation script', 'String', 'write
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 = 2; % this.nChannels < does not work?
channels=this.nChannels/this.nDirections;
if(channels<1)
warning('Number of channels per record was not detected correctly, assuming 2 channel records');
channels = 2;
end
% Content type switcher between time domain (ir) and frequency domain (dft)
% (requires different data functions)
......
function output = ita_headphone_equalization(HPTF,type)
% <ITA-Toolbox>
% This file is part of the application HRTF_Measurement for the ITA-Toolbox. All rights reserved.
% You can find the license for this m-file in the application folder.
% </ITA-Toolbox>
output = HPTF(1,1);
for cdx = 1:HPTF(1,1).nChannels
hp = merge(HPTF.ch(cdx));
switch type
case 'mean'
hp.freqData = mean(abs(hp.freqData),2);
case 'max'
hp.freqData = max(abs(hp.freqData),[],2);
case 'mSTD'
hp.freqData = mean(abs(hp.freqData),2) + 2*std(abs(hp.freqData),0,2);
otherwise
error('Unknown type');
end
%% Short Filter with no correction for low freqs and minimun phase
R = hp;
% R = ke4;
n = R.nSamples;
aux = max(abs(R.freqData),[],2);
% find first maximum and truncate low freq correction at this point
idx1 = find(R.freqVector > 100,1,'first');
idx2 = find(R.freqVector < 300,1,'last');
d_aux = diff(aux(idx1:idx2));
idx = find(diff(sign(d_aux)) ~= 0,1,'first');
aux(1:idx1+idx+1) = aux(idx1+idx+2);
aux(aux==0) = rand(1)*eps;
R.freqData = aux;
% do smoothing
%R = ita_smooth(R,'LogFreqOctave1',1/6,'Abs');
R = ita_invert_spk_regularization(R,[0 18000],'beta',.01);
% minimum phase
R = ita_time_shift(R,n/2,'samples');
R = ita_extend_dat(R,R.nSamples*2);
R = ita_time_shift(R,-n/2,'samples');
N = 2*n;
aux = R.timeData;
R.timeData = ifft(log(abs(fft(aux))));
T = 2^16;
if T >= N-3
T = N-3;
end
if rem(T-1,2) == 1
T = T+1;
end
u = [1; 2*ones((N-T-1)/2,1); cos(pi*(0:T-1)'/T)+1; zeros((N-T-1)/2,1)]; %pode colocar uma transicao mais suave no meio.
R.timeData = R.timeData.*u;
H = R;
H.timeData = ifft(exp(fft(R.timeData)),'symmetric');
H = ita_time_crop(H,[1 n],'samples');
% H = ita_time_window(H,[2^12 2^13],'samples');
% H = ita_time_window(H,[0.1 0.2],'time','dc',true);
H = ita_time_window(H,[2^9 2^10],'samples');
H.channelUnits = {''};
% It is necessary to correct the overall level of the filter
% we can either guaranty no gain
% H = H/max(abs(H.freqData));
% spl = ita_spk2level(H,3,'averaged');
%
% H = H/mean(spl.data);
output.freqData(:,cdx) = H.freqData;
end
% or guaranty that the average level of the signal is not altered
% this also means that the overall loudness of the signal will not be
% considerably altered.
% output = output/norm(output.rms);
\ No newline at end of file
......@@ -26,7 +26,7 @@ classdef itaMSTFinterleaved < itaMSTF
nWait % samples between two subsequent sweeps
end
properties(Dependent = false, Hidden = false, SetObservable = true, AbortSet = true)
repititions = 1; % how many times do you want the signal to be played?
repetitions = 1; % how many times do you want the signal to be played?
end
methods
......@@ -92,7 +92,7 @@ classdef itaMSTFinterleaved < itaMSTF
% properties.
addlistener(this,'twait','PreSet',@this.init);
addlistener(this,'outputChannels','PreSet',@this.init);
addlistener(this,'repititions','PreSet',@this.init);
addlistener(this,'repetitions','PreSet',@this.init);
end
%% INIT
......@@ -414,11 +414,11 @@ classdef itaMSTFinterleaved < itaMSTF
nSamples = result.nSamples;
% extend result
result.fftDegree = nSamples + nWaitSum * (this.repititions - 1);
result.fftDegree = nSamples + nWaitSum * (this.repetitions - 1);
timeData = single(result.time).';
idxx_init = (1:nSamples);
ita_verbose_info('itaMSTFinterleaved::appending time data.',1);
for idx = 2:this.repititions
for idx = 2:this.repetitions
idxx = idxx_init+nWaitSum * (idx-1);
timeData(:,idxx) = timeData(:,idxx) + singleTimeData;
end
......@@ -528,8 +528,8 @@ classdef itaMSTFinterleaved < itaMSTF
sArgs = struct('tIR', 0.005, 'tSpace', 0.002);
sArgs = ita_parse_arguments(sArgs, varargin);
reps = this.repititions;
this.repititions = 1;
reps = this.repetitions;
this.repetitions = 1;
% Measure Signal, Distortion and Noise:
SDN = this.run;
% Measure Distortion and Noise: (Turn off one loudspeaker each measurement)
......@@ -539,7 +539,7 @@ classdef itaMSTFinterleaved < itaMSTF
DN(idx) = temp.ch(idx); %#ok<AGROW>
end
this.pre_scaling = ones(1,numel(this.outputChannels));
this.repititions = reps;
this.repetitions = reps;
DN = merge(DN);
% Apply window:
......@@ -630,7 +630,7 @@ classdef itaMSTFinterleaved < itaMSTF
% jri: changed default output behaviour to multi instance
% this is done to have a consistent behaviour with
% each use case of the class
timeData = reshape(data.timeData(1:nSamplesWait*this.repititions*nOutputChannels, :) , nSamplesWait, data.nChannels*this.repititions*nOutputChannels);
timeData = reshape(data.timeData(1:nSamplesWait*this.repetitions*nOutputChannels, :) , nSamplesWait, data.nChannels*this.repetitions*nOutputChannels);
if (nOutputChannels > 1)
resultsMI = itaAudio(1,nOutputChannels);
for index = 1:nOutputChannels
......@@ -646,8 +646,8 @@ classdef itaMSTFinterleaved < itaMSTF
% no speed up possible... do it the slow way:
%--------------------------------------------------------------
nSamplesWait = repmat(nSamplesWait,this.repititions,1);
final_idx = numel(this.outputChannels) * this.repititions; % Compute total number of sweeps that has been played back.
nSamplesWait = repmat(nSamplesWait,this.repetitions,1);
final_idx = numel(this.outputChannels) * this.repetitions; % Compute total number of sweeps that has been played back.
% Prepare data for cropping. Extract as many samples from
% recording data as originally were in the excitation.
......@@ -674,7 +674,7 @@ classdef itaMSTFinterleaved < itaMSTF
currentData = data.';
end
for idx = ch_idx:nOutputChannels:(this.repititions * nOutputChannels)
for idx = ch_idx:nOutputChannels:(this.repetitions * nOutputChannels)
% Crop data to an interval equal to the number of the wait samples.
% Move this interval by the number of wait samples for every loop iteration.
......@@ -750,7 +750,7 @@ classdef itaMSTFinterleaved < itaMSTF
% properties to be saved during the savin process.
result = {'mTwait','mCommentData', 'repititions'};
result = {'mTwait','mCommentData', 'repetitions'};
end
......
......@@ -47,7 +47,11 @@ sArgs = struct('pos1_array','itaMicArray','pos2_f','numeric','pos3_steering_th',
%% do the calculation
% positions of array microphones
arrayPositions = array.cart;
if isempty(array.weights) || numel(array.weights) ~= array.nPoints
weights = array.w;
else
weights = array.weights;
end
% make a matrix with spherical coordinates for the unit sphere with
% given angular resolution
resolution = 1;
......@@ -67,7 +71,7 @@ k = 2*pi*f/double(ita_constants('c'));
v = squeeze(ita_beam_steeringVector(k,arrayPositions,scanPositions,sArgs.wavetype));
% ... and multiply with the manifold vector for the steering
% direction to get the beampattern
v_steer = ita_beam_steeringVector(k,arrayPositions,steer_vec,sArgs.wavetype).';
v_steer = weights(:).*ita_beam_steeringVector(k,arrayPositions,steer_vec,sArgs.wavetype).';
v = v'*v_steer./sum(abs(v_steer).^2);
B = reshape(v,numel(theta),numel(phi));
......@@ -75,7 +79,7 @@ B = B./max(abs(B(:))); % normalize to maximum
if ~strcmpi(sArgs.plotPlane,'none')
v = itaResult(v(:).'./max(abs(v(:))),f,'freq');
v.channelCoordinates = itaCoordinates(scanPositions.');
v.channelCoordinates = itaCoordinates(scanPositions);
if (strcmpi(sArgs.plotPlane,'xyz') || strcmpi(sArgs.plotPlane,'3d'))
switch sArgs.plotType
case 'lin'
......
function [ loudspeakerSignals ] = ita_decodeAmbisonics( Bformat, LoudspeakerPos, varargin )
%ITA_DECODEAMBISONICS Summary of this function goes here
% Detailed explanation goes here
opts.decoding='remax' % Decoding strategy (remax,inphase,plane)
opts = ita_parse_arguments(opts,varargin);
% Initializing further parameters
nmax=numel(Bformat);
% Weighting for Inphase/ReMax
g_re=ones(nmax,1); %init weighting factors for ReMax
g_in=g_re; %init weighting factors for InPhase
% ReMax Decoding
if(sum(strcmp(lower(opts.decoding),{'remax' 'both'})))
% Calculates the B-Format channel weights for ReMax-Decoding
% see J.Daniel Thesis p.312 A.67/A.68
weightsReMax=zeros(nmax,1); % init weights
% g_m=P_m(r_E) with n=0 for P
% r_E is largest root of P_(order+1)
syms x;
f=1/2^(TheOrder+1)/factorial(TheOrder+1)*diff(((x^2-1)^(TheOrder+1)),(TheOrder+1));%Legendre Polynom n=0 m=order+1
maxroot=max(eval(solve(f))); %find maximum root(Nullstelle)
for k=1:nmax
leggie=legendre(ita_sph_linear2degreeorder(k),abs(maxroot)); % g_m=P_m(r_E)
weightsReMax(k)=leggie(1);% pick n=0
end
g_re=weightsReMax;
end
% Inphase Decoding
if(sum(strcmp(lower(opts.decoding),{'inphase' 'both'})))
% Calculates the B-format channel weights for InPhase-Decoding
TheOrder = obj.ambGetOrder; % amb order from sim room
nmax = (TheOrder+1)^2; % number of channels
weightsInPhase=zeros(nmax,1);
N=numel(obj.monitorRoom.getSourcePosition)/3;% get number of loudspeaker