Micdelay

meanTp = round(mean(diff(rawTimes)))*1e-6;
fs = 1/meanTp(1);
%insert samplerange you want to crosscorrelate
a = 2000;
b = 20000;

%Removing DC components BEFORE crosscorrelation
m1 = allData.Mic1(a:b) -512;
m2 = allData.Mic2(a:b) -512;
m3 = allData.Mic3(a:b) -512;

%Crosscorrealtion between all three mics
[acor1, lag1] = xcorr(m1, m2);
[~, I_1] = max( abs(acor1));
lagDiff1 = lag1(I_1);
timeDiff1 = lagDiff1/fs;

[acor2, lag2] = xcorr(m1, m3);
[~, I2] = max( abs(acor2));
lagDiff2 = lag2(I2);
timeDiff2 = lagDiff2/fs;

[acor3, lag3] = xcorr(m2, m3);
[~, I3] = max( abs(acor3));
lagDiff3 = lag3(I3);
timeDiff3 = lagDiff3/fs;

%Incoming angle estimation
theta = atan(sqrt(3)*((timeDiff1 + timeDiff2)/(timeDiff1 - timeDiff2 - 2*timeDiff3)));
if (lagDiff3 < 0),
    theta = theta + pi;
elseif (lagDiff1 == lagDiff2 && lagDiff1 > 0),
        theta = -pi/2;
elseif (lagDiff1 == lagDiff2 && lagDiff1 < 0),
        theta = pi/2;
end

%plot relevant figures
figure;
subplot(3,1,1);
plot(lag1/fs, acor1);
title('Lag from m2 to m1');
xlabel(['sample difference: ' num2str(lagDiff1) ' Time difference: ' num2str(timeDiff1)]);

subplot(3,1,2);
plot(lag2/fs, acor2);
title('Lag from m3 to m1');
xlabel(['sample difference: ' num2str(lagDiff2) ' Time difference: ' num2str(timeDiff2)]);

subplot(3,1,3);
plot(lag3/fs, acor3);
title('Lag from m3 to m2');
xlabel(['sample difference: ' num2str(lagDiff3) ' Time difference: ' num2str(timeDiff3)]);

figure;
plot(allTimes.Mic1(a:b),m1, '-o',...
    allTimes.Mic2(a:b),m2, '-o',...
    allTimes.Mic3(a:b),m3, '-o'...
    );
title('Part of signals m1, m2 and m3 that are crosscorrelated');
xlabel('t [s]');
ylabel('10-bit Conversion value');
legend('Mic1','Mic2','Mic3');

figure;

plot([0 sqrt(3)/2],[1 -0.5], 'g--o',...
    [-sqrt(3)/2 0], [-0.5 1], 'b--o',...
    [sqrt(3)/2 -sqrt(3)/2], [-0.5 -0.5], 'r--o',...
    [0 10*cos(theta)], [0 10*sin(theta)]...
    );
axis([-1.5 1.5 -1.5 1.5]);
text(0, 1.2, 'm1');
text(-1, -0.7, 'm2');
text(1, -0.7, 'm3');
grid on;