Skip to content

GitLab

toolbox
toolbox
Commit 277a6781 authored by Philipp Marco Schäfer's avatar Philipp Marco Schäfer
Browse files
Options

removed application/Numeric/SH 

-seem to be unused files originating from Kunkemöller's Diplomarbeit
parent ba929e71
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 3254 deletions
applications/Numeric/SH/@itaAudioMIMO/itaAudioMIMO.m deleted 100644 → 0
View file @ ba929e71
classdef itaAudioMIMO < itaAudio
% <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>
properties(Access = private, Hidden = true)
mReceiverPosition = itaCoordinates();
mReceiverUpVector = itaCoordinates();
mReceiverViewVector = itaCoordinates();
mSourcePosition = itaCoordinates();
mSourceUpVector = itaCoordinates();
mSourceViewVector = itaCoordinates();
end
properties(Dependent = true, Hidden = false)
nSources
nReceivers
end
properties(Dependent = true, Hidden = true)
end
methods
function this = itaAudioMIMO(varargin)
this = this@itaAudio(varargin{1}(1));
[nReceivers, nSources] = size(varargin{1});
nSamplesOrBins = size(this.data, 1);
this.(this.domain) = zeros(nSamplesOrBins, nReceivers, nSources);
channelNames = cell(nReceivers, nSources);
for iReceiver = 1:nReceivers
for iSource = 1:nSources
if varargin{1}(iReceiver, iSource).nChannels > 1
error('nur ein channel pro itaAudio. muss man anpassen wenn man was anderes will')
end
this.(this.domain)(:, iReceiver, iSource) = varargin{1}(iReceiver, iSource).(this.domain);
channelNames(iReceiver, iSource) = varargin{1}(iReceiver, iSource).channelNames;
end
end
this.channelNames = channelNames;
end
function res = source(this, iSource)
res = this;
dataSize = [this.nReceivers this.nSources];
% use linear index because channelNames is always 1-D and data 2-D
linearIndex = sub2ind(dataSize, 1:dataSize(1), iSource*ones(1,dataSize(1)));
res.data = this.data(:, linearIndex);
res.channelNames = this.channelNames(linearIndex);
end
function res = receiver(this, iReceiver)
res = this;
dataSize = [this.nReceivers this.nSources];
% use linear index because channelNames is always 1-D and data 2-D
linearIndex = sub2ind(dataSize, iReceiver*ones(1,dataSize(2)), 1:dataSize(2));
res.data = this.data(:, linearIndex);
res.channelNames = this.channelNames(linearIndex);
end
function nSources = get.nSources(this)
nSources = size(this.(this.domain), 3);
end
function nReceivers = get.nReceivers(this)
nReceivers = size(this.(this.domain), 2);
end
end
end
\ No newline at end of file
applications/Numeric/SH/@itaSphSynthDirectivity/SHfilter2SHrir.m deleted 100644 → 0
View file @ ba929e71
function SHfilter2SHrir(this)
% This function extends the frequency range of the filters, calculated by
% itaSphSynthDirectivity.makeSHfilter, and convolve them with the
% measurement.
%
% The result is first saved in a big matrix 'rirData' and then exportet in
% itaAudio-format into the directory this.folder/SH_RIR
%
% see also: itaSphSynthDirectivity, makeSynthArray, makeSHfilter, convolve_itaBalloon_and_SHfilter
if ~strcmpi(this.outputFormat, 'itaAudio')
error(' ');
end
%check measurement data files
for idxD = 1:length(this.measurementDataFolder)
testFile = [this.measurementDataFolder{idxD} filesep this.filemask '1.ita'];
if ~exist(testFile,'file');
error(['There is no such file: "' testFile '.ita"']);
end
end
nmax_l = (this.nmax+1)^2;
actFF = [0 0]; %index Folder and File of the current measurementdata in the memory
freqOffset = length(this.array.freqVector(this.array.freqVector < this.freqVector(1)));
this.myWaitbar(this.nSpeakers+1, 'SHfilter2SHrir');
for idxS = 1:this.nSpeakers
this.myWaitbar(idxS);
% parse filter data 2 itaAudio
filter = itaAudio('samplingRate',this.array.samplingRate, 'signalType','energy',...
'dataType', this.precision);
filter.freqData = zeros(this.array.nBins, nmax_l, this.precision);
filter.freqData(freqOffset+(1:this.nBins),:) = ...
permute(this.mFilterData.get_data(1:this.nBins, idxS, 1:nmax_l), [1 3 2]);
filter = this.extendFreqRange(filter);
%read measurement data
TRC = this.speaker2idxTiltRotCh(idxS,:); %tilt & rotation angle, channel
FF = this.idxTiltRot2idxFolderFile{TRC(1),TRC(2)}; %folder file
if sum(actFF == FF) ~= 2
% only load if neccessary
data = ita_read([this.measurementDataFolder{FF(1)} filesep this.filemask int2str(FF(2)) '.ita']);
actFF = FF;
end
% initialize output's freqData
if idxS == 1
rirData = zeros(data(1).nBins, nmax_l, data(1).nChannels, this.precision);
nMic = data(1).nChannels;
end
%adapt length of filter and data
if data(1).nSamples ~= filter.nSamples
if data(1).nSamples < filter.nSamples
% that should never happen
error('sorry, I did not expect that, maybe you could code that?');
else
filter = filter.';
filter = ita_time_window(filter, round(filter.nSamples/2+[-0.005*filter.samplingRate 0]),'samples','symmetric');
filter = ita_extend_dat(filter, data(1).nSamples,'symmetric');
filter = filter';
end
end
% convolve
for idxM = 1:nMic
rirData(:,:,idxM) = rirData(:,:,idxM) + bsxfun(@times, cast(data(TRC(3)).freqData(:,idxM), this.precision), filter.freqData);
end
clear filter;
end
this.myWaitbar(this.nSpeakers+1);
% initialize output objects
if ~isdir([this.folder filesep 'SH_RIR']),
mkdir([this.folder filesep 'SH_RIR']);
end
for idxD = 1:nMic
out = itaAudio;
out.samplingRate = this.array.samplingRate;
out.signalType = 'energy';
out.dataType = this.precision;
out.freqData = rirData(:,:,idxD);
out.comment = ['SH - RIR (' data(1).channelNames{idxD} ')'];
out.userData = struct('freqRange', this.freqRange);
for idxC = 1:out.nChannels
[n m] = ita_sph_linear2degreeorder(idxC);
out.channelNames{idxC} = ['sph ' int2str(n) ', ' int2str(m)];
end
ita_write(out, [this.folder filesep 'SH_RIR' filesep 'SH_RIR_Mic' int2str(idxD) '.ita']);
end
this.myWaitbar([]);
applications/Numeric/SH/@itaSphSynthDirectivity/convolve_itaBalloon_and_SH_RIR.m deleted 100644 → 0
View file @ ba929e71
function out_m = convolve_itaBalloon_and_SH_RIR(this, balloon, varargin)
%
% Convolves a target-directivity (an itaBalloon-object) and the impulse
% responses of abstract directivities with the form of spherical harmonic
% basefunctions.
% see also: makeSHfilter and SHfilter2SHrir
%
% input:
% - balloon : target - directivity - function
%
% options:
% 'channels' you can choose one ore multiple channels (directivity of a
% multichannel itaBalloon will be sumed up).
% 'mpb_filter' result will be band widhth filtered by ita_mpb_filter
% 'rotate' here you can give a set of euler rotation angles to rotate the input balloon.
% The output will be an array of filters- one for each
% position
% 'rotate', {[orientation 1], [orientation 2], ... }
sArgs = struct('channels',1:balloon.nChannels,'rotate',zeros(1,3));
if nargin > 2
sArgs = ita_parse_arguments(sArgs, varargin);
end
if ~isdir([this.folder filesep 'SH_RIR']) || ~numel(dir([this.folder filesep 'SH_RIR' filesep '*.ita']))
error('First proceed itaSphSyntDirectivity.SHfilter2SHrir');
end
if ~iscell(sArgs.rotate), sArgs.rotate = {sArgs.rotate}; end
for idxR = 1:length(sArgs.rotate)
if size(sArgs.rotate{idxR},2)~=3, error('size(rotatate,2) != 3 (euler angle)'); end
end
% target's directivity
ao = balloon.idxCoefSH2itaAudio(1:this.mFilterData.dimension(3),'channels',sArgs.channels,'sum_channels', true);
ao.dataType = this.precision;
if balloon.fftDegree ~= this.array.fftDegree
error('I can not handle different fftDegrees. Please use itaBalloon.convert_fftDegree as work arround.');
end
% synthesised abstract room impulse responses
nMic = numel(dir([this.folder filesep 'SH_RIR' filesep 'SH_RIR_Mic*.ita']));
out = itaAudio(length(sArgs.rotate),nMic);
for idxM = 1:nMic
RIR = ita_read([this.folder filesep 'SH_RIR' filesep 'SH_RIR_Mic' int2str(idxM) '.ita']);
for idxR = 1:length(sArgs.rotate)
out(idxR, idxM) = itaAudio('dataType', this.precision, 'signalType', 'energy','samplingRate', this.array.samplingRate);
ao2 = ao;
if strcmpi(this.SHType, 'complex')
ao2.freqData = ita_sph_rotate_complex_valued_spherical_harmonics(ao.freqData.',sArgs.rotate{idxR}).';
elseif strcmpi(this.SHType, 'real')
ao2.freqData = ita_sph_rotate_real_valued_spherical_harmonics(ao.freqData.',sArgs.rotate{idxR}).';
else
error('unknown SHType');
end
% adapt data
if ao2.nSamples < RIR.nSamples
ao2 = ita_time_window(ao2, round(ao2.nSamples/2+[-0.005*ao2.samplingRate 0]),'samples','symmetric');
ao2 = ita_extend_dat(ao2, RIR.nSamples,'symmetric');
else
ao2 = ita_extract_dat(ao2, RIR.nSamples,'symmetric');
ao2 = ita_time_window(ao2, round(ao2.nSamples/2+[-0.005*ao2.samplingRate 0]),'samples','symmetric');
end
%convolve and add
out(idxR, idxM) = sum(ita_multiply_spk(ao2,RIR));
%adapt latency samples
if balloon.latencySamples ~= this.array.latencySamples
out(idxR, idxM) = ita_time_shift(out(idxR, idxM), balloon.latencySamples - this.array.latencySamples, 'samples');
end
out(idxR, idxM).channelNames{1} = '';
out(idxR, idxM).comment = ['synthesized RIR of ' balloon.name ' (mic: ' int2str(idxM) ')'];
end
end
out_m = itaAudio(length(sArgs.rotate),1);
for idxR = 1:length(sArgs.rotate)
out_m(idxR) = merge(out(idxR,:));
out_m(idxR).history = {' '};
end
\ No newline at end of file
applications/Numeric/SH/@itaSphSynthDirectivity/evaluate_synthesisError.m deleted 100644 → 0
View file @ ba929e71
function varargout = evaluate_synthesisError(this, varargin)
% returns and plots the relative root mean square error of the synthesis of
% abstract directivities with the surface of spherical harmonics.
%
% options: freq: vector of frequencies beeing evaluated
% degault : this.freqRange in 1/3-octave steps
% namx: maximum degree of evaluation
%
sArgs = struct('freq', [], 'nmax', this.nmax);
if nargin > 1
sArgs = ita_parse_arguments(sArgs, varargin);
end
if isempty(sArgs.freq)
sArgs.freq = ita_ANSI_center_frequencies(this.freqRange, 3);
end
RMS = zeros((sArgs.nmax+1)^2, length(sArgs.freq), this.precision);
filter = this.mFilterData.get_data(this.freq2idxFreq(sArgs.freq), 1:this.nSpeakers, 1:(sArgs.nmax+1)^2);
for idxF = 1:length(sArgs.freq)
RMS(:,idxF) = sum(abs(eye((sArgs.nmax+1)^2) - this.freq2coefSH_synthArray(sArgs.freq(idxF), 'nmax', sArgs.nmax)...
* permute(filter(idxF,:,:), [2 3 1])).^2,2);
end
RMS = sqrt(RMS);
ita_sph_plot_coefs_over_freq(RMS, sArgs.freq, 'type','max');
if nargout
varargout = {RMS};
else
varargout = {};
end
applications/Numeric/SH/@itaSphSynthDirectivity/extendFreqRange.m deleted 100644 → 0
View file @ ba929e71
function ao = extendFreqRange(this, ao)
%
% this function smoothes and extends the a calculated frequency response
% using rational fit (must be improoved) or polynomes
%
overlap = 1/3; % overlap of sections
fitSecDegree = [3 18];
%% initialize
freqData = ao.freqData;
outData = zeros(size(freqData), this.precision);
% estimate frequency range
channel = 1;
while ~sum(abs(freqData(:,channel))) && channel < ao.nChannels
channel = channel + 1;
end
idx_FR = [find(abs(freqData(:,channel))>0, 1, 'first'),...
find(abs(freqData(:,channel))>0, 1, 'last')];
% mean amplitude
meanAmp = mean(abs(freqData(idx_FR(1):idx_FR(2),:)), 1);
if idx_FR(1) > 1
% indicees to extend the first part of the frequency range ...
ids_fit = ao.freq2index(ao.freqVector(idx_FR(1)*[1 2^overlap]));
nAddPs = floor(min(length(ids_fit), ids_fit(1)/2)/2)*2;
ids_all = [(1:nAddPs).'; ids_fit];
newIds = ids_all(1):ids_all(end);
extend_low = true;
fitSecDegree(1) = min(fitSecDegree(1), length(ids_all)-1);
else
extend_low = false;
ids_fit = [1 1];
end
if idx_FR(2) < ao.nBins
% indicees to extend the end of the frequency range
ide_fit = ao.freq2index(ao.freqVector(idx_FR(2)*[2^(-overlap) 1]));
nAddPe = length(ide_fit);
p1 = 0.25*ide_fit(end) + 0.75*ao.nBins; pe = ao.nBins;
addPe = (round(p1) : round((pe-p1)/(nAddPe)) : round(pe)).';
addPe(end) = ao.nBins;
ide_all = [ide_fit; addPe];
newIde = ide_all(1):ide_all(end);
extend_high = true;
fitSecDegree(2) = min(fitSecDegree(2), length(ide_all)-5);
else
extend_high = false;
ide_fit = [1 1]*ao.nBins;
end
%% xover
xover = zeros(size(freqData,1),2, this.precision);
xover(1:ids_fit(1) ,1) = 1;
xover(ids_fit(2):ide_fit(1),2) = 1;
xover(ide_fit(2):end, 3) = 1;
if extend_low
% first X
lX = ids_fit(2)-ids_fit(1);
xover(ids_fit(1):ids_fit(2),1) = 1 - (0:lX)/lX;
xover(ids_fit(1):ids_fit(2),2) = (0:lX)/lX;
end
if extend_high
% second X
lX = ide_fit(2)-ide_fit(1);
xover(ide_fit(1):ide_fit(2),2) = 1 - (0:lX)/lX;
xover(ide_fit(1):ide_fit(2),3) = (0:lX)/lX;
end
nBins = ao.nBins;
for idxC = 1:ao.nChannels
newStuff = zeros(nBins,2, this.precision);
if extend_low
% extend the first part of the frequency range ...
datas = freqData(ids_fit, idxC);
[P S muX] = polyfit(ids_all, [repmat(meanAmp(idxC), nAddPs, 1); datas] , fitSecDegree(1)); %#ok<ASGLU>
newStuff(newIds,1) = polyval(P, (newIds-muX(1))/muX(2));
end
if extend_high
% extend the second part of the frequency range ...
datae = freqData(ide_fit, idxC);
[P S muX] = polyfit(ide_all, [datae; repmat(meanAmp(idxC), length(addPe), 1)], fitSecDegree(2)); %#ok<ASGLU>
newStuff(newIde,2) = polyval(P, (newIde-muX(1))/muX(2));
end
%% xover ber Abschnitte
outData(:,idxC) = sum([newStuff(:,1) freqData(:,idxC) newStuff(:,2)] .* xover, 2);
end
ao.freqData = outData;
applications/Numeric/SH/@itaSphSynthDirectivity/getPositions.m deleted 100644 → 0
View file @ ba929e71
function this = getPositions(this, varargin)
% reads rotationAngle out of a set of measurement data set by
% measurementDataFolder = {'folder1', 'folder2',...}
%
% you can use 'ita_roomacoustics_correlation_coefficient_longtime to select good
% speaker positions
%
% settings:
% 'nPos' : maximum number of rotations per folder to build a synthSuperSpeaker
% 'corcoef' : the result of
% "ita_roomacoustics_correlation_coefficient_longtime", that
% determines the best correlating measurement positions
sArgs = struct('corcoef',[], 'nPos',[]);
%% init
check_empty(this)
if nargin > 1
sArgs = ita_parse_arguments(sArgs,varargin);
end
if isempty(this.measurementDataFolder)
error('give me some "measurementDataFolder"');
end
if ~isempty(sArgs.corcoef)
% select speaker positions via correlation analyis
[dummy idxPositions] = sort(sArgs.corcoef.timeData,'descend'); %#ok<ASGLU>
if ~isempty(sArgs.nPos)
idxPositions = idxPositions(1:sArgs.nPos);
end
else
idxPositions = []; lastIndex = 0;
for idxF = 1:length(this.measurementDataFolder)
nFiles = numel(dir([this.measurementDataFolder{idxF} filesep this.filemask '*.ita']));
if ~nFiles, error(['can not find no file "' this.measurementDataFolder{idxF} filesep this.filemask '*.ita"']); end
% select speaker positions via a maximum number of positions
if ~isempty(sArgs.nPos)
idxPositions = [idxPositions, lastIndex + (1:floor(nFiles/sArgs.nPos):nFiles)]; %#ok<AGROW>
if length(idxPositions) > sArgs.nPos
idxPositions = idxPositions(1:sArgs.nPos);
end
if length(idxPositions) ~= sArgs.nPos
disp('check');
end
% take all measured data
else
idxPositions = [idxPositions, lastIndex + (1:nFiles)]; %#ok<AGROW>
end
lastIndex = lastIndex + nFiles;
end
end
% routing : position 2 folder/file
idxPos2idxFolderFile = zeros(0,2);
for idxF = 1:length(this.measurementDataFolder)
nFiles = numel(dir([this.measurementDataFolder{idxF} filesep this.filemask '*.ita']));
if ~nFiles, error(['can not find no file "' this.measurementDataFolder{idxF} filesep this.filemask '*.ita"']); end
idxPos2idxFolderFile = [idxPos2idxFolderFile; ones(nFiles,1)*idxF (1:nFiles).']; %#ok<AGROW>
end
idxPos2idxFolderFile = idxPos2idxFolderFile(idxPositions,:);
%% read rotation angle out of measurement data
% (*) Im Moment gilt noch: idxTilt = idxFolder.
% Spter sollte / knnte man auch idxTilt mal aus den Messdaten auslesen (?)
nTilt = length(this.measurementDataFolder);
nRotMax = 0;
for idxT = 1:nTilt
nRotMax = max(nRotMax, length(idxPos2idxFolderFile(idxPos2idxFolderFile(:,1) == idxT, 2)));
end
this.idxTiltRot2idxFolderFile = cell(nTilt, nRotMax);
this.rotationAngle = cell(nTilt,1);
%get rotation angles from channel data.coordinates
for idxP = 1:length(idxPositions)
idxFolder = idxPos2idxFolderFile(idxP,1);
idxFile = idxPos2idxFolderFile(idxP,2);
data = ita_read([this.measurementDataFolder{idxFolder} filesep this.filemask int2str(idxFile) '.ita']);
if ~this.rotationAngle_counterClockWise
phi = mod(10*pi - data(1).channelCoordinates.phi(1), 2*pi);
else
phi = mod(data(1).channelCoordinates.phi(1), 2*pi);
end
% (*)
idxTilt = idxFolder;
this.rotationAngle{idxTilt} = [this.rotationAngle{idxTilt} phi];
this.idxTiltRot2idxFolderFile{idxTilt, length(this.rotationAngle{idxTilt})} = [idxFolder idxFile];
end
\ No newline at end of file
applications/Numeric/SH/@itaSphSynthDirectivity/ita_test_itaSphSynthDirectivity.m deleted 100644 → 0
View file @ ba929e71
%% Short demo for itaSphSynthDirectivity class
%
%
% Das hier dient vorlufig auch als Demo fr diese Klasse. Ist im Prinzip
% bisschen wir itaBalloon: Einen Haufen settings geeignet setzen, play
% drcken, warten, warten und dann ein paar Auswertungsroutinen genieen
%
%
%
%% settings
this = itaSphSynthDirectivity;
this.folder = [homeFolder '\evaluation\superDode'];
this.name = 'superDode';
this.array = ita_read_itaBalloon([balloonHomeFolder '\Dode_Mid2_c\DODE_MID']);
% basic settings
this.freqRange = [100 16000];
this.arrayNmax = 30; % maximale Ordnung, bis zu der das synthetisch erweiterete Array berechnet wird
this.arrayChannels = 1:12; % die Channels des itaBalloons, die in der Messung angewendet wurden
this.nmax = 30; % maximale Ordnung, bis zu der alles andere Berechnet wird
this.precision = 'single';
% important stuff:
this.measurementDataFolder = ... % Pade der Messdaten
{[homeFolder filesep measurementFolder{1} '\data_p'], ...
[homeFolder filesep measurementFolder{2} '\data_p'], ...
[homeFolder filesep measurementFolder{3} '\data_p']};
this.tiltAngle = ... % zugeordnete Euler - Kippwinkel
{[0 180 0; 38.5+80 16.2 0]* pi/180, ...
[0 180 0; 38.5+0 16.2 0] * pi/180, ...
[0 180 0; 38.5+40 16.2 0]* pi/180};
this.filemask = 'dode'; % filemask
this.rotationAngle_counterClockWise = false; % angles by itaItalian : dreht im Urzeigersinn !!
this.getPositions('nPos',20); % liest this.rotationAngle aus Messdaten, setzt rooting channel2
save(this);
%% --- Der Abschnitt dauert bisschen (trotz intensiven Gebrauchs des profilers...)!!
% engster Flaschenhals : fft in 'SHfilter2SHrir'
%
% Fr diese Funktionen gibt es jeweils aufschlussreiche Dokumentation
%
disp('makin awesome enlarged array');
this.makeSynthArray; % die SH-Koeffizienten des synthetischen Arrays werden berechnet (Die groe D-Matrix)
disp('makin awesome synthesis filter');
this.makeSHfilter; % Invertierung der D-Matrix, regularisiert. + Umsortieren der Vorzugsrichtung:
% Bislang waren die Daten nach Frequenzen
% unterteilt, jetzt nach Einzellautsprechern
disp('makin awesome synthesized room impulse responses');
this.SHfilter2SHrir; % Die Filter werden mit da Messung gefaltet
% Quelle synthetisieren
target = ita_read_itaBalloon([balloonHomeFolder '\CUBEd\CUBE']);
ao = this.convolve_itaBalloon_and_SH_RIR(target);
%% --- Auswertung -----------------
% coeffizienten des erweiterten Arrays:
coef = this.freq2coefSH_synthArray(2000); % -> see doc
% 'mein' graph
RMS = this.evaluate_synthesisError; % -> see doc
% filter anschauen (index von zu synthetisierender sphericsal harmonic 2
% filter)
ao = this.idxSH2filter(1);
%% Vergleich mit Messung;
target = ita_read_itaBalloon([balloonHomeFolder '\CUBEd\CUBE']);
orgFiles = [homeFolder filesep measurementFolder{7} '\data_p\cube'];
org = itaAudio(numel(dir([orgFiles '*'])),1);