
% roomAdapter.m 
%
% Dylan Menzies 2017-2020
% Copyright University of Southampton 
% 
% Modify target reverb response to generate a play back response,
% that produces the target approximately when played into a room
% 'room in room' effect
%
%


% clearall;
close all;
load('rng1');

path(pathdef);

rng_save = rng;
rng(rng1);

%%


fs = 44100;

% Time duration of early response, s
tEarly = 0.020; %0.020 0.009; .,.
NI_ = floor(tEarly * fs);
NI = NI_;   % inlcudes added zeros
% NI = NI_* 2;  

nSamplesPlot = 10000;

Nf = 10; % order of filters (even number)
NbandsPerOctave = 3;

[filterBank, fmb] = filterBank_create(Nf, NbandsPerOctave, fs);   % fmb are the band centre frequencies
Nb = length(filterBank); % number of bands, set in filter_bank

% Check behaviour of filterBank, and find compensation factors
% [IR, fs_] = audioread('35x25y.wav');
% IR = IR(1:NI_)';
IR = [1];
IRb = filterBank_apply(IR, filterBank);
DIRb = calc_densities(IRb, IRb, fmb, fs) % (+) 
bandDensityOverlapFactor = 1.013; % Estimte from (+)
IR_ = sum(IRb);

IRb_ = filterBank_apply(IR_, filterBank);
DIRb_ = calc_densities(IRb_, IRb_, fmb, fs); 
DIRb_ ./ DIRb   % (*)
bandSumOverlapFactor = sqrt(1.136); % Estimate from (*)

% figure(1)
% plot([1:length(IR)]/fs, IR)
% figure(2)
% plot(abs(fft(IR)))

IRenergy = IR * IR' % expect some losses - bands don't cover full band


%% Test band density multiplication
% Db=0;
% n=0;
% %%
% % rng(.1);
% Nsamples = 100000;
% It = rand(1,Nsamples)-0.5;
% Ir = rand(1,Nsamples)-0.5;
% Itb = filterBank_apply(It, filterBank);
% Irb = filterBank_apply(Ir, filterBank);
% ItIrb = filterBank_apply(convfft(It, Ir), filterBank);
% It_Irb = filterBank_apply((It + Ir), filterBank);
% DItb = calc_densities(Itb, Itb, fmb, fs) / bandDensityOverlapFactor;
% DIrb = calc_densities(Irb, Irb, fmb, fs) / bandDensityOverlapFactor;
% DItIrb = calc_densities(ItIrb, ItIrb, fmb, fs) / bandDensityOverlapFactor;
% DIt_Irb = calc_densities(It_Irb, It_Irb, fmb, fs) / bandDensityOverlapFactor;
% % Db = Db + DItb + DIrb;
% % n = n+1;
% % Db/n
% % figure
% % plot(DItb)
% % plot(DItIrb ./ DItb ./ DIrb)
% plot(DIt_Irb ./ (DItb + DIrb) )
% % plot(DIt_Irb ./ 1.66e+04 )
% hold on
%%

RT60offset = 5000; % For RT60 : offset to start point for analysis in samples.

maxgain = 20; % Limit equalisation gains.



%%

% TARGET REVERB, I early, V late

% Synthesized target
It = [1]; 

% Band levels .,.
b = ones(1, Nb) /36;      %38 30 24 10      %[1.2:-1/(Nb-1):0.2]; %0.2*ones(1, Nb);  %[1:-1/(Nb-1):0];

% Band decays
beta = 3*[1:3/(Nb-1):4]; % 3*[1:3/(Nb-1):4]; %4*[1:19/(Nb-1):20]; %4*[1:3/(Nb-1):4]; %2*[1:7/(Nb-1):8]; % 4* 8*  [1:3/(Nb-1):4] %4.0*ones(1, Nb); %[1:9/(Nb-1):10];   4.0*ones(1, Nb);
beta_t = beta';

tPre = 0.020;
Npre = floor(tPre * fs);
Vt = [ zeros(1,Npre), make_reverb(beta, b, 2.0, fs, filterBank, fmb) ];

%plot(abs(fft(Vt))); %check how band synthesis links together

IVt = Vt + [It, zeros(1, length(Vt)-length(It))];


%% Recorded target - comment out if synthesized
target_recorded = 1;
%addpath('../../rooms/4. AIR/staircase/BRIRs/')
addpath('rooms')
addpath('rooms/C4DM_classroom')
 [IVt, fs_] = audioread('0.wav');
% [IVt, fs_] = audioread('alcuin.wav');
% [IVt, fs_] = audioread('basement.wav');
% [IVt, fs_] = audioread('stairway.wav');
% [IVt, fs_] = audioread('25x15y.wav');
% [IVt, fs_] = audioread('35x25y.wav');
% addpath('../../rooms/2. BBCusability/RIRs/pressure/low5.0/')
% [IVt, fs_] = audioread('C.wav');
%IVt = IVt(floor(0.010 * fs):end, 1)';  % remove lead space
IVt = IVt(:,1)';
Vt = [ zeros(1,NI), IVt(NI_+1:length(IVt)) ];
It = IVt(1:NI_);
%%% Modify gain of late reponse
Vt = Vt * 2;
IVt = [ It, Vt(NI+1:length(Vt)) ];



% figure
% plot([1:nSamplesPlot]/fs, IVt(1:nSamplesPlot))
% title('Target Impulse Response')
% xlabel('Time (s)')




%%

% ROOM REVERB

% Synthesized room
Ir = [1]; 
% bigger -> easier to compensate, but less interesting. 
% smaller -> compensation eventually breaks down, tradeoff required.

% (solvable for room with lower level and higher decay rate than target)

% Band levels .,.
b = ones(1, Nb) /28; %1/31 = 0.032  1/28=0.036  %31 28      %30 24 20  15 no break down  %12 small breakdown     %[1.2:-1/(Nb-1):0.2];
%b([1:10]) = 1.5*b([1:10]);    % boost bass to make more noticeable
%b(Nb) = 0;
%b(6:8) = 1.3;

%b([Nb-13:Nb]) = 0.5*ones(1,14);

