% === ryx_RIR_corrected2_stereo.m === 

% Дана програма забезпечує: 1) уточнене знаходження піків ВКФ з метою
% отримання якісної усередненої ВКФ; таке уточнення відбувається пошуком
% максимумів піків 
% 2) корекція форми усередненого піку ВКФ врахуванням АЧХ гучномовця
% 3) при корекції придушуються високі частоти
% 4) обчислюється 2 канали

RIRs = zeros(65536,2);
% ВКФ як оцінка імп.хар-ки гучномовця
filename_x = input('Вкажіть ім"я файлу зразка (без лапок):   ', 's');
[x,Fs]=audioread(filename_x);
%x = [0; x];
filename_y = input('Вкажіть ім"я файлу запису (без лапок):   ', 's');
[y,fs]=audioread(filename_y);
%y = [0; y];
for chan = 1:2
    ryx_rir=xcorr(y(:,chan),x);
%ryx_rir_logabs=20*log10(abs(ryx_rir/max(abs(ryx_rir))));

% побудова графіків
ty=1/fs:1/fs:length(y)/fs;
%figure;plot(ty,y);
%xlabel('Time (s)'); ylabel('Level'); title('Signal'); grid on;
tx=1/fs:1/fs:length(x)/fs;
%figure;plot(tx,x);
%xlabel('Time (s)'); ylabel('Level'); title('Reference'); grid on;
%tr1=1/fs:1/fs:length(ryx_rir_logabs)/fs;
%figure;plot(tr1,ryx_rir_logabs);
%xlabel('Time (s)'); ylabel('Level'); title('Ryx log abs'); grid on;

% === Нарізаємо ВКФ на сегменти ===
% шукаємо перший сплеск
max_rir = max(abs(ryx_rir));
tresh = 0.7*max_rir;
k = find(diff(sign(abs(ryx_rir)-tresh)),1);    % номер відліку першого сплеску

%% === Початок уточнення знаходження піків ===
% нарізаємо ВКФ на сегменти - стовпці матриці
%n_pik = input('Вкажіть кількість піків ВКФ:   ');
n_pik = 17;
k1=k-20;                          % початок зони пошуку 1-го піка
k41 = k+20;                 % кінець зони із 41 точок
[k1max,i1max] = max(abs(ryx_rir(k1:k41)));   % максимум 1-го піка та його індекс
% пошук максимумів піків
mm_rir = zeros(41,n_pik);
for nn = 1:n_pik
    mm_rir_frag = ryx_rir((k1+(nn-1)*65536):(k41+(nn-1)*65536));
    mm_rir(:,nn) = mm_rir_frag;
end
[k1maxm,i1maxm] = max(abs(mm_rir));   % максимуми піків та їх індекси

m_rir = zeros(65536,n_pik);
for nn = 1:n_pik
    m_rir_frag = ryx_rir((k+i1maxm(nn)-24+(nn-1)*65536):(k+i1maxm(nn)-25+nn*65536));
    m_rir(:,nn) = m_rir_frag;
end

%% === Кінець уточнення знаходження піків ===

% перевірка правильності вибору n-го піка
tpik = 1/fs:1/fs:65536/fs;
for nn = 1:1
    %figure; plot(tpik(1:200),m_rir(1:200,nn));
   % title(['CCR peak No.' num2str(nn)]);
end
% перевірка результату усереднення піків ВКФ
m_rir(:,1) = [];
mean_rir = mean(m_rir,2);
%figure; plot(tpik(1:200),mean_rir(1:200,1));
%title('Averaged CCR');

%% корекція форми усередненого піка ВКФ
fCCR = fft(mean_rir, 65536);
H = load('fresp_Genius_delf50Hz_Type3.mat');
%H = load('b_proba.mat');
fresp1 = H.fresp;
if size(fresp1,2) == 1;
    fresp = fresp1;
else
    fresp = fresp1(:,chan);
end
fresp_ = flipud(fresp(2:end-1));
full_fresp = [fresp; fresp_];
% НЧ-фільтр - спроба 1
f1 = 1;          % гранична частота НЧ фільтру для придушення ВЧ
Nfft = 65536; Fs = 44100;
k = round((Fs/2-f1)/(Fs/Nfft)); n0 = Nfft/(2*k);
fragcos = ones(Nfft,1);
for r = Nfft/2+1-k:Nfft/2+1+k
    fragcos(r,1) = 0.5*(1-cos(2*pi/Nfft*(r-Nfft/2-1)*n0));
end
% НЧ-фільтр - спроба 2
f2 = 16000;                      % гранична частота НЧ фільтру для придушення ВЧ
dfhn = round(2*f2/(Fs/Nfft));       % ширина вікна Ханна у вибірках
if mod(dfhn,2)==0
    dfhn = dfhn-1;
end;
hn = hann(dfhn);
fragcos2 = zeros(Nfft,1);
fragcos2(1:dfhn) = hn(1:dfhn);
fragcos2 = circshift(fragcos2, -(dfhn-1)/2);

%RIR = real(ifft(fCCR./full_fresp.*fragcos));
RIR = real(ifft(fCCR./full_fresp.*fragcos2));
RIR = RIR/max(abs(RIR));

RIRs(:,chan) = -RIR(:,1);
% нейтралізація фазових спотворень гучномовця
%ph = load('phresp_Genius', 'phresp');
%phz = smooth(unwrap(ph.phresp));
%RIR = real(ifft(fCCR./full_fresp.*fragcos.*exp(-j*phz)));

end

%% графіки піків 
figure;plot(tpik(1:fs/10),RIRs(1:fs/10,1)); 
title('Left');
xlabel('Time (s)'); ylabel('Level')
figure; plot(tpik(1:fs/10),RIRs(1:fs/10,2));
title( 'Right');
xlabel('Time (s)'); ylabel('Level')

%  запис RIRs на диск як wav-файл
filename_RIRs = input('Вкажіть ім"я файлу ІХ (без лапок):   ', 's');
audiowrite(filename_RIRs,RIRs,fs,'BitsPerSample',24);