function [inf, out] = relanoiborra2019_mfbtd(in,lmf,umf,style,fs)
%RELANOIBORRA2019_MFBTD - applies modulation filterbank	
%
%   Usage: [inf, out] = mfbtd(in,lmf,umf,style,fs)
%
%   Input parameters:
%     in    : input column vector.
%     lmf   : lowest modulation-filter center frequency,
%             if 0 the output of a 2nd-order
%             Butterworth lowpass filter is additionally computed.
%     umf   : highest modulation-filter center frequency,
%             for typical applications choose umf = 1500.
%             If lmf = umf only the output of a single filter
%             at lmf is computed.
%     style : determines fc of the lowpass filter: 1 = 2.5 Hz, 2 = 7.5 Hz.
%     fs    : sampling rate in Hz,
%             should be greater than 8000 Hz to avoid aliasing errors.
% 
%   Output parameters:
%     inf		   : center frequencies of the modulation filters
%     out          : each column of martrix out contains the output of
%                    a single modulation filter
%
%   applies a modulation filterbank as used in Relano-Iborra et al. 2019
%
%   Url: http://amtoolbox.org/amt-1.1.0/doc/modelstages/relanoiborra2019_mfbtd.php

% Copyright (C) 2009-2021 Piotr Majdak, Clara Hollomey, and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.1.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%   #Author: Stephan Ewert and Torsten Dau, Universitaet Oldenburg
%   #Author: Clara Hollomey (2021): adapted to the AMT
%   #Author: Piotr Majdak (2021): adapted to the AMT 1.0


if fs < 8000
   warning('sample rate lower than 8000 Hz')
end

Q = 2;
bw = 5;
ex=(1+1/(2*Q))/(1-1/(2*Q));
%lpcut = 150;
%ex=((3+sqrt(5))/2)^(1/den);

if lmf == 0
   sw = 1;
	switch style
      case 1
         startmf = 5;
         b2lpcut = 2.5;
   		wb2lp = 2*pi*b2lpcut/fs;
			[b2,a2] = local_solp(wb2lp,1/sqrt(2));
   	case 2
   		startmf = 10;
   		b2lpcut = 7.5;
      	wb2lp = 2*pi*b2lpcut/fs;
      	[b2,a2] = local_solp(wb2lp,1/sqrt(2));
   end
elseif lmf > 0
   startmf = lmf;
   sw = 2;
end

if lmf == umf
   mf = startmf;
   sw = 0;
   if lmf == 0
      sw =3;
   end
else
   if startmf >= 10
      tmp = fix(log(umf/startmf)/log(ex));
   	mf = 0:tmp;
   	mf = ex.^mf;
      mf = startmf*mf;
   else
      tmp = fix((min(umf,10) - startmf)/bw); %changed
      tmp = 0:tmp;
      mf = startmf + 5*tmp;
      tmp2 = (mf(end)+bw/2)/(1-1/(2*Q));
      tmp = fix(log(umf/tmp2)/log(ex));
   	tmp = 0:tmp;
      tmp = ex.^tmp;
      mf=[mf tmp2*tmp];
   end
end

% 150 Hz LP
[b1,a1] = butter(1,150/(fs/2));
outtmp = filter(b1,a1,in);
% size(outtmp)
% pause

% No LP
% outtmp=in; % no LP

switch sw
   case 0									% only one modulation filter
      w0 = 2*pi*mf/fs;
      if mf < 10
   		[b3,a3] = efilt(w0,2*pi*bw/fs);
      else
         [b3,a3] = efilt(w0,w0/Q);
      end
%       out = 2*filter(b3,a3,in);
      out = 2*filter(b3,a3,outtmp);
      inf = mf;
	case 1									% lowpass and modulation filter(s)
      out = zeros(size(in,1),size(in,2),length(mf)+1);
% 		out(:,:,1) = filter(b2,a2,in);
        out(:,:,1) = filter(b2,a2,outtmp);
		for i=1:length(mf)
         w0 = 2*pi*mf(i)/fs;
         if mf(i) < 10
   			[b3,a3] = local_efilt(w0,2*pi*bw/fs);
         else
            [b3,a3] = local_efilt(w0,w0/Q);
         end
%          out(:,:,i+1) = 2*filter(b3,a3,in);
            out(:,:,i+1) = 2*filter(b3,a3,outtmp);
      end
      inf = [0 mf];
   case 2									% only modulation filters
      out = zeros(length(in),length(mf));
      for i=1:length(mf)
         w0 = 2*pi*mf(i)/fs;
         if mf(i) < 10
   			[b3,a3] = local_efilt(w0,2*pi*bw/fs);
         else
            [b3,a3] = local_efilt(w0,w0/Q);
         end
%    		out(:,i) = 2*filter(b3,a3,in);
            out(:,i) = 2*filter(b3,a3,outtmp);
		end
      inf = mf;
   case 3									% only lowpass
%       out = filter(b2,a2,in);
        out = filter(b2,a2,outtmp);
      inf = 0;
end
   
% subfunctions

% complex frequency shifted first order lowpass
function [b,a] = local_efilt(w0,bw)

e0 = exp(-bw/2);

b = 1 - e0;
a = [1, -e0*exp(1i*w0)];

% second order Butterworth lowpass filter
function [b,a] = local_solp(w0,Q)

W0 = tan(w0/2);

b = [1; 2; 1];
a = [1 + 1/(Q*W0) + 1/W0^2; -2/W0^2 + 2; 1 - 1/(Q*W0) + 1/W0^2];

b = b/a(1);
a = a/a(1);

% % first order lowpass filter (from mfb2.m - MJ)
% function [b,a] = local_folp(w0)
% 
% W0 = tan(w0/2);
% 
% b = [W0, W0]/(1 + W0);
% a = [1,(W0 - 1)/(1 + W0)];
%