Commit 32d8c056 authored by Joanna Luberadzka's avatar Joanna Luberadzka

add first content

parents
function [probability DOAazimuth time] = DOASVMloc(sig,fs)
% Inputs:
% sig: 4ch x nsampl input signal
% channel 1: front left BTE
% channel 2: front right BTE
% channel 3: middle left BTE
% channel 4: middle right BTE
% fs: fs of input signal
%
% Outputs:
% probability: 72 x T DOA probability map
% DOAazimuth: directions in , corrsponding to the first dimension of probability
% time: time in s
%
% HK 21-01-2016
% hendrik.kayser@uni-oldenburg.de
% load DOASVM model
sSensors = 'frontmiddle';
sModeltype = 'HRIR';
LocModel = load(['LocModel_' sModeltype '_' sSensors]);
DOAazimuth = LocModel.az;
time = 0:1./LocModel.fs.*LocModel.wndlen./2: size(sig,2)./LocModel.fs -1./LocModel.fs.*LocModel.wndlen;
% rearrange LocModel parameters
for IDXpos = 1: size(LocModel.oModel,2)
LocModel.w(IDXpos,:) = LocModel.oModel(IDXpos).stModel.w(1:end-1);
LocModel.b(IDXpos) = LocModel.oModel(IDXpos).stModel.w(end);
LocModel.x(IDXpos) = LocModel.oModel(IDXpos).data(1);
LocModel.y(IDXpos) = LocModel.oModel(IDXpos).data(2);
end
% STFT from Multi-Channel BSS Locate toolbox R. Lebarbenchon and E. Camberlein)
% http://bass-db.gforge.inria.fr/bss_locate/
% Resample to model fs
if fs ~= LocModel.fs
sig = resample(sig.',LocModel.fs,fs).';
end
X = MBSS_stft_multi(sig,LocModel.wndlen);
% Feature extraction
mFeature = [];
for IDXpair = 1:size(LocModel.vPairs,2)
G_length = LocModel.wndlen./2+1;
ifftwin = hanning_(2.*G_length);
ifftwin = ifftwin(G_length+1:end);
G = squeeze(conj(X(:,:,LocModel.vPairs(1,IDXpair))) .* X(:,:,LocModel.vPairs(2,IDXpair)));
% PHAT weighting
G = G./max(abs(G),eps);
% Apply ifft window
G = G .* ifftwin(:,ones(1,size(G,2))) .*LocModel.nUpsample;
% zero-padding according to upsampling factor
G = [G;zeros(floor((LocModel.nUpsample-1).*LocModel.wndlen./2),size(G,2))];
% reconstruct mirror frequencies, assumes even wndlen
G = [G; conj(G(end-1:-1:2,:))];
% ifft and swap halves of csd to compensate for fftshift
vGCC = ifftshift(ifft(G,[], 1,'symmetric'),1);
% Use only relevant part of GCC-PHAT function
vDelaybins = size(vGCC,1)./2-LocModel.maxLag(IDXpair)+1:size(vGCC,1)./2+LocModel.maxLag(IDXpair)+1;
mFeature = [mFeature; vGCC(vDelaybins,:)];
end
% Compute DOA probability map
for IDXt = 1 : size(mFeature,2)
% Apply linear SVM model (one for each direction) to feature vector
c = (LocModel.w * mFeature(:,IDXt)) + LocModel.b.';
% map to probability using a sigmoid transformation
probability(:,IDXt) = 1 ./ (1 + exp(-(LocModel.x.' + c .* LocModel.y.')));
end
function [probability DOAazimuth time] = DOASVMloc1(sig,fs)
% Inputs:
% sig: 4ch x nsampl input signal
% channel 1: front left BTE
% channel 2: front right BTE
% channel 3: middle left BTE
% channel 4: middle right BTE
% fs: fs of input signal
%
% Outputs:
% probability: 72 x T DOA probability map
% DOAazimuth: directions in �, corrsponding to the first dimension of probability
% time: time in s
%
% HK 21-01-2016
% hendrik.kayser@uni-oldenburg.de
% load DOASVM model
sSensors = 'frontmiddle';
sModeltype = 'HRIR';
LocModel = load(['LocModel_' sModeltype '_' sSensors]);
DOAazimuth = LocModel.az;
time = 0:1./LocModel.fs.*LocModel.wndlen./2: size(sig,2)./LocModel.fs -1./LocModel.fs.*LocModel.wndlen;
% STFT from Multi-Channel BSS Locate toolbox R. Lebarbenchon and E. Camberlein)
% http://bass-db.gforge.inria.fr/bss_locate/
% Resample to model fs
if fs ~= LocModel.fs
sig = resample(sig.',LocModel.fs,fs).';
end
X = MBSS_stft_multi(sig,LocModel.wndlen);
% Feature extraction
mFeature = [];
for IDXpair = 1:size(LocModel.vPairs,2)
G_length = LocModel.wndlen./2+1;
ifftwin = hanning_(2.*G_length);
ifftwin = ifftwin(G_length+1:end);
G = squeeze(conj(X(:,:,LocModel.vPairs(1,IDXpair))) .* X(:,:,LocModel.vPairs(2,IDXpair)));
% PHAT weighting
G = G./max(abs(G),eps);
% Apply ifft window
G = G .* ifftwin(:,ones(1,size(G,2))) .*LocModel.nUpsample;
% zero-padding according to upsampling factor
G = [G;zeros(floor((LocModel.nUpsample-1).*LocModel.wndlen./2),size(G,2))];
% reconstruct mirror frequencies, assumes even wndlen
G = [G; conj(G(end-1:-1:2,:))];
% ifft and swap halves of csd to compensate for fftshift
vGCC = ifftshift(ifft(G,[], 1,'symmetric'),1);
% Use only relevant part of GCC-PHAT function
vDelaybins = size(vGCC,1)./2-LocModel.maxLag(IDXpair)+1:size(vGCC,1)./2+LocModel.maxLag(IDXpair)+1;
mFeature = [mFeature; vGCC(vDelaybins,:)];
end
% Compute DOA probability map
for IDXt = 1 : size(mFeature,2)
% Apply linear SVM model (one for each direction) to feature vector
c = (LocModel.w * mFeature(:,IDXt)) + LocModel.b.';
% map to probability using a sigmoid transformation
probability(:,IDXt) = 1 ./ (1 + exp(-(LocModel.x.' + c .* LocModel.y.')));
end
function [probability DOAazimuth time] = DOASVMloc(sig,fs)
% Inputs:
% sig: 4ch x nsampl input signal
% channel 1: front left BTE
% channel 2: front right BTE
% channel 3: middle left BTE
% channel 4: middle right BTE
% fs: fs of input signal
%
% Outputs:
% probability: 72 x T DOA probability map
% DOAazimuth: directions in �, corrsponding to the first dimension of probability
% time: time in s
%
% HK 21-01-2016
% hendrik.kayser@uni-oldenburg.de
% load DOASVM model
sSensors = 'frontmiddle';
sModeltype = 'HRIR';
LocModel = load(['LocModel_' sModeltype '_' sSensors]);
DOAazimuth = LocModel.az;
time = 0:1./LocModel.fs.*LocModel.wndlen./2: size(sig,2)./LocModel.fs -1./LocModel.fs.*LocModel.wndlen;
% STFT from Multi-Channel BSS Locate toolbox R. Lebarbenchon and E. Camberlein)
% http://bass-db.gforge.inria.fr/bss_locate/
% Resample to model fs
if fs ~= LocModel.fs
sig = resample(sig.',LocModel.fs,fs).';
end
X = MBSS_stft_multi(sig,LocModel.wndlen);
% Feature extraction
mFeature = [];
for IDXpair = 1:size(LocModel.vPairs,2)
G_length = LocModel.wndlen./2+1;
ifftwin = hanning_(2.*G_length);
ifftwin = ifftwin(G_length+1:end);
G = squeeze(conj(X(:,:,LocModel.vPairs(1,IDXpair))) .* X(:,:,LocModel.vPairs(2,IDXpair)));
% PHAT weighting
G = G./max(abs(G),eps);
% Apply ifft window
G = G .* ifftwin(:,ones(1,size(G,2))) .*LocModel.nUpsample;
% zero-padding according to upsampling factor
G = [G;zeros(floor((LocModel.nUpsample-1).*LocModel.wndlen./2),size(G,2))];
% reconstruct mirror frequencies, assumes even wndlen
G = [G; conj(G(end-1:-1:2,:))];
% ifft and swap halves of csd to compensate for fftshift
vGCC = ifftshift(ifft(G,[], 1,'symmetric'),1);
% Use only relevant part of GCC-PHAT function
vDelaybins = size(vGCC,1)./2-LocModel.maxLag(IDXpair)+1:size(vGCC,1)./2+LocModel.maxLag(IDXpair)+1;
mFeature = [mFeature; vGCC(vDelaybins,:)];
end
% Compute DOA probability map
for IDXt = 1 : size(mFeature,2)
% Apply linear SVM model (one for each direction) to feature vector
c = (LocModel.w * mFeature(:,IDXt)) + LocModel.b.';
% map to probability using a sigmoid transformation
probability(:,IDXt) = 1 ./ (1 + exp(-(LocModel.x.' + c .* LocModel.y.')));
end
function X=MBSS_stft_multi(x,wlen)
% fichier : MBSS_stft_multi.m
% Multichannel short-time Fourier transform (STFT) using
% half-overlapping sine windows.
%
% X=MBSS_stft_multi(x,wlen)
%
% Inputs:
% x: nchan x nsampl matrix containing nchan time-domain mixture signals
% with nsampl samples
% wlen: window length (default: 1024 samples or 64ms at 16 kHz, which is
% optimal for speech source separation via binary time-frequency masking)
%
% Output:
% X: nbin x nfram x nchan matrix containing the STFT coefficients with nbin
% frequency bins and nfram time frames
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2008 Emmanuel Vincent
% Copyright 2015 Ewen Camberlein and Romain Lebarbenchon
% This software is distributed under the terms of the GNU Public License
% version 3 (http://www.gnu.org/licenses/gpl.txt)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Difference with 2008 copyrigth:
% - signal is truncated instead of using zero padding to respect window
% length used (wlen)
%%% Errors and warnings %%%
if nargin<1, error('Not enough input arguments.'); end
if nargin<2, wlen=1024; end
[nchan,nsampl]=size(x);
if nchan>nsampl, error('The signals must be within rows.'); end
if wlen~=4*floor(wlen/4), error('The window length must be a multiple of 4.'); end
if mod(nsampl,wlen/2) ~= 0, nsampl = floor(nsampl/wlen*2)*wlen/2; x = x(:,1:nsampl); fprintf('Input signal truncated at %d samples with respect to window length used\n',nsampl); end
%%% Computing STFT coefficients %%%
% Defining sine window
win=sin((.5:wlen-.5)/wlen*pi).';
nfram=nsampl/wlen*2-1;
nbin=wlen/2+1;
X = complex(double(zeros(nbin,nfram,nchan)));
for i=1:nchan,
for t=0:nfram-1,
% Framing
frame=x(i,t*wlen/2+1:t*wlen/2+wlen).'.*win;
% FFT
fframe=fft(frame);
X(:,t+1,i)=fframe(1:nbin);
end
end
return;
\ No newline at end of file
% anechoic moving speech source
load testSig-180_175.mat
[probability DOAazimuth time] = DOASVMloc1(sig,fs);
figure
imagesc(time,DOAazimuth,probability)
xlabel('Time (s)')
ylabel('DOA (�)')
function w = hanning_(n);
%
% function w = hanning_(n);
%
% returns the n-point Hanning window in a column vector.
% Formular is 0.5 - 0.5*cos(2*pi*(0:N-1)'/N) where N
% equals fix(n).
%
% See Also : hanning (sig. proc. toolbox)
%
% author / date : jens-e. appell / 1.95
%
N = fix(n);
w = .5*(1-cos(2*pi*(0:N-1)'/N));
return;
%%-------------------------------------------------------------------------
%%
%% Copyright (C) 1995 Jens-E. Appell, Carl-von-Ossietzky-Universitat
%%
%% Permission to use, copy, and distribute this software/file and its
%% documentation for any purpose without permission by the author
%% is strictly forbidden.
%%
%% Permission to modify the software is granted, but not the right to
%% distribute the modified code.
%%
%% This software is provided "as is" without express or implied warranty.
%%
%%
%% AUTHOR
%%
%% Jens-E. Appell
%% Carl-von-Ossietzky-Universitat
%% Fachbereich 8, AG Medizinische Physik
%% 26111 Oldenburg
%% Germany
%%
%% e-mail: jens@hinz.physik.uni-oldenburg.de
%%
%%-------------------------------------------------------------------------
% If def.remoteAccessEnable == 1, and this function is existing and in the MATLAB search path,
% than it is called after each retrieved response or simulation step and when the current run
% is finished. It can execute arbitrary code. Usually it is used to let the operator interrupt
% a measurement session or to send messages to the subject.
% The given example interrupts the measurement session after the current run is finished.
% The subject is ask to confirm and to contact the operator. The function might be copied
% via network access to the folder from where the measurement is started at any time.
global work
% send message to response window when run is terminated
if ( work.terminate == 1 )
% open the response window if not existing
if isempty( findobj('Tag','afc_win') )
afc_win('open');
end
h=findobj('Tag','afc_win');
hm = findobj('Tag','afc_message');
% display a message
string_msg = {'Experiment interrupted by Operator.', ...
'Press any button to end.', ...
'Please contact your Operator'};
set(hm,'string', string_msg );
% let the window accept any button ans set userdata to 0 when pressed
set(h,'Userdata',4);
% wait for any button to be pressed
waitfor(h,'UserData',0);
% tell afc to exit
work.abortall = 1;
work.terminate = 1;
% close the window if winEnable is off
if ( def.afcwinEnable == 0 )
afc_win('close');
end
end
% revision 0.94.4 beta, modified 19/04/02
function extern_exampleTDT(action)
global def
global work
switch action
case 'open'
% call a specific hardware here
% TDTInit(work.userpar1);
% just play some zeros to get the DAC running
if ( def.soundEnable > 0 )
afc_sound('play_warmup_zeros');
end
%disp('hardware open');
case 'initialize'
%disp('hardware initialize');
case 'reset'
%disp('hardware reset');
case 'close'
%disp('hardware close');
% shut the hardware down
% TDTClose(work.userpar1);
otherwise
disp('hardware unknown action');
end
% eof
\ No newline at end of file
% revision 1.40.1
function returnValue = sound_exampleWavplay( action )
% SE 16.03.2014
% custom sound function interface
% must support the following actions
% (leave empty if nothing should be done):
% 'init'
% 'close'
% 'play_warmup_zeros'
% 'bgloopwav_restart'
% 'play'
% 'isSoundmex'
% 'isSoundmexMarking'
% 'blockWhilePlaying'
% 'stop'
global work
global def
returnValue = [];
% switch the wanted action
switch ( action )
case 'init'
% some specifics for audioplayer
if ( work.matlabVersionRelease > 7.1 )
def.markIntervalDelay = 0.5; % audioplayer not suited for button marking
else
def.markIntervalDelay = 0;
end
def.deviceID = -1;
case 'close'
% clear audioplayer handle in base
evalin('base','clear afc_ap_handle');
case 'play_warmup_zeros'
ap_handle = audioplayer(zeros(100,2), def.samplerate, def.bits, def.deviceID);
assignin('base','afc_ap_handle',ap_handle);
evalin('base','play(afc_ap_handle);');
case 'bgloopwav_restart'
case 'play'
if ( def.checkOutputClip > 0 )
clipped = max(max(abs(work.out)));
if ( clipped > 1 )
switch def.checkOutputClip
case 1
%warning('AFC:clipped', 'Output stimulus clipped');
warning('Output stimulus clipped');
case 2
%error('AFC:clipped', 'Output stimulus clipped');
error('Output stimulus clipped');
end
end
end
% 1.001 use actual output sig len for blocking
work.blockButtonTime = size(work.out,1)/def.samplerate+0.1; % plus 0.1 sec safety margin
ap_handle = audioplayer(work.out, def.samplerate, def.bits, def.deviceID);
assignin('base','afc_ap_handle',ap_handle);
evalin('base','play(afc_ap_handle);');
if ( def.markIntervalDelay > 0 )
pause( def.markIntervalDelay )
end
case 'isSoundmex'
returnValue = 0;
case 'isSoundmexMarking'
returnValue = 0;
case 'blockWhilePlaying'
% needed for secure blocking in release 6.5
pause(0.01);
case 'stop'
otherwise
warning('sound_exampleAudioplayer: unknown action specified');
end
% eof
\ No newline at end of file
% revision 1.40.1
function returnValue = sound_exampleWavplay( action )
% SE 16.03.2014
% custom sound function interface
% must support the following actions
% (leave empty if nothing should be done):
% 'init'
% 'close'
% 'play_warmup_zeros'
% 'bgloopwav_restart'
% 'play'
% 'isSoundmex'
% 'isSoundmexMarking'
% 'blockWhilePlaying'
% 'stop'
global work
global def
returnValue = [];
% switch the wanted action
switch ( action )
case 'init'
case 'close'
case 'play_warmup_zeros'
wavplay(zeros(100,2),def.samplerate,'async');
case 'bgloopwav_restart'
case 'play'
if ( def.checkOutputClip > 0 )
clipped = max(max(abs(work.out)));
if ( clipped > 1 )
switch def.checkOutputClip
case 1
%warning('AFC:clipped', 'Output stimulus clipped');
warning('Output stimulus clipped');
case 2
%error('AFC:clipped', 'Output stimulus clipped');
error('Output stimulus clipped');
end
end
end
% 1.001 use actual output sig len for blocking
work.blockButtonTime = size(work.out,1)/def.samplerate+0.1; % plus 0.1 sec safety margin
wavplay(work.out,def.samplerate,'async');
if ( def.markIntervalDelay > 0 )
pause( def.markIntervalDelay )
end
case 'isSoundmex'
returnValue = 0;
case 'isSoundmexMarking'
returnValue = 0;
case 'blockWhilePlaying'
% needed for secure blocking in release 6.5
pause(0.01);
case 'stop'
otherwise
warning('sound_exampleSoundmexpro: unknown action specified');
end
% eof
\ No newline at end of file
%
% afc - top level function for alternative forced choice procedure -
%
% Usage: afc( commandName, [experiment/session], [subject], [userpar1, ..., userparN, condition], [key1, value1])
%
% commandName = string specifying mode of operation: 'help', 'about',
% 'setup', 'version', 'main', 'session', 'demo', 'selftest'
%
% if commandName == 'main'
%
% experiment = string containing the name of the experiment to be performed
% subject = string containing the name of the subject
% userpar(s) = optional string(s) containing user defined parameter(s)
% condition = optional string containing experimental condition,
% always the last optional argument after subject and before optional key1
% key1 = optional string 'runs', expects value1 as next argument
% value1 = optional integer following key1, number of runs before AFC stops
%
% example: afc('main','example','ab','cd1')
% starts the experiment example for subject ab and condition 1.
%
% if commandName == 'session'
%
% session = string containing the name of the session to be performed
%
% example: afc('session','example')
% starts the session defined in ASCII file session_example.dat
%
% other commandNames
%'help' opens this help
%'setup' opens configuration window
%'about' opens about window
%'demo' starts demo
%'selftest' starts selftest
%'version' displays version
%
% See also help afc_main, afc_session, example_cfg, example_set, example_user
%
%------------------------------------------------------------------------------
% AFC for Mathwork’s MATLAB
%
% Version 1.40.0
%
% Author(s): Stephan Ewert
%
% Copyright (c) 1999-2014, Stephan Ewert.
% All rights reserved.
%
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0).
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-nd/4.0/ or send a letter to Creative
% Commons, 444 Castro Street, Suite 900, Mountain View, California, 94041, USA.
%------------------------------------------------------------------------------
% last modified 19-04-2013 13:02:51
% $Revision: 1.20.0 beta$ $Date: 06.06.2011 13:41 $
function afc(varargin)
% figure out version
%mstr = version;
%matlabVersion = str2num(mstr(1));
if ( isempty(varargin) | strcmp( varargin{1}, 'help' ) )
% r1400 remove about evalin('base','help afc; afc(''about'');');