diff --git a/Code/capon.m b/Code/capon.m
index 5d5dc9ea66a468284f8c046b10652af309d30aa3..089ee1c259a43a2ea0bfa7d6f368fc3f4fd3f764 100644
--- a/Code/capon.m
+++ b/Code/capon.m
@@ -6,6 +6,8 @@ function P = capon(Y, doa, pos, f)
% doa Directions of Arrival
% pos Positions of sensors
% f Frequency of signal
+%
+% Author Clas Veib�ck
n = size(Y,2);
c = 343;
diff --git a/Code/common/alignSignal.m b/Code/common/alignSignal.m
index a580d03ed928dbdd19b3bc0f9f43de89d2051797..46837e18a905e73b77242bbae8c1fcb12a221c9d 100644
--- a/Code/common/alignSignal.m
+++ b/Code/common/alignSignal.m
@@ -1,9 +1,11 @@
function [aligned, prepad, postpad] = alignSignal(signal,T, trim)
%ALIGNSIGNAL Align signals in discrete time units.
-% Shifts the signals to align according to a given offset for each signal
-% signal The signals to align
-% T Sample offset for each signal
-% trim Trim the end points to retain signal length (default false)
+% Shifts the signals to align according to a given offset for each signal
+% signal The signals to align
+% T Sample offset for each signal
+% trim Trim the end points to retain signal length (default false)
+%
+% Author Clas Veib�ck
N = size(signal,1);
n = length(T);
diff --git a/Code/common/correlationMatrix.m b/Code/common/correlationMatrix.m
index 90848ba2e23f6415afdeb2e780b1146588406e4c..6c6004775313d7984c7fad49dc46666be707c4e9 100644
--- a/Code/common/correlationMatrix.m
+++ b/Code/common/correlationMatrix.m
@@ -1,6 +1,11 @@
function [corr, lags] = correlationMatrix(Y, maxlag, atLags)
%CORRELATIONMATRIX Find best correlation between all combinations of
%columns
+% Y Signals
+% maxlag Maximum number of lags to consider
+% atLags Return correlation matrix for given lags(leave empty otherwise)
+%
+% Author Clas Veib�ck
[N, n] = size(Y);
diff --git a/Code/common/delayAndSumSignal.m b/Code/common/delayAndSumSignal.m
index 58e765cb1bb5e911a379da2a31ff9b1ce3cd330b..7db6e08cbd3eab0e7eccf949afbf025764db54f5 100644
--- a/Code/common/delayAndSumSignal.m
+++ b/Code/common/delayAndSumSignal.m
@@ -1,9 +1,11 @@
function S = delayAndSumSignal(Y, ds)
%DELAYANDSUMSIGNAL Delay and sum signal in continuous time. Units are
%samples, which can be fractional.
-% Gridded interpolant with spline is used.
-% Y Signals to delay and sum
-% ds Number of samples to delay for each signal before summing
+% Gridded interpolant with spline is used.
+% Y Signals to delay and sum
+% ds Number of samples to delay for each signal before summing
+%
+% Author Clas Veib�ck
N= size(Y,1);
n = size(ds,1);
diff --git a/Code/common/intFFT.m b/Code/common/intFFT.m
index edcaee1c44ec9c506dd5a7c645c401b70f4064a0..53edd4831b1e4a20a655f05612037d218ab92dd7 100644
--- a/Code/common/intFFT.m
+++ b/Code/common/intFFT.m
@@ -1,9 +1,11 @@
function sx = intFFT(x, fs)
%INTFFT Integrate a vector in the frequency domain
-% Integrates a vector sampled at the rate fs in the frequency domain.
-% Columnwise if signal is matrix.
-% x Signal vector
-% fs Sampling frequency (fs = 1)
+% Integrates a vector sampled at the rate fs in the frequency domain.
+% Columnwise if signal is matrix.
+% x Signal vector
+% fs Sampling frequency (fs = 1)
+%
+% Author Clas Veib�ck
if nargin < 2
fs = 1;
diff --git a/Code/common/showResults.m b/Code/common/showResults.m
index 72cbefc9d2b6bbd41906c5a89178db0d949fbd51..1e5bdc212d7a546c0ced06074af5d7a9a76df0d3 100644
--- a/Code/common/showResults.m
+++ b/Code/common/showResults.m
@@ -1,13 +1,15 @@
function angles = showResults(angle, time, P, name, doa)
%SHOWRESULTS Displays the results for the segments in time and the grid
%points.
-% An image is shown with the maximum angle for each segment overlayed, as
-% well as a moving median.
-% angle Angles in grid (radians)
-% time Start and end time (uniform sampling assumed in between)
-% P Result to be display (high value -> better match)
-% name Name to display on plot
-% doa True angles (optional)
+% An image is shown with the maximum angle for each segment overlayed, as
+% well as a moving median.
+% angle Angles in grid (radians)
+% time Start and end time (uniform sampling assumed in between)
+% P Result to be display (high value -> better match)
+% name Name to display on plot
+% doa True angles (optional)
+%
+% Author Clas Veib�ck
if nargin < 5; doa = []; end
diff --git a/Code/delayAndSum.m b/Code/delayAndSum.m
index 49d55ee90339eb3b8c60e220395af5e938549ef7..39a2e8f842700f7deb28269c4334f4cc04b1ce66 100644
--- a/Code/delayAndSum.m
+++ b/Code/delayAndSum.m
@@ -1,13 +1,15 @@
function [S, P] = delayAndSum(Y, doa, pos, fs, domain, dir)
%DELAYANDSUM Broadband far-field Delay and Sum beamformer and Doa
%estimator
-% Performs DoA estimation using the Delay and sum method,
+% Performs beamforming and DoA estimation using the Delay and sum method,
% which is broadband, data-independent and far-field.
% Y Signal recordings
% doa Directions of Arrival
% pos Positions of sensors
% fs Sampling frequency of signal
% domain Domain of computation: 'temporal', 'spectral' (default)
+%
+% Author Clas Veib�ck
n = size(Y, 2);
N = size(Y, 1);
diff --git a/Code/diffLinDoA.m b/Code/diffLinDoA.m
index 321e8da9d3fb7a7a98fa0cac686dab13fc7be66d..251e371a5b1ca128d0a749103d187a78da15def6 100644
--- a/Code/diffLinDoA.m
+++ b/Code/diffLinDoA.m
@@ -1,11 +1,13 @@
function [S, P] = diffLinDoA(Y, doa, pos, ndeg)
%DIFFLINDOA Broadband far-field diffLinDoA beamformer and DoA estimator
-% Performs DoA estimation using the diffLinDoA method,
+% Performs beamforming and DoA estimation using the diffLinDoA method,
% which is broadband, data-dependent and far-field.
% Y Signal recordings
% doa Directions of Arrival
% pos Positions of sensors
% ndeg Degree of Taylor polynomial
+%
+% Author Clas Veib�ck
c = 343;
[n, N] = size(doa);
diff --git a/Code/linDoA.m b/Code/linDoA.m
index 784db73a099c075437e6f7c98585bddcf81e13de..78e38e02947589165bce3af0a16f13ab10ab5564 100644
--- a/Code/linDoA.m
+++ b/Code/linDoA.m
@@ -1,11 +1,13 @@
function [S, MSE, rho2, xc,epss] = linDoA(Y, doa, pos, ndeg)
%LINDOA Broadband far-field LinDoA beamformer and DoA estimator
-% Performs DoA estimation using the LinDoA method,
+% Performs beamforming and DoA estimation using the LinDoA method,
% which is broadband, data-dependent and far-field.
% Y Signal recordings
% doa Directions of Arrival
% pos Positions of sensors
% ndeg Degree of Taylor polynomial
+%
+% Author Clas Veib�ck
c = 343;
diff --git a/Code/mccc.m b/Code/mccc.m
index 806130e13f031b9e70c68e34cd4b0a7eb3324862..f55ee4b6c3820042be658a72a5ee0f7323c6b93f 100644
--- a/Code/mccc.m
+++ b/Code/mccc.m
@@ -6,6 +6,8 @@ function [R, P] = mccc(Y, doa, pos, fs)
% doa Directions of Arrival
% pos Positions of sensors
% cost Cost function to employ: 'MCCC'(default), 'SLPM', 'BMUSIC'
+%
+% Author Clas Veib�ck
[N, n] = size(Y);
c = 343;
@@ -32,15 +34,5 @@ function [R, P] = mccc(Y, doa, pos, fs)
R(:,:,i) = Ss'*Ss;
P(i) = -logDet(R(:,:,i), 'chol');
end
-
-% S = zeros(N,n);
-% for i=1:numel(doa)
-%
-% for j=1:n
-% S(:,j) = F{j}(ds(j,i)+(1:N)');
-% end
-% R(:,:,i) = S'*S;
-%
-% P(i) = -logDet(R(:,:,i), 'chol');
-% end
+
end
diff --git a/Code/music.m b/Code/music.m
index 5020e6dffcef54a2a63bcbf375216b956399bcf0..1c23fe0ae3e4c10e3a67fe524ab0fffa40a85000 100644
--- a/Code/music.m
+++ b/Code/music.m
@@ -7,6 +7,8 @@ function P = music(Y, doa, pos, f, nsignal)
% pos Positions of sensors
% f Frequency of signal
% nsignal Number of sources (optional)
+%
+% Author Clas Veib�ck
n = size(Y, 2);
c = 343;
diff --git a/Code/mvdr.m b/Code/mvdr.m
index bd971e2ccda45d80ddb2a21326302a8b976387b3..306dd5428f4056fc17588135f1ccfb93b3b7a40a 100644
--- a/Code/mvdr.m
+++ b/Code/mvdr.m
@@ -6,6 +6,8 @@ function [S, P] = mvdr(Y, doa, pos, fs)
% doa Directions of Arrival
% pos Positions of sensors
% f Frequency of signal
+%
+% Author Clas Veib�ck
[N, n] = size(Y);
c = 343;
diff --git a/Code/ncc.m b/Code/ncc.m
index e74fcd22c7f32c8d394ffbc45f2b608c34222600..b40476451a7a6a5e7cf6a446d8e3a5c48b841b09 100644
--- a/Code/ncc.m
+++ b/Code/ncc.m
@@ -1,11 +1,13 @@
function [S, P, lags, quality] = ncc(Y, doa, pos, fs)
%NCC Broadband far-field beamformer (NCC)
-% Performs DoA estimation and beamforming using the NCC method, which is
+% Performs beamforming and DoA estimation using the NCC method, which is
% broadband, data-dependent and far-field.
% Y Signal recordings
% doa Directions of Arrival
% pos Positions of sensors
% fs Sampling frequency of signal
+%
+% Author Clas Veib�ck
n = size(Y, 2);
N = size(Y, 1);
diff --git a/Code/stationaryRTS.m b/Code/stationaryRTS.m
index b18acfc1857030b8998ae21e53e3d7cffa7957b1..2a02f6400c021b0ec533f267635fafbc5d87de7e 100644
--- a/Code/stationaryRTS.m
+++ b/Code/stationaryRTS.m
@@ -1,4 +1,15 @@
function [ xn, Pn ] = stationaryRTS(Y, H, R, F, Q)
+%STATIONARYRTS Stationary Rauch-Tung-Striebel smoother
+% Computes the smoothing estimate using a stationary Rauch-Tung-Striebel
+% smoother. The covariances and gains are assumed stationary and are
+% pre-computed for improved computational performance.
+% Y Signal recordings
+% H Observation matrix
+% R Observation noise covariance
+% F Transition matrix
+% Q Process noise covariance
+%
+% Author Clas Veib�ck
[m, N] = size (Y);
[m1, n] = size(H);
@@ -26,6 +37,7 @@ for i = 1:20
erro = errn;
end
Bf = eye(size(F))-Kf*H;
+A = Bf * F;
KY = Kf*Y;
Ks = Pff * F' / Ppf;
Pn = Pff;
@@ -38,21 +50,25 @@ for i = 1:20
erro = errn;
end
-x = zeros(n,N); % x_{k|k}
-xp = zeros(n,N); % x_{k|k-1}
-xn = x;
+x = ltitr(A,eye(2),[KY(:,2:end) [0;0]]',KY(:,1))';
-x(:,1) = Kf * Y(:,1);
-for k = 2:N
- xp(:,k) = F * x(:,k-1);
- x(:,k) = Bf * xp(:,k) + KY(:,k);
-end
-
-xn(:,N) = x(:,N);
+xf = [fliplr(x(:,1:end-1)) zeros(n,1)];
+u = xf - (Ks * F) * xf;
+xn = fliplr(ltitr(Ks,eye(2),u',x(:,N))');
-for k = N:-1:2
- xn(:,k-1) = x(:,k-1) + Ks * (xn(:,k) - xp(:,k));
-end
+% Actual implementation
+% x = zeros(n,N); % x_{k|k}
+% xp = zeros(n,N); % x_{k|k-1}
+% xn = x;
+% x(:,1) = KY(:,1);
+% for k = 2:N
+% x(:,k) = A * x(:,k-1) + KY(:,k);
+% end
+% xp = F * x;
+% xn(:,N) = x(:,N);
+% for k = N:-1:2
+% xn(:,k-1) = x(:,k-1) + Ks * (xn(:,k) - xp(:,k-1));
+% end
end
diff --git a/Code/tcLinDoA.m b/Code/tcLinDoA.m
index 2248f3ccd0836888b409308018c445ba396bc026..06d6cd25c244477d9d0dd282a9613863b7a56ae3 100644
--- a/Code/tcLinDoA.m
+++ b/Code/tcLinDoA.m
@@ -1,13 +1,15 @@
function [S, MSE, kurtosis, xc] = tcLinDoA(Y,doa,pos, ndeg, fs, directTime)
%TCLINDOA Broadband far-field TCLinDoA beamformer and DoA estimator
-% Performs DoA estimation using the time-constrained LinDoA method,
-% which is broadband, data-dependent and far-field.
+% Performs beamforming and DoA estimation using the time-constrained LinDoA
+% method, which is broadband, data-dependent and far-field.
% Y Signal recordings
% doa Directions of Arrival
% pos Positions of sensors
% ndeg Degree of Taylor polynomial
% fs Sampling rate
% directTime Flags that doa contains time delays rather than directions
+%
+% Author Clas Veib�ck
c = 343;
@@ -44,9 +46,7 @@ function [S, MSE, kurtosis, xc] = tcLinDoA(Y,doa,pos, ndeg, fs, directTime)
H = T;
R = eye(m);
F = kron( eye(n),[1 1/fs; 0 1]);
- Q = kron(eye(n),[0 0; 0 1e9]); % Assumes second order
- % F = kron( eye(length(Th)),[0 0; 0 0]);
- % Q = kron(eye(n),[1e9 0; 0 1e9]); % Assumes second order
+ Q = kron(eye(n),[0 0; 0 1]); % Assumes second order
[St, ~] = stationaryRTS(Y', H, R, F, Q);
S(:,:,i) = St';
diff --git a/README.md b/README.md
index da0a8f0955f2fd69a136de05a9bb05cdf1c8117c..d032d3b7640d8471dc7be32f09c10a85524155e5 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@ All methods work for single sources, and some also work for multiple sources. Pr
The following algorithms are implemented
- LinDoA Beamformer and Doa Estimator
-- Differentiated LinDoA (Diff-LinDoA) Beamformer and Doa Estimator
+- Differential LinDoA (Diff-LinDoA) Beamformer and Doa Estimator
- Time-constrained LinDoA (TC-LinDoA) Beamformer and Doa Estimator
- Delay and Sum (DaS) Beamformer and DoA estimator
- Minimum Variance Distortionless Response (MVDR) Beamformer and Doa Estimator