%b([1:10]) = 0.5*ones(1,10);


% Band decays .,.
beta = 4*[1:5/(Nb-1):6]; %4* %4*[1:5/(Nb-1):6]; % %8*[1:31/(Nb-1):32]; %10.0*ones(1, Nb);  % 4+20*[0:Nb-1]/(Nb-1);  %10.0*ones(1, Nb); 
% beta([Nb-5:Nb]) = 200*ones(1,6);
% beta(Nb) = 200;
beta_r = beta';

tPre = 0.005; %0.005;
Npre = floor(tPre * fs);
Vr = [ zeros(1,Npre), make_reverb(beta, b, 2.0, fs, filterBank, fmb) ];
% Vr = [0 0 0 0];
IVr = Vr + [Ir, zeros(1, length(Vr)-length(Ir))];


%% Recorded room - comment out if synthesized
room_recorded = 1;
%addpath('../../rooms/2. BBCusability/RIRs/pressure/low5.0/')
addpath('rooms')
 [IVr, fs_] = audioread('C.wav');
% [IVr, fs_] = audioread('genesis6.wav');
%[IVr, fs_] = audioread('basement.wav');
IVr = IVr(:, 1)';
%IVr = IVt; %%%
Vr = [ zeros(1,NI), IVr(NI_+1:length(IVr)) ];
Ir = IVr(1:NI_);

% Modify gain of late response
Vr = Vr *2; %*0.01;% / 2.5; %.,.
IVr = [ Ir, Vr(NI+1:length(Vr)) ];

figure
plot([1:nSamplesPlot]/fs, IVr(1:nSamplesPlot))
title('Room Impulse Response')
xlabel('Time (s)')

%%


% Decompose into signal bands
IVtb = filterBank_apply(IVt, filterBank);
Itb = filterBank_apply(It, filterBank);
Vtb = filterBank_apply(Vt, filterBank);

IVrb = filterBank_apply(IVr, filterBank);
Irb = filterBank_apply(Ir, filterBank);
Vrb = filterBank_apply(Vr, filterBank);


ItIrb = filterBank_apply(convfft(It, Ir), filterBank);
VtIVrb = filterBank_apply(convfft(Vt, IVr), filterBank);
ItVrb = filterBank_apply(convfft(It, Vr), filterBank);

VtIrb = filterBank_apply(convfft(Vt, Ir), filterBank);
VtVrb = filterBank_apply(convfft(Vt, Vr), filterBank);


% Calculate energy spectral densities in bands
DIVtb = calc_densities(IVtb, IVtb, fmb, fs) / bandDensityOverlapFactor;
DItb = calc_densities(Itb, Itb, fmb, fs) / bandDensityOverlapFactor;
DVtb = calc_densities(Vtb, Vtb, fmb, fs) / bandDensityOverlapFactor;

DIVtb ./ (DItb + DVtb);  % check energy conserved by split
% Low freqencies inaccurate because early window too short

DIVrb = calc_densities(IVrb, IVrb, fmb, fs) / bandDensityOverlapFactor;
DIrb = calc_densities(Irb, Irb, fmb, fs) / bandDensityOverlapFactor;
DVrb = calc_densities(Vrb, Vrb, fmb, fs) / bandDensityOverlapFactor;

DIVrb ./ (DIrb + DVrb); %??

% Clarity ratio prediction calc. based on initial parameters
Rtb = DVtb ./ DItb;
Rrb = DVrb ./ DIrb;
ERtr_tb = 1 ./ ( 1 + Rrb + Rrb ./ Rtb ); % expected ratio tr / t


%%

%%% Uncomment for recorded impulse responses
if exist('target_recorded')
% RT60offset = length(Vt)/2
for n = 1:Nb
    beta_t(n) = 3000/t60(Vtb(n,:), fs, 0, RT60offset) * log(10);
end
end

%%% Uncomment for recorded impulse responses
if exist('room_recorded')
% RT60offset = length(Vrb(1,:))/2
for n = 1:Nb
    beta_r(n) = 3000/t60(Vrb(n,:), fs, 0, RT60offset) * log(10);
end
end


% For precise cross-density calculations
DItIrb = calc_densities(ItIrb, ItIrb, fmb, fs) / bandDensityOverlapFactor;
DVtIVrb = calc_densities(VtIVrb, VtIVrb, fmb, fs) / bandDensityOverlapFactor;
DItVrb = calc_densities(ItVrb, ItVrb, fmb, fs) / bandDensityOverlapFactor;
DVtIVr_ItVrb = calc_densities(VtIVrb, ItVrb, fmb, fs) / bandDensityOverlapFactor;

DVtIrb = calc_densities(VtIrb, VtIrb, fmb, fs) / bandDensityOverlapFactor;
DVtVrb = calc_densities(VtVrb, VtVrb, fmb, fs) / bandDensityOverlapFactor;



% Check statistical assumption for band calculation .,.
% DItIrb_ = DItIrb ./ DItb ./ DIrb
% DVtIVrb_ = DVtIVrb ./ DVtb ./ DIVrb
% DItVrb_ = DItVrb ./ DItb ./ DVrb
% DVtIrb_ = DVtIrb ./ DVtb ./ DIrb
% DVtVrb_ = DVtVrb ./ DVtb ./ DVrb

%%% test bands disjoint
% Vr * Vr'    
% Vr_ = sum(Vrb, 1);
% Vr_ * Vr_'
% 
% Vtr = convfft(Vt, Vr);  % test 
% Vtrb = filter_bank(Vtr, fs);
% DVtrb = calc_densities(Vtrb, fmb, fs);
% DVtb .* DVrb  ./  DVtrb
% s



%%

% EQUALISATION


% DIrb(DIrb == 0) = 1; % prevent divsion by 0
% DIrb(DIrb < 0.1) = 0.1; % limit EQ


% Calculate late and early playback reverb band densities

% Initial guess at direct playback:
DIpb = DItb ./ DIrb;

% Calculate resulting fraction of reproduced late energy due to early playback - late room:
DVprbInDVtb = DIpb .* DVrb ./ DVtb;

% Tradeoffs if exact solution not possible. 
%.,.

