function [referenceMatchingRms,refm_cpl_fft,refs_cpl_fft]=MEMsumref(reference,Feqz)%#codgen
% Initial Paremeters
% Feqz=44100;
fftSize=4096;
% convert windows size
timeArea=15;
areaSize = timeArea * Feqz;
refsiz = size(reference, 1);
%% Matching RMS
% setup ref segment
refm= (reference(:, 1) + reference(:, 2)) ./ 2;
refs= (reference(:, 1) - reference(:, 2)) ./ 2;
% sizes of parts
refN = fix(refsiz / areaSize) + 1;
refPn = fix(refsiz / refN);
% get rms of ref
refm_RMS = zeros(1, refN);
for j = 1 : 1 : refN
    refm_RMS(j) = rms(refm( (refPn * (j - 1) + 1) : (refPn*j) ));
end
% May add array preallocation here in future
% position beyond the average RMS mid
ref_RMSB_P= find(refm_RMS>=rms(refm_RMS));
ref_extra=refm_RMS>=rms(refm_RMS);
refm_an = zeros( sum(ref_extra, 2), refPn);
refs_an = zeros( sum(ref_extra, 2), refPn);
for i = 1 : 1 : sum(ref_extra, 2)
    refm_an(i, :) = refm( (refPn * (ref_RMSB_P(i) - 1) + 1) : (refPn * ref_RMSB_P(i)) );
    refs_an(i, :) = refs( (refPn * (ref_RMSB_P(i) - 1) + 1) : (refPn * ref_RMSB_P(i)) );
end
ref_an_rms = zeros(1,sum(ref_extra, 2));
for i = 1 : 1 :sum(ref_extra, 2)
    ref_an_rms(i) = refm_RMS(ref_RMSB_P(i));
end
referenceMatchingRms = rms(ref_an_rms);
%% Matching FR
refm_fft_N = fix(refPn / fftSize);
refs_fft_N = refm_fft_N;
%----refm MID channel----
refm_av_fft = zeros(sum(ref_extra, 2), fftSize);
for j = 1 : 1 : sum(ref_extra, 2)
    refm_fft = zeros(refm_fft_N, fftSize);
    for i = 1 : 1 : refm_fft_N
        refm_fft(i, :) = abs(fft(refm_an(j, (fftSize * (i - 1) + 1) : (fftSize * i) ))/fftSize);
    end
    if size(refm_fft, 1) == 1
        refm_av_fft(j, :) = refm_fft(1,:);
    else
        refm_av_fft(j, :) = sum(refm_fft) ./ refm_fft_N;
    end   
end
% ref mid match reference mid
refm_av_fft(refm_av_fft <  0.0001) = 0.0001;

% average FFT of reference mid channel
if size(refm_av_fft, 1) == 1
    refm_cpl_fft = refm_av_fft;
else
    refm_cpl_fft = sum(refm_av_fft) ./ size(refm_an, 1);
end
%----refs SID channel----
refs_av_fft = zeros(size(refs_an, 1), fftSize);
for j = 1 : 1 : size(refs_an, 1)
    refs_fft = zeros(refs_fft_N, fftSize);
    for i = 1 : 1 : refs_fft_N
        refs_fft(i, :) = abs(fft(refs_an(j, (fftSize * (i - 1) + 1) : (fftSize * i) ))/fftSize);
    end
    if size(refs_fft, 1) == 1
        refs_av_fft(j, :) = refs_fft(1,:);
    else
        refs_av_fft(j, :) = sum(refs_fft) ./ refs_fft_N;
    end 
end
% average FFT of reference SID channel
if size(refs_av_fft, 1) == 1
    refs_cpl_fft = refs_av_fft;
else
    refs_cpl_fft = sum(refs_av_fft) ./ size(refs_an, 1);
end
% disp('Calculating side matching FFT...')
refs_cpl_fft(refs_cpl_fft < 0.0001) =0.0001;
end