function [outsig, mfc, step] = king2019_modfilterbank(insig,fs,varargin)
%KING2019_MODFILTERBANK modulation filterbank used by King et al. 2019
%
%   Usage: [outsig, mfc, step] = king2019_modfilterbank(insig,fs,varargin)
% 
%   Input parameters:
%     insig: Input signal(s)
%     fs     : Sampling rate in Hz
%     fc     : Centre frequencies of the input signals
%     mflow  : minimum modulation centre frequency in Hz
%     mfhigh : maximum modulation centre frequency in Hz
%     N      : Number of logarithmically-spaced (between fmin and fmax) 
%              modulation filters (default N = 10)
%     Qfactor: Quality factor of the filters (default Qfactor = 1).
%
%   Output parameters:
%     outsig: Modulation filtered signals 
%     mfc   : Centre frequencies of the modulation filters
%     step  : Contains some intermediate outputs.
%
%   KING2019_MODFILTERBANK calculates the modulation filterbank used by King et al. 2019
%
%   Url: http://amtoolbox.org/amt-1.3.0/doc/modelstages/king2019_modfilterbank.php

%   #Author: Leo Varnet and Andrew King (2020)
%   #Author: Alejandro Osses (2020) Original implementation for the AMT
%   #Author: Clara Hollomey (2021) Adapted for AMT
%   #Author: Piotr Majdak (2021) Further adaptations to AMT 1.0

% This file is licensed unter the GNU General Public License (GPL) either 
% version 3 of the license, or any later version as published by the Free Software 
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and 
% at <https://www.gnu.org/licenses/gpl-3.0.html>. 
% You can redistribute this file and/or modify it under the terms of the GPLv3. 
% This file is distributed without any warranty; without even the implied warranty 
% of merchantability or fitness for a particular purpose. 

definput.keyvals.mfc=[];
[flags,kv]=ltfatarghelper({},definput,varargin);

% first order modulation Butterworth lowpass filter with a cut-off
% frequency of 150 Hz. This is to remove all modulation frequencies
% above 150 Hz. The motivation behind this filter can be found in kohlrausch2000
if flags.do_LP_150_Hz % modbank_LPfilter
    [b_highest,a_highest] = butter(1,150/(fs/2));
    insig = filter(b_highest,a_highest,insig);
end

% Parameters modulation filter:
mflow           = kv.mflow; 
mfhigh          = kv.mfhigh; 
modbank_Nmod    = kv.modbank_Nmod;
modbank_Qfactor = kv.modbank_Qfactor;
    
% -------------------------------------------------------------------------
% -- 1. modulation_filterbank.m
if isempty(modbank_Nmod)
    step_mfc = (sqrt(4*modbank_Qfactor^2+1)+1)/(sqrt(4*modbank_Qfactor^2+1)-1);
    logfmc = log(mflow):log(step_mfc):log(mfhigh);
    
    modbank_Nmod = length(logfmc);
else
    logfmc = linspace(log(mflow), log(mfhigh), modbank_Nmod); %log(fmin):log((sqrt(4*Qfactor^2+1)+1)/(sqrt(4*Qfactor^2+1)-1)):log(fmax);
end
mfc = exp(logfmc);

for ichan = 1:modbank_Nmod
    flim(ichan,:) = mfc(ichan)*sqrt(4+1/modbank_Qfactor^2)/2 +  [-1 +1]*mfc(ichan)/modbank_Qfactor/2; %sqrt((fmc(ichan)/Qfactor)^2+8*fmc(ichan))/2 + [-1 +1]*(fmc(ichan)/(2*Qfactor));% [fmc(ichan)*((sqrt(5)-(1/Qfactor))/2) fc(ichan)*((sqrt(5)+(1/Qfactor))/2)];
    [b(ichan,:),a(ichan,:)] = butter(2,2*[flim(ichan,:)]/fs); 
end

% -------------------------------------------------------------------------
% 2. apply_filterbank.m

% E_mod = apply_filterbank(BB, AA, inoutsig, cfg.modbank_filtfilt);

Nchannels=size(b,1);
if size(a,1)~= Nchannels
    error('Number of lines in ''a'' and ''b'' must be equal');
end

Nsamples=size(insig,1);
Nsignals=size(insig,2);

outsig=zeros(Nsamples,Nsignals,Nchannels);

for ii=1:Nchannels
    outsig(:,:,ii) = filter(b(ii,:),a(ii,:),insig);
end

if nargout >= 3 % Intermediate outputs
    step.a = a;
    step.b = b;
    step.fs_design = fs;
    
    step.E_mod = outsig;
    step.fmc = mfc;
end
    
% -------------------------------------------------------------------------
% 3. Phase insensitivity
% ---
if flags.do_phase_insens_hilbert
  amt_disp('  Phase insensitivity (hilbert)',flags.disp);
    
    phase_insens_cut = kv.phase_insens_cut; % Hz
    
    N_fc = size(outsig,2);
    for i=1:modbank_Nmod
        if mfc(i)>phase_insens_cut
            for j=1:N_fc
                % hilbert
                outsig_env = local_hilbert_extraction(squeeze(outsig(:,j,i)));
                scal_factor = rms(outsig(:,j,i))/ rms(outsig_env); 
                outsig(:,j,i) = outsig_env * scal_factor;
            end
       end
    end
    
    if nargout>=3
        step.E_phase_ins = outsig;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function outsig_env = local_hilbert_extraction(D)
% Author: Leo Varnet 2016 - last modified 11/10/2018

if isreal(D)
    hilbert_responses = hilbert(D);
else
    hilbert_responses = D;
end

outsig_env = abs(hilbert_responses);