%Reduce early playback to prevent excess late reproduced energy, (at
%expense of early reproduced):
gamma0 = 1.0 %1.0 0.8 1  % Max limit of IpxVr contribution in Vt, for any room decay.

i = find( DVprbInDVtb > gamma0 );
DIpb(i) = gamma0 * (DVtb(i)) ./ (DVrb(i));

% 
% 
% % Reduce early playback to promote contribution from late playback to late reproduced, when
% % late target decay time is significantly longer than late room decay time (overwrites excess prevention above)
% %
beta_rt_min = 1.2;    % Min amount by which room decay time is shorter than target rate for alpha limit to apply,
                      % to reduce influence of late room in late reproduced.
gamma1 = 0.9 %1.0 0.8 %0.6 0.7 0.5 1.0    % Max limit of IpxVr contribution in Vt, when room decay is fast.

i = find( beta_r ./ beta_t > beta_rt_min  &  DVprbInDVtb > gamma1 );
DIpb(i) = gamma1 * (DVtb(i)) ./ (DVrb(i));



% 
% 
% 
% 
% eta = 0.0; % Proportion of missing early reproduced energy to be produced using late energy
%? prob better to normalise overall energy in band - 
%%%
% DVtb_ = DVtb  + eta * (DItb - DIpb .* DIrb);       % Transfer any energy missing from early playback - early room to late target.
DVpb = max(0,  DVtb - DIpb .* DVrb ) ./ DIVrb;    % Calculate late playback
% DVpb = max(0,  DVtb - DIpb .* DVrb ) ./ (DIrb + DVrb);    % Calculate late playback .,.
% 
% % Prediction of total reproduced energy density: should be near 1.0 in each
% ( DVpb .* DIVrb  +  DIpb .* DIVrb ) ./ DIVtb
% 
% Check predicted early energy
% ( DIpb .* DIrb  ) ./ DItb
% 
% Check predicted late energy
% ( DVpb .* DIVrb  +  DIpb .* DVrb  ) ./ DVtb
%
% 
% beta_rt_max = 0.5;
% % Zero playback for bands where the room decays significantly slower than
% % target
% i = find( beta_r ./ beta_t < beta_rt_max );
% % i = find( beta_r ./ beta_t < beta_rt_max &  DVrbInDVtb(i) > 1);
% DIpb(i) = 0;
% %DVpb(i) = 0;
% 
% 
% DIpb(i) = DIVtb(i) ./ ( DItb(i) .* DIVrb(i) );


%% Normalise predicted total reproduced energy to target
f = DIVtb ./ ( (DIpb + DVpb) .* DIVrb ); 
% f = (DItb + DVtb) ./ ( (DIpb + DVpb) .* (DIrb + DVrb) ); 
DIpb = DIpb .* f;
DVpb = DVpb .* f;


%%%DVpb =  DVtb ./ DVrb;    % simple late only
% ( (DIpb .* DVrb) + DVpb .* (DIrb + DVrb) ) ./ DVtb   % DVprb / DVtb




% Clarity ratio prediction calc. based on initial parameters
Rpb = DVpb ./ DIpb;
ERpr_tb = Rtb ./ ( Rpb .* (1 + Rrb) + Rrb );




%% Object-based equalisation

% Early eqalisation
%Ipb = Itb ./ sqrt(DIrb);
Ipb = Itb .* min(maxgain, sqrt( abs(DIpb ./ DItb) ) );
Ip = sum(Ipb, 1) / bandSumOverlapFactor;  % total impulse = sum across bands  %! with correction for band overlap

% Late eqalisation
Vpb = Vtb .* min(maxgain, sqrt( abs(DVpb ./ DVtb) ) );
Vp = sum(Vpb, 1) / bandSumOverlapFactor;    
%Vp = Vp([Npre+1:length(Vp)]);  %%% Remove leading space, NI samples, to offset convolution delay
%Vp = [ zeros(1, 1000), Vp ];    % Reduce padding to compensate for swell delay caused by Vr 

IVp = Vp + [Ip, zeros(1, length(Vp)-length(Ip))];




%% Precision equalisation, based on cross-density

