DOASVMloc1.m 2.34 KB
Newer Older
Joanna Luberadzka's avatar
Joanna Luberadzka committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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