clear;

[p, fs] = audioread('ket_zaj_all_001_direct.wav');
M = size(p, 2);

N = size(p,1);
t = (0 : N-1)' / fs;

%% mixdown
f0 = 1000;
om0 = 2*pi*f0;
mx = exp(-1i*om0*t);
pmx = bsxfun(@times, mx, p);

%%
fc = 5;
[b, a] = butter(3, fc/(fs/2));
P = filter(b, a, pmx);

%% R - CSM
% idx = [36 : -1 : 25, 1 : 12 ];%48 : -1 : 37, 13 : 24];
R = zeros(M, M);
tlim = [3 4];
idx = round(tlim * fs);
for i = idx(1) : idx(2)
    R = R + P(i,:).' * conj(P(i,:));
end
R = R / (idx(2)-idx(1) + 1);