% gIb = sqrt(DItb ./ DItIrb);
% 
% % gIb = 2*ones(23,1);
% % gIb = [1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1]';
% 
% a = DVtIVrb;
% b = 2* gIb .* DVtIVr_ItVrb;
% c = gIb .^2 .* DItVrb - DVtb;
% 
% gVb = real(  ( -b + sqrt(b.^2 - 4*a.*c) )  ./ (2*a)  );
% 
% % check
% %( DItb ./ DItIrb  ) ./ (gIb .^2) 
% %( gVb .^2 .* DVtIVrb + 2*gVb.*gIb .* DVtIVr_ItVrb + gIb .^2 .* DItVrb ) ./ DVtb
% 
% Ip = (gIb' * Itb) / bandSumOverlapFactor;
% Vp = (gVb' * Vtb) / bandSumOverlapFactor;
% 
% IVp = [Ip, zeros(1, length(Vp)-length(Ip))] + Vp;


%% Channel equalisation 0 - Simple overall equalisation.

% Separate I, V soln should improve on this!!
% IVp0b = IVtb ./ sqrt( abs(DIVrb) ); % Compensate for eq of whole room
% IVp0 = sum(IVp0b, 1) / bandSumOverlapFactor;
% Ip0 = IVp0(1:NI);   %! cuts off bass response
% Vp0 = [zeros(1,NI), IVp0([NI+1:length(IVp0)]) ];

Ip0b = Itb .* min(maxgain, 1./sqrt(DIVrb)); % Compensate for eq of whole room
Ip0 = sum(Ip0b, 1) / bandSumOverlapFactor;
Vp0b = Vtb .* min(maxgain, 1./sqrt(DIVrb)); % Compensate for eq of whole room
Vp0 = sum(Vp0b, 1) / bandSumOverlapFactor;
IVp0 = Vp0 + [Ip0, zeros(1, length(Vp0)-length(Ip0))];

% IVp0b = IVtb .* min(maxgain, 1./sqrt(DIVrb)); % Compensate for eq of whole room
% IVp0 = sum(IVp0b, 1) / bandSumOverlapFactor;


%% Channel equalisation 1 - Equalise so that early response is flat

% IVp1b = IVtb ./ sqrt( abs(DIrb) ); % Compensate for eq of early room
% IVp1 = sum(IVp1b, 1) / bandSumOverlapFactor;
% Ip1 = IVp1(1:NI);   %! cuts off bass response
% Vp1 = [zeros(1,NI), IVp1([NI+1:length(IVp1)]) ];

% Ip1b = Itb .* min(maxgain, 1./sqrt(DIrb)); % Compensate for eq of early room
% Ip1 = sum(Ip1b, 1) / bandSumOverlapFactor;
% Vp1b = Vtb .* min(maxgain, 1./sqrt(DIrb)); % Compensate for eq of early room
% Vp1 = sum(Vp1b, 1) / bandSumOverlapFactor;
% IVp1 = Vp1 + [Ip1, zeros(1, length(Vp1)-length(Ip1))];

% IVp1b = IVtb .* min(maxgain, 1./sqrt(DIVrb)); % Compensate for eq of early room
% IVp1 = sum(IVp1b, 1) / bandSumOverlapFactor;

%%%
%%% Uncomment for density analysis
%%%
% IVpb = IVp0b;
% IVp = IVp0;
% Ip = Ip0;

% Test filter analysis and recomposition
% IVpb = filterBank_apply(IVp, filterBank);
% DIVpb = calc_densities(IVpb, IVpb, fmb, fs);

% figure
% plot(IVp(1:nSamplesPlot))
% title('playback')


%% Convolve target and playback with room

% IVtr = convfft(IVt, IVr);
% Itr = convfft(It, Ir);
% Vtr = IVtr - [Itr, zeros(1, length(IVtr) - length(Itr))];
%Vtr = convfft(Vt, Vr);
% Itr = IVtr(1:NI);
% Vtr = [ zeros(1,NI), IVtr(NI+1:length(IVtr)) ];

IVpr = convfft(IVp, IVr);
Ipr = convfft(Ip, Ir);  %! Note at low frequency Ipr will continue past NI window due to band EQ in Ip: Need better EQ method.
Vpr = IVpr - [Ipr, zeros(1, length(IVpr) - length(Ipr))];
% Ipr = IVpr(1:NI);
% Vpr = [ zeros(1,NI), IVpr(NI+1:length(IVpr)) ];

IVp0r = convfft(IVp0, IVr); % Channel
Ip0r = convfft(Ip0, Ir); % Channel
Vp0r = IVp0r - [Ip0r, zeros(1, length(IVp0r) - length(Ip0r))];
% Ip0r = IVp0r(1:NI);
% Vp0r = [ zeros(1,NI), IVp0r(NI+1:length(IVp0r)) ];

% IVp1r = convfft(IVp1, IVr); % Channel
% Ip1r = convfft(Ip1, Ir); % Channel
% Vp1r = IVp1r - [Ip1r, zeros(1, length(IVp1r) - length(Ip1r))];

% Ip1r = IVp1r(1:NI);
% Vp1r = [ zeros(1,NI), IVp1r(NI+1:length(IVp1r)) ];

% figure
% plot(IVpr(1:nSamplesPlot))
% title('playback into room')


%%
% close all
% n = 30000;
% figure
% plot(IVt(1:n))
% figure
% plot(IVr(1:n))
% figure
% plot(IVtr(1:n))
% s

%% Reform impulse check

Ipb_ = filterBank_apply(Ip, filterBank);    % reanalyse
DIpb_ = calc_densities(Ipb_, Ipb_, fmb, fs) / bandDensityOverlapFactor;

%% Find densities of playback-through-room impulse

IVprb = filterBank_apply(IVpr, filterBank);
Iprb = filterBank_apply(Ipr, filterBank);
Vprb = filterBank_apply(Vpr, filterBank);
IVp0rb = filterBank_apply(IVp0r, filterBank);
Ip0rb = filterBank_apply(Ip0r, filterBank);
Vp0rb = filterBank_apply(Vp0r, filterBank);
% IVp1rb = filterBank_apply(IVp1r, filterBank);
% Ip1rb = filterBank_apply(Ip1r, filterBank);
% Ip1b = filterBank_apply(Ip1, filterBank);
% Vp1rb = filterBank_apply(Vp1r, filterBank);
% IVtrb = filterBank_apply(IVtr, filterBank);
% Itrb = filterBank_apply(Itr, filterBank);
% Vtrb = filterBank_apply(Vtr, filterBank);

DIVprb = calc_densities(IVprb, IVprb, fmb, fs) / bandDensityOverlapFactor;
DIprb = calc_densities(Iprb, Iprb, fmb, fs) / bandDensityOverlapFactor;
DVprb = calc_densities(Vprb, Vprb, fmb, fs) / bandDensityOverlapFactor;
DIVp0rb = calc_densities(IVp0rb, IVp0rb, fmb, fs) / bandDensityOverlapFactor;
DIp0rb = calc_densities(Ip0rb, Ip0rb, fmb, fs) / bandDensityOverlapFactor;
DVp0rb = calc_densities(Vp0rb, Vp0rb, fmb, fs) / bandDensityOverlapFactor;
% DIVp1rb = calc_densities(IVp1rb, IVp1rb, fmb, fs) / bandDensityOverlapFactor;
% DIp1rb = calc_densities(Ip1rb, Ip1rb, fmb, fs) / bandDensityOverlapFactor;
% DIp1b = calc_densities(Ip1b, Ip1b, fmb, fs) / bandDensityOverlapFactor;
% DVp1rb = calc_densities(Vp1rb, Vp1rb, fmb, fs) / bandDensityOverlapFactor;
% DIVtrb = calc_densities(IVtrb, IVtrb, fmb, fs) / bandDensityOverlapFactor;
% DItrb = calc_densities(Itrb, Itrb, fmb, fs) / bandDensityOverlapFactor;
% DVtrb = calc_densities(Vtrb, Vtrb, fmb, fs) / bandDensityOverlapFactor;

% Calculate direct to late energy ratios
% measured from input
Rtb = DItb ./ DVtb;
Rrb = DIrb ./ DVrb;
% calculated
Rpb = DIpb ./ DVpb;
% measured from convolution
Rprb = DIprb ./ DVprb;
Rp0rb = DIp0rb ./ DVp0rb;
% Rp1rb = DIp1rb ./ DVp1rb;
% Rtrb = DItrb ./ DVtrb;

 

% Rpb = DVpb ./ DIpb;
% Rprb = DVprb ./ DIprb;
% Rp0rb = DVp0rb ./ DIp0rb;
% Rrb = DVrb ./ DIrb;
% Rtrb = DVtrb ./ DItrb;
% Rtb = DVtb ./ DItb; 

% Calculate ratio of energy ratios
Rpr_tb = Rprb ./ Rtb;
Rp0r_tb = Rp0rb ./ Rtb;
% Rp1r_tb = Rp1rb ./ Rtb;
% Rtr_tb = Rtrb ./ Rtb;
Rt_rb = Rtb ./ Rrb;
Rr_tb = Rrb ./ Rtb;


% Compare energy densities of target against the combined playback and room
% DIVprb ./ DIVtb
% DIprb ./ DItb
% DVprb ./ DVtb


%%

figure
semilogx(fmb', DIVprb ./ DIVtb, '-o')
xlim([fmb(1) fmb(end)])
ylim([0.0,2.0])
hold on
semilogx(fmb', DIVp0rb ./ DIVtb, '-x')
% semilogx(fmb',DIVp1rb ./ DIVtb, '-+')
%semilogx(fmb',DIVrb ./ DIVtb, 'g-^')
semilogx(fmb', ones(1,Nb), '--', 'color', 'black')
title({'Ratio of total densities'})
% h_legend = legend({'object\_playback-room / target', 'channel\_playback-room / target', 'I\_channel\_playback-room / target'}, 'Location', 'Southeast');
h_legend = legend({'object\_playback-room / target', 'channel\_playback-room / target'}, 'Location', 'Southeast');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')



% figure
% % semilogx(fmb',DIprb, '-o')
% xlim([fmb(1) fmb(end)])
% ylim([0 inf])
% 
% % semilogx(fmb',DIp0rb, '-x')
% % semilogx(fmb',DIp1rb, '-+')
% % semilogx(fmb',DItb, '-v')
% % hold on;
% % semilogx(fmb',DVtb, '-^')
% % 
% semilogx(fmb',DItb0, '-v', 'color', 'blue')
% hold on;
% semilogx(fmb',DVtb0, '-^', 'color', 'red')
% % semilogx(fmb',DItb1, '-v', 'color', 'blue')
% % semilogx(fmb',DVtb1, '-^', 'color', 'red')
% semilogx(fmb',DItb2, '-v', 'color', 'blue')
% hold on
% semilogx(fmb',DVtb2, '-^', 'color', 'red')
% semilogx(fmb',DItb3, '-v', 'color', 'blue')
% semilogx(fmb',DVtb3, '-^', 'color', 'red')
% semilogx(fmb',DItb4, '-v', 'color', 'blue')
% semilogx(fmb',DVtb4, '-^', 'color', 'red')
% 
% % semilogx(fmb',DIrb, '-^')
% % title('Early band density')
% % title('Target densities')
% title('Densities')
% % h_legend = legend({'object\_playback-room', 'channel\_playback-room', 'I\_channel\_playback-room', 'target', 'room'} , 'Location', 'Southeast');
% % h_legend = legend({'target', 'room'} , 'Location', 'Northeast');
% h_legend = legend({'early', 'late'} , 'Location', 'Northeast');
% set(h_legend,'FontSize',14);
% set(gca,'fontsize', 14);
% xlabel('Band centre frequency (Hz)')



figure
% semilogx(fmb',DIprb, '-o')
xlim([fmb(1) fmb(end)])
ylim([0 inf])
semilogx(fmb',DItb, '-v', 'color', 'blue')
hold on;
semilogx(fmb',DVtb, '-^', 'color', 'red')
title('Target densities')
% h_legend = legend({'object\_playback-room', 'channel\_playback-room', 'I\_channel\_playback-room', 'target', 'room'} , 'Location', 'Southeast');
% h_legend = legend({'target', 'room'} , 'Location', 'Northeast');
h_legend = legend({'early', 'late'} , 'Location', 'Northeast');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')



figure
% semilogx(fmb',DVprb, '-o')
xlim([fmb(1) fmb(end)])
% ylim([0 DVp0rb(1)]); %inf %ylim([0 DVtb(1)*2.0])

% semilogx(fmb',DVp0rb, '-x')
% semilogx(fmb',DVp1rb, '-+')
% semilogx(fmb',DVtb, '-v')
semilogx(fmb',DIrb, '-v')
hold on;
semilogx(fmb',DVrb, '-^')
title('Room densities')
% h_legend = legend({'object\_playback-room', 'channel\_playback-room', 'I\_channel\_playback-room', 'target', 'room'} , 'Location', 'Northeast');
% h_legend = legend({'target', 'room'} , 'Location', 'Northeast');
h_legend = legend({'early', 'late'} , 'Location', 'Northeast');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')



figure
semilogx(fmb', DIpb .* DIrb, '-x')
hold on;
semilogx(fmb', DIpb .* DVrb + DVpb .* (DIrb + DVrb), '-+')
% semilogx(fmb', DVpb .* DIrb, '-+')
xlim([fmb(1) fmb(end)])
ylim([0 inf])
% semilogx(fmb',DIp0rb, '-x')
% semilogx(fmb',DIp1rb, '-+')
% semilogx(fmb',DItb, '-v')

% semilogx(fmb',DVtb, '-^')
% semilogx(fmb',DIrb, '-^')
% title('Early band density')
title('Predicted object\_playback-room densities')
% h_legend = legend({'object\_playback-room', 'channel\_playback-room', 'I\_channel\_playback-room', 'target', 'room'} , 'Location', 'Southeast');
% h_legend = legend({'target', 'room'} , 'Location', 'Northeast');
h_legend = legend({'early', 'late'} , 'Location', 'Northeast');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')


figure
semilogx(fmb',DIprb, '-x')
hold on;
semilogx(fmb',DVprb, '-+')
xlim([fmb(1) fmb(end)])
ylim([0 inf])
% semilogx(fmb',DIp0rb, '-x')
% semilogx(fmb',DIp1rb, '-+')
% semilogx(fmb',DItb, '-v')

% semilogx(fmb',DVtb, '-^')
% semilogx(fmb',DIrb, '-^')
% title('Early band density')
title('Actual object\_playback-room densities')
% h_legend = legend({'object\_playback-room', 'channel\_playback-room', 'I\_channel\_playback-room', 'target', 'room'} , 'Location', 'Southeast');
% h_legend = legend({'target', 'room'} , 'Location', 'Northeast');
h_legend = legend({'early', 'late'} , 'Location', 'Northeast');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')


% figure
% semilogx(fmb',DItb)
% xlim([fmb(1) fmb(end)])
% hold on;
% semilogx(fmb',DIprb_)
% title('Early band densities: target and channel-playback-room')



% figure
% semilogx(fmb', Rpb)
% xlim([fmb(1) fmb(end)])
% title('Ratio of early-to-late energy in playback')

% figure
% semilogx(fmb', Rr_tb)
% xlim([fmb(1) fmb(end)])
% ylim([0.0 1.0])
% title('Ratio of early-to-late density ratios, room to target')

% figure
% semilogx(fmb', Rpr_tb)
% xlim([fmb(1) fmb(end)])
% ylim([0.0 2.0])
% title('Ratio of early-to-late energy ratios, playback-room to target')
% 
% figure
% semilogx(fmb', Rtr_tb)
% xlim([fmb(1) fmb(end)])
% %ylim([0.0 2.0])
% title('Ratio of early-to-late energy ratios, target-room to target')


figure
%semilogx(fmb', Rpr_tb, '-o')
semilogx(fmb', (DIpb .* DIrb ./ (DVpb .* DIVrb + DIpb .* DVrb)) ./ (DItb ./ DVtb), '-o')
% semilogx(fmb', (DIpb .* DIrb ./ (DVpb .* (DIrb+DVrb) + DIpb .* DVrb)) ./ (DItb ./ DVtb), '-o')
xlim([fmb(1) fmb(end)])
ylim([0.0 2.0])
hold on
% semilogx(fmb', ERpr_tb, '-.')
semilogx(fmb', (DItb .* DIrb ./ (DVtb .* DIVrb + DItb .* DVrb)) ./ (DItb ./ DVtb), '-x')
% semilogx(fmb', (DItb .* DIrb ./ (DVtb .* (DIrb+DVrb) + DItb .* DVrb)) ./ (DItb ./ DVtb), '-x')
% % semilogx(fmb', Rp0r_tb, '-x')
% semilogx(fmb', (DIp0b .* DIrb ./ DVp0b ./ DVrb) ./ (DItb ./ DVtb), '-x')
% semilogx(fmb', Rp1r_tb, '-+')
% semilogx(fmb', Rtr_tb, '-^')
% semilogx(fmb', ERtr_tb, '--')
semilogx(fmb', (DIrb ./ DVrb) ./ (DItb ./ DVtb) , '-^')
semilogx(fmb', ones(1,Nb), '--', 'color', 'black')
title('Predicted ratio of early / late density ratios');
% h_legend = legend('object\_playback-room / target', 'object\_playback-room / target (expected)', 'channel\_playback-room / target', 'I\_channel\_playback-room / target', 'target-room / target ', 'target-room / target (expected)','room / target');
h_legend = legend('object\_playback-room / target', 'channel\_playback-room / target', 'room / target');
%set(h_legend,'FontSize',14);
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')



figure
semilogx(fmb', Rpr_tb, '-o')
xlim([fmb(1) fmb(end)])
ylim([0.0 2.0])
hold on
% semilogx(fmb', ERpr_tb, '-.')
semilogx(fmb', Rp0r_tb, '-x')
% semilogx(fmb', Rp1r_tb, '-+')
% semilogx(fmb', Rtr_tb, '-^')
% semilogx(fmb', ERtr_tb, '--')
semilogx(fmb', (DIrb ./ DVrb) ./ (DItb ./ DVtb) , '-^')
semilogx(fmb', ones(1,Nb), '--', 'color', 'black')
title('Actual ratio of early / late density ratios');
% h_legend = legend('object\_playback-room / target', 'object\_playback-room / target (expected)', 'channel\_playback-room / target', 'I\_channel\_playback-room / target', 'target-room / target ', 'target-room / target (expected)','room / target');
h_legend = legend('object\_playback-room / target', 'channel\_playback-room / target',  'room / target');
%set(h_legend,'FontSize',14);
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')




%,,,
figure
semilogx(fmb', (DIpb .* DIrb) ./ DItb, '-o')
xlim([fmb(1) fmb(end)])
ylim([0.0 2.0])
hold on
semilogx(fmb', (DVpb .* (DIrb + DVrb) + DIpb .* DVrb) ./  DVtb, '-x')
% semilogx(fmb', (DVpb .* DIVrb + DIpb .* DVrb) ./  DVtb, '-x')
% ( DVpb .* DIVrb  +  DIpb .* DVrb  ) ./ DVtb
semilogx(fmb', (DIrb ./ DVrb) ./ (DItb ./ DVtb) , '-s')
semilogx(fmb', ones(1,Nb), '--', 'color', 'black')
title('Predicted density ratios');
h_legend = legend('early object\_playback-room / target', 'late object\_playback-room / target',  'room / target');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')

figure
semilogx(fmb', DIprb ./ DItb, '-o')
xlim([fmb(1) fmb(end)])
ylim([0.0 2.0])
hold on
semilogx(fmb', DVprb ./  DVtb, '-x')

semilogx(fmb', (DIrb ./ DVrb) ./ (DItb ./ DVtb) , '-s')
semilogx(fmb', ones(1,Nb), '--', 'color', 'black')
title('Actual density ratios');
h_legend = legend('early object\_playback-room / target', 'late object\_playback-room / target',  'room / target');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')

%%

% Late decay time comparison RT60

% Skip some lower bands, because RT60 function fails when the response is
% too small.

% RT60offset = length(Vt)/2
for n = 1:Nb
    Vprb_RT60(n) = t60(Vprb(n,:), fs, 0, RT60offset);
end

for n = 1:Nb
    Vp0rb_RT60(n) = t60(Vp0rb(n,:), fs, 0, RT60offset);
end

% for n = 1:Nb
%     Vp1rb_RT60(n) = t60(Vp1rb(n,:), fs, 0, RT60offset);
% end

% for n = 1:Nb
%     Vtrb_RT60(n) = t60(Vtrb(n,:), fs, 0, RT60offset);
% end

% RT60offset = length(Vr)/2
for n = 1:Nb
    Vrb_RT60(n) = t60(Vrb(n,:), fs, 0, RT60offset);
    Vrb_RT60(b == 0) = 0;   % prevent small signal giving a valid RT60
end

% RT60offset = length(Vt)/2
for n = 1:Nb
    Vtb_RT60(n) = t60(Vtb(n,:), fs, 0, RT60offset);
end

cmap = colormap(lines);
figure
semilogx(fmb', Vtb_RT60'/1000, '-v', 'color', cmap(4,:))
xlim([fmb(1) fmb(end)])
ylim([0 5])
hold on
semilogx(fmb', Vrb_RT60'/1000, '-^', 'color', cmap(5,:))
% semilogx(fmb', Vprb_RT60'/1000, '-o')
% semilogx(fmb', Vp0rb_RT60'/1000, '-x')
title('RT60 (s)');
% h_legend = legend('target', 'room', 'object\_playback-room', 'channel\_playback-room');
h_legend = legend('target', 'room');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);


figure
semilogx(fmb', Vprb_RT60' ./ Vtb_RT60', '-o')
xlim([fmb(1) fmb(end)])
ylim([0 2])
hold on
semilogx(fmb', Vp0rb_RT60' ./ Vtb_RT60', '-x')
% semilogx(fmb', Vp1rb_RT60' ./ Vtb_RT60', '-+')
% semilogx(fmb',Vtrb_RT60 ./ Vtb_RT60, '-^')
% semilogx(fmb',Vrb_RT60 ./ Vtb_RT60, '-d')
semilogx(fmb', ones(1,Nb), '--', 'color', 'black')
title('Ratio of RT60s');
% h_legend = legend('object\_playback-room / target', 'channel\_playback-room / target', 'I\_channel\_playback-room / target', 'target-room / target', 'room / target');
h_legend = legend('object\_playback-room / target', 'channel\_playback-room / target');
set(h_legend,'FontSize',14);
set(gca,'fontsize', 14);
xlabel('Band centre frequency (Hz)')


% figure
% semilogx(fmb',Vprb_Vtb_RT60_ratio)
% xlim([fmb(1) fmb(end)])
% title('RT60 ratio, playback-room to target')



%%

vol = 1.0;

% display('Room')
% sound(IVr /max(abs(IVr))/2 *vol, fs)
% 
% input('');
% display('Target')
% sound(IVt /max(abs(IVt))/2 *vol, fs)
% 
% input('');
% display('Playback into Room')
% sound(IVpr /max(abs(IVpr))/2 *vol, fs)
% 
% input('');
% display('Channel playback into Room')
% sound(IVpr_ /max(abs(IVpr_))/2 *vol, fs)
% 
% input('');
% display('Target into Room')
% sound(IVtr /max(abs(IVtr))/2 *vol, fs)
% 
% input('');
% display('Playback')
% sound(IVp /max(abs(IVp))/2 *vol, fs)
% 
% s

%%


addpath('../../rooms/')
[W, fs_] = audioread('womanShort.wav');
W = W';



%%

WIr = convfft(Ir, W);
WIt = convfft(It, W);
WIpr = convfft(Ipr, W);
WIp0r = convfft(Ip0r, W);
% WIp1r = convfft(Ip1r, W);
% WItr = convfft(Itr, W);
WIp = convfft(Ip, W);

WVr = convfft(Vr, W);
WVt = convfft(Vt, W);
WVpr = convfft(Vpr, W);
WVp0r = convfft(Vp0r, W);
% WVtr = convfft(Vtr, W);
WVp = convfft(Vp, W);

WIVr = convfft(IVr, W);
WIVt = convfft(IVt, W);
WIVpr = convfft(IVpr, W);
WIVp0r = convfft(IVp0r, W);
% WIVp1r = convfft(IVp1r, W);
% WIVtr = convfft(IVtr, W);
WIVp = convfft(IVp, W);



%%

vol = 1.0;

% display('Direct')
% sound(W /max(abs(W))/2 *vol, fs)
% 
% input('');
% display('Room')
% sound(WIVr /max(abs(WIVr))/2 *vol, fs)
% 
% input('');
% display('Target')
% sound(WIt /max(abs(WIt))/2 *vol, fs)
% 
% input('');
% display('Playback into Room')
% sound(WIpr /max(abs(WIpr))/2 *vol, fs)
% 
% input('');
% display('Channel playback into Room')
% sound(WIp0r /max(abs(WIp0r))/2 *vol, fs)
% 
% 
% input('');
% display('Target into Room')
% sound(WItr /max(abs(WItr))/2 *vol, fs)
% 
% 
% 
% display('Target')
% sound(WVt /max(abs(WVt))/2 *vol, fs)
% 
% input('');
% display('Playback into Room')
% sound(WVpr /max(abs(WVpr))/2 *vol, fs)
% 
% input('');
% display('Channel playback into Room')
% sound(WVp0r /max(abs(WVp0r))/2 *vol, fs)
% 
% input('');
% display('Target into Room')
% sound(WVtr /max(abs(WVtr))/2 *vol, fs)

WW = W * W';



input('');
display('Target')
s = WIVt * sqrt( WW / (WIVt*WIVt') ) *vol;
sound(s, fs);
audiowrite('1. target.wav', s, fs)

input('');
display('Room')
s = WIVr * sqrt( WW / (WIVr*WIVr') ) *vol;
sound(s, fs);
audiowrite('2. room.wav', s, fs)

% input('');
% display('Target into Room')
% s = WIVtr * sqrt( WW / (WIVtr*WIVtr') ) *vol;
% sound(s, fs);
% audiowrite('3. target-room.wav', s, fs)

input('');
display('Channel-based playback into Room')
s = WIVp0r * sqrt( WW / (WIVp0r*WIVp0r') ) *vol;
sound(s, fs);
audiowrite('4. channel-room.wav', s, fs)

% input('');
% display('I-Channel playback into Room')
% s = WIVp1r * sqrt( WW / (WIVp1r*WIVp1r') ) *vol;
% sound(s, fs);
% audiowrite('5. I_channel_playback-room.wav', s, fs)

input('');
display('Object-based playback into Room')
s = WIVpr * sqrt( WW / (WIVpr*WIVpr') ) *vol;
sound(s, fs);
audiowrite('6. playback-room.wav', s, fs)




% 
% input('');
% display('Playback')
% sound(WIVp /max(abs(WIVp))/2 *vol, fs)

% input('');
% display('Direct')
% sound(W /max(abs(W))/2 *vol, fs)


sss



%%

vol = 0.5;

input('');
display('Room')
sound(Ir /max(abs(Ir))/2 *vol, fs)

input('');
display('Target')
sound(It /max(abs(It))/2 *vol, fs)

input('');
display('Playback into Room')
sound(Ipr /max(abs(Ipr))/2 *vol, fs)

input('');
display('Target into Room')
sound(Itr /max(abs(Itr))/2 *vol, fs)

input('');
display('Playback')
sound(Ip /max(abs(Ip))/2 *vol, fs)

%%

vol = 0.5;

input('');
display('Room')
sound(Vr /max(abs(Vr))/2 *vol, fs)

input('');
display('Target')
sound(Vt /max(abs(Vt))/2 *vol, fs)

input('');
display('Playback into Room')
sound(Vpr /max(abs(Vpr))/2 *vol, fs)

input('');
display('Target into Room')
sound(Vtr /max(abs(Vtr))/2 *vol, fs)

input('');
display('Playback')
sound(Vp /max(abs(Vp))/2 *vol, fs)


%%

s


%{

% Test energy of convolution.
% Simulate fractional bandwidth M/N

N = 10000;
M = 1000;

% Select band location.
A = N-M+1;  %4*M;
B = N; %A + M - 1;

n1 = rand(N,1)-.5;
n2 = rand(N,1)-.5;
f1 = fft(n1);
f2 = fft(n2);
f1 = f1(A:B);   % select band
f2 = f2(A:B);

e1 = f1' * f1 / N;
e2 = f2' * f2 / N;

% e1 = n1' * n1
% e2 = n2' * n2

nn = conv(n1,n2);
ff = fft(nn);
ff = ff(A:B);

ee = ff' * ff / N;

%ee = nn' * nn

display('N/M')
N/M

display('ee / (e1 * e2)')
ee / e1 / e2

%}


%% Band filter

% fs = 100000; %44100;
% bw = 1/3; % in octaves
% 
% fMin = 100; % Hz
% fMax = 20000;
% 
% octs = log2(fMax/fMin);
% bmax = floor(octs/bw)-1;
% 
% fc = fMin * 2.^( (1:bmax) * bw ); % centre frequencies
% % fc = fl *2^(-bw/2); % lower cutoffs
% % fu = fl *2^(+bw); % upper cutoffs
% 
% numBands = length(fc);
% 
% b = cell(numBands,1);
% a = cell(numBands,1);
% 
% figure
% 
% noise = rand(1, 10000)-0.5;
% noiseMod = zeros(1, 10000);
% for n = 1:length(fc)
% 
%     [b{n}, a{n}] = oct3dsgn(fc(n), fs, 2); 
%     %[b{n}, a{n}] = butter(2, [fl(n) fu(n)]/(fs/2), 'bandpass');
%     [h,f] = freqz(b{n}, a{n}, 16384, fs);
% 
%     hold on;
%     plot(f, 20*log10(abs(h)) );
%     
%     
%     noiseMod = noiseMod + filter(b{n}, a{n}, noise);
% end
% set(gca, 'XScale', 'log')
% ylim([-50 0])
% 
% figure;
% amp = abs(fft(noiseMod));
% plot(amp(1:5000))

s

%%

Fs = 44100;

BW = '1/3 octave';
N = 6;
F0 = 1000;
Nx = 10000;


oneThirdOctaveFilter = octaveFilter('FilterOrder', N, ...
    'CenterFrequency', F0, 'Bandwidth', BW, 'SampleRate', Fs);
F0 = getANSICenterFrequencies(oneThirdOctaveFilter);
F0(F0<20) = [];
F0(F0>20e3) = [];
Nfc = length(F0);
for i=1:Nfc
    oneThirdOctaveFilterBank{i} = octaveFilter('FilterOrder', N, ...
        'CenterFrequency', F0(i), 'Bandwidth', BW, 'SampleRate', Fs); %#ok
end




scope2 = dsp.SpectrumAnalyzer('SpectralAverages',50,'SampleRate',Fs,...
    'PlotAsTwoSidedSpectrum',false,'FrequencyScale','Log',...
    'RBWSource','Property','RBW',2000);
yw = zeros(Nx,Nfc);
tic
while toc < 15
    % Run for 15 seconds
    xw = randn(Nx,1);
    for i=1:Nfc
        oneThirdOctaveFilter = oneThirdOctaveFilterBank{i};
        yw(:,i) = oneThirdOctaveFilter(xw);
    end
    
    plot(abs(fft( sum(yw,2) )));
    
    scope2( [xw, yw, sum(yw,2)] );
    
    %scope2(sum(yw,2));
end



s

%%


fb = filterBank_create(fs);

Vt = 0.2 * (2*rand(1, 10000) -1);
Vtb = filterBank(Vt, fb);  % test re-composition of bands
Vt_ = sum(Vtb, 1); 
figure
plot(abs(fft(Vt)));
figure
plot(abs(fft(Vt_)));
sound(Vt, fs)
input('');
sound(Vt_, fs)

