Commit 023f7650 authored by VINCENT Emmanuel's avatar VINCENT Emmanuel
Browse files

Version 2 - March 21, 2012

parent 0c774357
The Matlab M-files bss_locate_clus.m, bss_locate_spec.m, qstft.m and stft_multi.m are authored by Charles Blandin and/or Emmanuel Vincent and distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
The Matlab MAT-files in the directory "filters" are authored by Charles Blandin and distributed under the the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 2.0 license (http://creativecommons.org/licenses/by-nc-sa/2.0/).
The WAV audio files music_s1.wav, music_s2.wav, music_s3.wav, music_s4.wav and music_s5.wav are authored by Shannon Hurley (http://amiestreet.com/artist/shannon-hurley/) and distributed under the the terms of the Creative Commons Attribution-NonCommercial 3.0 license (http://creativecommons.org/licenses/by-nc/3.0/).
The WAV audio file music_s6.wav is authored by The Ultimate NZ Tour (http://newzealandeducated.com/contest/ultimate-nz-tour-i/) and distributed under the the terms of the Creative Commons Attribution-Noncommercial-ShareAlike 3.0 license (http://creativecommons.org/licenses/by-nc-sa/3.0/).
All other WAV audio files are distributed under the the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 2.0 license (http://creativecommons.org/licenses/by-nc-sa/2.0/).
# bss_locate
This repository contains version 2 of BSS Locate dated from March 21, 2012, and previously available at http://bass-db.gforge.inria.fr/bss_locate/.
# Purpose
BSS Locate is a Matlab toolbox that estimates the **Time Differences of Arrival (TDOAs)** of **multiple sources** in a **stereo audio signal** recorded by a pair of **omnidirectional microphones**. The number of sources and the distance between the two microphones are assumed to be known but are not restricted to a given range.
BSS Locate implements 12 different source localization methods:
- 8 methods based on the detection of the peaks of an angular spectrum (GCC-PHAT, GCC-NONLIN, MUSIC and several SNR-based spectra)
- 4 methods based on clustering of the time-frequency bins
GCC-NONLIN provides the best results over a wide range of configurations.
# Usage
Two categories of source localization methods are implemented:
- angular spectrum-based methods (bss_locate_spec.m)
- clustering-based methods (bss_locate_clus.m)
In order to apply the most accurate method to date (GCC-NONLIN), just type
tau = bss_locate_spec(x, fs, d, nsrc)
where x is an nsampl x 2 matrix containing a stereo signal, fs is the sampling frequency in Hz, d is the microphone spacing in meters, and nsrc is the number of sources. The output tau is then the 1 x nsrc vector of estimated TDOAs in seconds.
For information about the syntax, type
help bss_locate_spec
help bss_locate_clus
# Reproductibility
The toolbox comes with the source signals and the mixing filters used for evaluation in [[1]](#1).
The filters are named <J>src_<RT>_<d>_<r>_config<k>.mat where <J> is the
number of sources, <RT> the reverberation time, <d> the microphone spacing,
<r> the distance between the sources and the center of the microphone pair and
<k> the index of one configuration satisfying the above characteristics.
To generate the signals used for evaluation, add the directories sources/ and filters/ to the Matlab path, and type
x = mix(fname, stype);
where fname is one of the mixing filter filenames in the directory filters/
and stype is 'female', 'male' or 'music'
To evaluate the results of your algorithm, type
taue = true_tdoas(fname);
[R, P, F, Acc] = eval_tdoas(taue, tau, d, 5);
# Reference
<a id="1">[1]</a> C. Blandin, A. Ozerov, and E. Vincent, ["Multi-source TDOA estimation in reverberant audio using angular spectra and clustering"](https://hal.inria.fr/inria-00630994v2/document), *Signal Processing* 92, pp. 1950-1960, 2012.
# License
The Matlab M-files bss_locate_clus.m, bss_locate_spec.m, qstft.m and stft_multi.m are authored by Charles Blandin and/or Emmanuel Vincent and distributed under the terms of the [GNU Public License version 3](http://www.gnu.org/licenses/gpl.txt).
The Matlab MAT-files in the directory "filters" are authored by Charles Blandin and distributed under the the terms of the [Creative Commons Attribution-NonCommercial-ShareAlike 2.0 license](http://creativecommons.org/licenses/by-nc-sa/2.0/).
The WAV audio files music_s1.wav, music_s2.wav, music_s3.wav, music_s4.wav and music_s5.wav are authored by [Shannon Hurley](http://amiestreet.com/artist/shannon-hurley/) and distributed under the the terms of the [Creative Commons Attribution-NonCommercial 3.0 license](http://creativecommons.org/licenses/by-nc/3.0/).
The WAV audio file music_s6.wav is authored by [The Ultimate NZ Tour](http://newzealandeducated.com/contest/ultimate-nz-tour-i/) and distributed under the the terms of the [Creative Commons Attribution-Noncommercial-ShareAlike 3.0 license](http://creativecommons.org/licenses/by-nc-sa/3.0/).
All other WAV audio files are distributed under the the terms of the [Creative Commons Attribution-NonCommercial-ShareAlike 2.0 license](http://creativecommons.org/licenses/by-nc-sa/2.0/).
\ No newline at end of file
function [tau, crit] = bss_locate_clus(x, fs, d, tau_init, method, noise, tau_grid)
% BSS_LOCATE_CLUS Estimation of the source TDOAs in a stereo convolutive
% mixture using a clustering method
%
% [tau, crit] = bss_locate_clus(x, fs, d, tau_init, method, tau_grid)
% [tau, crit] = bss_locate_clus(x, fs, d, tau_init, method, noise, tau_grid)
%
% Inputs:
% x: nsampl x 2 matrix containing a stereo mixture signal
% fs: sampling frequency in Hz
% d: microphone spacing in meters ([] for Sawada)
% tau_init: 1 x nsrc vector of initial TDOAs in seconds
% method: clustering method: 'Sawada', 'Izumi', 'EM-predom' (default),
% 'EM-multi'
% noise: noise variance tying for Izumi and EM-multi: 'TF', 'F', 'const'
% (default) or 'no'
% tau_grid: 1 x ngrid vector of possible TDOAs in seconds (default: 181
% values linearly spaced between -d/343 and d/343)
%
% Outputs:
% tau: 1 x nsrc vector of estimated TDOAs in seconds
% crit: value of the objective function (log-likelihood or negative
% squared Euclidean distance)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2010-2011 Emmanuel Vincent and Charles Blandin
% This software is distributed under the terms of the GNU Public License
% version 3 (http://www.gnu.org/licenses/gpl.txt)
% If you find it useful, please cite the following reference:
% Charles Blandin, Emmanuel Vincent and Alexey Ozerov, "Multi-source TDOA
% estimation in reverberant audio using angular spectra and clustering",
% Signal Processing 92, pp. 1950-1960, 2012.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Errors and default values %%%
if nargin<4, error('Not enough input arguments.'); end
[nsampl,nchan] = size(x);
if nchan>nsampl, error('The input signal must be in columns.'); end
if nchan~=2, error('The input signal must be stereo.'); end
if nargin < 5, method = 'EM-predom'; end
if ~any(strcmp(method, {'Sawada' 'Izumi' 'EM-predom' 'EM-multi'})), error('Unknown clustering method.'); end
if any(strcmp(method, {'Izumi' 'EM-multi'})),
if nargin < 6, noise = 'const'; end
if ~any(strcmp(noise, {'TF' 'F' 'const' 'no'})), error('Unknown noise tying setting.'); end
if nargin < 7, tau_grid = linspace(-d/343, d/343, 181); end
else
if nargin < 6,
tau_grid = linspace(-d/343, d/343, 181);
else
tau_grid = noise;
end
end
%%% Time-frequency transform %%%
wlen = 1024;
if isempty(strfind(method, 'EM')),
% Linear transform
X = stft_multi(x.',wlen);
X = X(2:end,:,:);
else
% Quadratic transform
hatRxx = qstft(x,fs,wlen);
hatRxx = permute(hatRxx(:,:,2:end,:),[3 4 1 2]);
end
f = fs/wlen*(1:wlen/2).';
%%% Clustering %%%
switch method
case 'Sawada'
[tau, crit] = clus_sawada(X, f, tau_init, tau_grid);
case 'Izumi'
[tau, crit] = clus_izumi(X, f, d, tau_init, noise, tau_grid);
case 'EM-predom'
[tau, crit] = em_predom(hatRxx, f, d, tau_init, tau_grid);
case 'EM-multi'
[tau, crit] = em_multi(hatRxx, f, d, tau_init, noise, tau_grid);
end
return;
function [tau, crit] = clus_sawada(X, f, tau_init, tau_grid)
% CLUS_SAWADA Estimation of the source TDOAs in a stereo convolutive
% mixture using the clustering algorithm in Section V.A of
% H. Sawada, S. Araki, R. Mukai, S. Makino, "Grouping separated frequency
% components by estimating propagation model parameters in frequency-domain
% blind source separation", IEEE Transactions on Audio, Speech, and
% Language Processing, 15(5):1592–1604, 2007.
%
%
% tau = clus_sawada(X, f, tau_init, tau_grid)
%
% Inputs:
% X : nbin x nfram x 2 matrix containing the STFT coefficients of the input
% signal in all time-frequency bins
% f: nbin x 1 vector containing the center frequency of each frequency bin
% in Hz
% tau_init: 1 x nsrc vector of initial TDOAs in seconds
% tau_grid: 1 x ngrid vector of possible TDOAs in seconds
%
% Output:
% tau: 1 x nsrc vector of estimated TDOAs in seconds
% crit: value of the negative squared Euclidean distance
[nbin,nfram] = size(X(:,:,1));
nsrc = length(tau_init);
ngrid = length(tau_grid);
%%% Input normalization %%%
% Phase normalization
Xtilde = zeros(nbin,nfram,2);
Xtilde(:,:,1) = abs(X(:,:,1));
Xtilde(:,:,2) = X(:,:,2).* conj(X(:,:,1)) ./ Xtilde(:,:,1);
% Amplitude normalization
norm = sqrt(abs(Xtilde(:,:,1)).^2 + abs(Xtilde(:,:,2)).^2);
Xtilde = Xtilde./repmat(norm,[1 1 2]);
%%% Clustering %%%
tau = tau_init;
maxiter = 100;
iter = 1;
converged = false;
while ((iter <= maxiter) && ~converged),
% Computing the theoretical frequency response in (17)
D = zeros(nbin,2,nsrc);
D(:,1,:) = sqrt(.5)*ones(nbin,1,nsrc);
D(:,2,:) = sqrt(.5)*exp(-2*1i*pi*f*tau);
% Clustering in (48)
crit_clus = zeros(nbin,nfram,nsrc);
for j=1:nsrc,
crit_clus(:,:,j) = abs(Xtilde(:,:,1)-repmat(D(:,1,j),1,nfram)).^2 + abs(Xtilde(:,:,2)-repmat(D(:,2,j),1,nfram)).^2;
end
[crit_clus, C] = min(crit_clus, [], 3);
% Updating the TDOAs in (49)
prevtau = tau;
T = real(repmat(Xtilde(:,:,2),[1,1,ngrid]).*permute(repmat(exp(2*1i*pi*f*tau_grid),[1,1,nfram]),[1,3,2]));
for j = 1:nsrc,
crit_tau = sum(sum(repmat(double(C==j),[1,1,ngrid]).*T,1),2);
[crit_tau, ind] = max(crit_tau);
tau(j) = tau_grid(ind);
end
% TDOA convergence check
converged = all(tau == prevtau);
iter = iter+1;
end
crit = -sum(sum(crit_clus)) / (nbin*nfram);
return;
function [tau, loglik] = clus_izumi(X, f, d, tau_init, noise, tau_grid)
% CLUS_IZUMI Estimation of the source TDOAs in a stereo convolutive
% mixture using the clustering algorithm in
% Y. Izumi, N. Ono, S. Sagayama, "Sparseness-based 2ch BSS using the EM
% algorithm in reverberant environment", in Proc IEEE Workshop on
% Applications of Signal Processing to Audio and Acoustics (WASPAA), pp.
% 147–150, 2007.
%
% tau = clus_izumi(X, f, d, tau_init, tau_grid)
%
% Inputs:
% X : nbin x nfram x 2 matrix containing the STFT coefficients of the input
% signal in all time-frequency bins
% f: nbin x 1 vector containing the center frequency of each frequency bin
% in Hz
% d: microphone spacing in meters
% tau_init: 1 x nsrc vector of initial TDOAs in seconds
% noise: noise variance tying: 'TF', 'F', 'const' or 'no'
% tau_grid: 1 x ngrid vector of possible TDOAs in seconds
%
% Outputs:
% tau: 1 x nsrc vector of estimated TDOAs in seconds
% loglik: value of the log-likelihood
[nbin,nfram] = size(X(:,:,1));
nsrc = length(tau_init);
ngrid = length(tau_grid);
R11 = abs(X(:,:,1)).^2;
R12 = X(:,:,1) .* conj(X(:,:,2));
R22 = abs(X(:,:,2)).^2;
%%% Initialization %%%
tau = tau_init;
if strcmp(noise, 'no'),
vb = 1e-8 * ones(nbin,nfram);
else
vb = ones(nbin,nfram);
end
c = 343;
SINC = repmat(sinc(2*f*d/c),1,nfram);
SINC2 = SINC.^2;
maxiter = 100;
iter = 1;
converged = false;
loglik=-inf;
while ((iter <= maxiter) && ~converged),
%%% E-Step %%%
% Computing the log-likelihood in (7) (erroneous equation in the paper)
M = zeros(nbin,nfram,nsrc);
for j = 1:nsrc,
EXP = repmat(exp(-2*1i*pi*tau(j)*f),1,nfram);
M(:,:,j) = (R11+R22-2*real(R12.*SINC))./(1-SINC2) - ((R11 + R22).*(ones(nbin,nfram) + SINC2 - 2*real(EXP).*SINC) + 2*real(R12.*(EXP - 2*SINC + SINC2.*conj(EXP)))) ./ (2*(1 - real(EXP).*SINC).*(1-SINC2));
end
logPx = repmat(-log(pi^2*vb.^2.*(1-SINC2)),[1 1 nsrc]) - M./repmat(vb,[1 1 nsrc]);
% Log-likelihood convergence check in (9)
prevloglik = loglik;
logPxmax = max(logPx,[],3);
gamma = exp(logPx-repmat(logPxmax,[1 1 nsrc]));
loglik = sum(sum(log(1/nsrc*sum(gamma,3))+logPxmax)) / (nbin*nfram);
converged = (loglik-prevloglik < 1e-2);
% Computing the source posterior probabilities in (15)
gamma = gamma ./ repmat(sum(gamma,3),[1 1 nsrc]);
%%% M-step %%%
% Updating the noise variance in (17) (erroneous equation in the paper)
switch noise
case 'TF'
vb = max(.5 * sum(gamma.*M,3), realmin);
case 'F'
vb = repmat(max(.5 * mean(sum(gamma.*M,3),2), realmin),1,nfram);
case 'const'
vb = repmat(max(.5 * mean(mean(sum(gamma.*M,3))), realmin),nbin,nfram);
end
if converged,
prevtau=tau;
% Updating the TDOA in (16)
logPx = zeros(nbin,nfram,ngrid);
for ind = 1:ngrid,
EXP = repmat(exp(-2*1i*pi*tau_grid(ind)*f),1,nfram);
M = (R11+R22-2*real(R12.*SINC))./(1-SINC2) - ((R11 + R22).*(ones(nbin,nfram) + SINC2 - 2*real(EXP).*SINC) + 2*real(R12.*(EXP - 2*SINC + SINC2.*conj(EXP)))) ./ (2*(1 - real(EXP).*SINC).*(1-SINC2));
logPx(:,:,ind) = -log(pi^2*vb.^2.*(1-SINC2)) - M./vb;
end
for j = 1:nsrc,
crit_tau = zeros(ngrid,1);
for ind = 1:ngrid,
crit_tau(ind) = sum(sum(gamma(:,:,j) .* logPx(:,:,ind)));
end
[crit_tau, ind] = max(crit_tau);
tau(j) = tau_grid(ind);
end
% TDOA convergence check
converged = all(tau == prevtau);
end
iter = iter+1;
end
return;
function [tau, loglik] = em_predom(hatRxx, f, d, tau_init, tau_grid)
% EM_PREDOM Estimation of the source TDOAs in a stereo convolutive
% mixture by EM-based clustering assuming a single predominant source in
% each time-frequency bin
%
% tau = em_predom(hatRxx, f, d, tau_init, tau_grid)
%
% Inputs:
% hatRxx : nbin x nfram x 2 x 2 array containing the spatial covariance
% matrices of the input signal in all time-frequency bins
% f: nbin x 1 vector containing the center frequency of each frequency bin
% in Hz
% d: microphone spacing in meters
% tau_init: 1 x nsrc vector of initial TDOAs in seconds
% tau_grid: 1 x ngrid vector of possible TDOAs in seconds
%
% Outputs:
% tau: 1 x nsrc vector of estimated TDOAs in seconds
% loglik: value of the log-likelihood
[nbin,nfram] = size(hatRxx(:,:,1,1));
nsrc = length(tau_init);
ngrid = length(tau_grid);
R11 = real(hatRxx(:,:,1,1));
R12 = hatRxx(:,:,1,2);
R21 = hatRxx(:,:,2,1);
R22 = real(hatRxx(:,:,2,2));
c = 343;
SINC = sinc(2*f*d/c);
SINC2 = SINC.^2;
%%% Computing the source and noise variance and the log-likelihood in all directions %%%
logPx = zeros(nbin,nfram,ngrid);
for ind = 1:ngrid,
% Computing inv(A) = [invA11 invA12; conj(invA11) -invA12]
EXP = exp(-2*1i*pi*tau_grid(ind)*f);
P = SINC .* EXP;
invA11 = sqrt(.5)./(1-real(P)).*(1-conj(P));
invA12 = -(1-P)./(SINC-EXP).*invA11;
% Computing inv(Lambda) = [.5 invL12; 0 invL22]
DEN = .5./(1-2*real(P)+SINC2);
invL12 = (SINC2-1).*DEN;
invL22 = 2*(1-real(P)).*DEN;
% Computing vs and vb without nonnegativity constraint
ARA1 = repmat(abs(invA11).^2,1,nfram).*R11 + repmat(abs(invA12).^2,1,nfram).*R22;
ARA2 = ARA1 - 2 * real(repmat(invA11.*invA12,1,nfram).*R21);
ARA1 = ARA1 + 2 * real(repmat(invA11.*conj(invA12),1,nfram).*R12);
vs = .5*ARA1 + repmat(invL12,1,nfram).*ARA2;
vb = repmat(invL22,1,nfram).*ARA2;
% Enforcing the nonnegativity constraint (on vs and vb)
neg = (vs < 0) | (vb < 0);
vs(neg) = 0;
vb0 = repmat(.5./(1-SINC2),1,nfram).*(R11 + R22 - 2*real(R21).*repmat(SINC,1,nfram));
vb(neg) = vb0(neg);
% Computing the log-likelihood
detRxx = vb.^2 .* repmat(1-SINC2,1,nfram) + 2*vs.*vb.*repmat(1-SINC.*real(EXP),1,nfram);
logPx(:,:,ind) = -log(pi^2*detRxx) - ((R11+R22).*(vs+vb) - 2*real(R12.*(vs.*repmat(EXP,1,nfram) + vb.*repmat(SINC,1,nfram))))./detRxx;
end
tau_ind = round(interp1(tau_grid,1:ngrid,tau_init,'linear'));
maxiter = 100;
iter = 1;
converged = false;
while ((iter <= maxiter) && ~converged),
%%% E-Step %%%
% Computing gamma
logPxmax = max(logPx(:,:,tau_ind),[],3);
gamma = exp(logPx(:,:,tau_ind)-repmat(logPxmax,[1 1 nsrc]));
gamma = gamma ./ repmat(sum(gamma,3),[1 1 nsrc]);
% Log-likelihood (unnecessary)
% loglik = sum(sum(log(1/nsrc*sum(exp(logPx(:,:,tau_ind)),3))))/(nbin*nfram);
%%% M-step %%%
% Updating tau
prevtau_ind=tau_ind;
for j = 1:nsrc,
crit_tau = zeros(ngrid,1);
for ind = 1:ngrid,
crit_tau(ind) = sum(sum(gamma(:,:,j) .* logPx(:,:,ind)));
end
[crit_tau, tau_ind(j)] = max(crit_tau);
end
% TDOA convergence check
converged = all(tau_ind == prevtau_ind);
iter = iter+1;
end
tau = tau_grid(tau_ind);
return;
function [tau, loglik] = em_multi(hatRxx, f, d, tau_init, noise, tau_grid)
% EM_MULTI Estimation of the source TDOAs in a stereo convolutive
% mixture by EM-based clustering assuming multiple sources in each
% time-frequency bin
%
% tau = em_multi(hatRxx, f, d, tau_init, noise, tau_grid)
%
% Inputs:
% hatRxx : nbin x nfram x 2 x 2 array containing the spatial covariance
% matrices of the input signal in all time-frequency bins
% f: nbin x 1 vector containing the center frequency of each frequency bin
% in Hz
% d: microphone spacing in meters
% tau_init: 1 x nsrc vector of initial TDOAs in seconds
% noise: noise variance tying: 'TF', 'F', 'const' or 'no'
% tau_grid: 1 x ngrid vector of possible TDOAs in seconds
%
% Outputs:
% tau: 1 x nsrc vector of estimated TDOAs in seconds
% loglik: value of the log-likelihood
[nbin,nfram] = size(hatRxx(:,:,1,1));
nsrc = length(tau_init);
%%% Initialization %%%
% Defining Psi and inv(Psi)
c = 343;
Psi=ones(nbin,2,2);
Psi(:,1,2)=sinc(2*f*d/c);
Psi(:,2,1)=Psi(:,1,2);
detPsi=1-Psi(:,1,2).^2;
invPsi=zeros(nbin,2,2);
invPsi(:,1,1)=1./detPsi;
invPsi(:,1,2)=-Psi(:,1,2)./detPsi;
invPsi(:,2,1)=invPsi(:,1,2);
invPsi(:,2,2)=invPsi(:,1,1);
% Initializing other variables
tau = tau_init;
D = zeros(nbin,2,nsrc);
D(:,1,:) = ones(nbin,1,nsrc);
D(:,2,:) = exp(-2*1i*pi*f*tau);
vs = 1e-3 * ones(nbin,nfram,nsrc);
if strcmp(noise, 'no'),
vb = 1e-8 * ones(nbin,nfram);
else
vb = ones(nbin,nfram);
end
Gs = zeros(nbin,nfram,nsrc,2);
hatRxs = zeros(nbin,nfram,2,nsrc);
hatRss = zeros(nbin,nfram,nsrc,nsrc);
maxiter = 100;
iter = 1;
converged = false;
loglik=-inf;
while ((iter <= maxiter) && ~converged),
%%% E-Step %%%
% Computing Rxx
Rxx = repmat(vb,[1 1 2 2]).*repmat(reshape(Psi,nbin,1,2,2),[1 nfram 1 1]);
for j = 1:nsrc,
Rxx(:,:,1,1) = Rxx(:,:,1,1) + repmat(abs(D(:,1,j)).^2,1,nfram) .* vs(:,:,j);
Rxx(:,:,1,2) = Rxx(:,:,1,2) + repmat(D(:,1,j) .* conj(D(:,2,j)),1,nfram) .* vs(:,:,j);
Rxx(:,:,2,2) = Rxx(:,:,2,2) + repmat(abs(D(:,2,j)).^2,1,nfram) .* vs(:,:,j);
end
Rxx(:,:,2,1) = conj(Rxx(:,:,1,2));
% Computing det(Rxx) and inv(Rxx)
detRxx = real(Rxx(:,:,1,1) .* Rxx(:,:,2,2) - abs(Rxx(:,:,1,2)).^2);
invRxx(:,:,1,1) = Rxx(:,:,2,2) ./ detRxx;
invRxx(:,:,1,2) = - Rxx(:,:,1,2) ./ detRxx;
invRxx(:,:,2,1) = conj(invRxx(:,:,1,2));
invRxx(:,:,2,2) = Rxx(:,:,1,1) ./ detRxx;
% Log-likelihood convergence test
prevloglik = loglik;
loglik = -sum(sum(real(log(detRxx * pi^2) + invRxx(:,:,1,1) .* hatRxx(:,:,1,1) + invRxx(:,:,2,2) .* hatRxx(:,:,2,2) + 2 * invRxx(:,:,1,2) .* hatRxx(:,:,2,1)) )) / (nbin * nfram);
converged = (loglik-prevloglik < 1e-2);
% Computing Gs and hatRxs
for j = 1:nsrc,
Gs(:,:,j,1) = (repmat(conj(D(:,1,j)),1,nfram) .* invRxx(:,:,1,1) + repmat(conj(D(:,2,j)),1,nfram) .* invRxx(:,:,2,1)) .* vs(:,:,j);
Gs(:,:,j,2) = (repmat(conj(D(:,1,j)),1,nfram) .* invRxx(:,:,1,2) + repmat(conj(D(:,2,j)),1,nfram) .* invRxx(:,:,2,2)) .* vs(:,:,j);
hatRxs(:,:,1,j) = hatRxx(:,:,1,1) .* conj(Gs(:,:,j,1)) + hatRxx(:,:,1,2) .* conj(Gs(:,:,j,2));
hatRxs(:,:,2,j) = hatRxx(:,:,2,1) .* conj(Gs(:,:,j,1)) + hatRxx(:,:,2,2) .* conj(Gs(:,:,j,2));
end
% Computing hatRss
for j1 = 1:nsrc,
for j2 = j1:nsrc,
hatRss(:,:,j1,j2) = Gs(:,:,j1,1) .* (hatRxs(:,:,1,j2) - repmat(D(:,1,j2),1,nfram) .* vs(:,:,j2)) + Gs(:,:,j1,2) .* (hatRxs(:,:,2,j2) - repmat(D(:,2,j2),1,nfram) .* vs(:,:,j2));
end
for j2 = 1:j1-1,
hatRss(:,:,j1,j2) = conj(hatRss(:,:,j2,j1));
end
hatRss(:,:,j1,j1) = hatRss(:,:,j1,j1) + vs(:,:,j1);
end
%%% M-step %%%
% Updating vs
for j = 1:nsrc,
vs(:,:,j) = real(hatRss(:,:,j,j));
end
% Updating vb
E = hatRxx;
for j1 = 1:nsrc,
E(:,:,1,1) = E(:,:,1,1) - 2 * real(repmat(D(:,1,j1),1,nfram) .* conj(hatRxs(:,:,1,j1))) + repmat(abs(D(:,1,j1)).^2,1,nfram) .* hatRss(:,:,j1,j1);
E(:,:,1,2) = E(:,:,1,2) - repmat(D(:,1,j1),1,nfram) .* conj(hatRxs(:,:,2,j1)) - repmat(conj(D(:,2,j1)),1,nfram) .* hatRxs(:,:,1,j1) + repmat(D(:,1,j1),1,nfram) .* hatRss(:,:,j1,j1) .* repmat(conj(D(:,2,j1)),1,nfram);
E(:,:,2,2) = E(:,:,2,2) - 2 * real(repmat(D(:,2,j1),1,nfram) .* conj(hatRxs(:,:,2,j1))) + repmat(abs(D(:,2,j1)).^2,1,nfram) .* hatRss(:,:,j1,j1);
for j2= 1:j1-1,
E(:,:,1,1) = E(:,:,1,1) + 2 * real(repmat(D(:,1,j1),1,nfram) .* hatRss(:,:,j1,j2) .* repmat(conj(D(:,1,j2)),1,nfram));
E(:,:,1,2) = E(:,:,1,2) + repmat(D(:,1,j1),1,nfram) .* hatRss(:,:,j1,j2) .* repmat(conj(D(:,2,j2)),1,nfram) + repmat(D(:,1,j2),1,nfram) .* hatRss(:,:,j2,j1) .* repmat(conj(D(:,2,j1)),1,nfram);
E(:,:,2,2) = E(:,:,2,2) + 2 * real(repmat(D(:,2,j1),1,nfram) .* hatRss(:,:,j1,j2) .* repmat(conj(D(:,2,j2)),1,nfram));
end
end
E(:,:,2,1) = conj(E(:,:,1,2));
M = real(E(:,:,1,1) .* repmat(invPsi(:,1,1),1,nfram) + 2 * E(:,:,1,2) .* repmat(invPsi(:,2,1),1,nfram) + E(:,:,2,2) .* repmat(invPsi(:,2,2),1,nfram));
switch noise
case 'TF'
vb = max(.5 * M, realmin);
case 'F'
vb = repmat(max(.5 * mean(M,2), realmin),1,nfram);
case 'const'
vb = repmat(max(.5 * mean(mean(M)), realmin),nbin,nfram);
end
if converged,
prevtau=tau;
% Updating tau
for j = 1:nsrc,
crit_tau = zeros(ngrid,1);
for ind = 1:ngrid,
D(:,2,j) = exp(-2*1i*pi*f*tau_grid(ind));
E = hatRxx;
for j1 = 1:nsrc,
E(:,:,1,1) = E(:,:,1,1) - 2 * real(repmat(D(:,1,j1),1,nfram) .* conj(hatRxs(:,:,1,j1))) + repmat(abs(D(:,1,j1)).^2,1,nfram) .* hatRss(:,:,j1,j1);
E(:,:,1,2) = E(:,:,1,2) - repmat(D(:,1,j1),1,nfram) .* conj(hatRxs(:,:,2,j1)) - repmat(conj(D(:,2,j1)),1,nfram) .* hatRxs(:,:,1,j1) + repmat(D(:,1,j1),1,nfram) .* hatRss(:,:,j1,j1) .* repmat(conj(D(:,2,j1)),1,nfram);
E(:,:,2,2) = E(:,:,2,2) - 2 * real(repmat(D(:,2,