function y = multiband(x,Mode,LFC,HFC,CT1, CS1,AT1,RT1,CT2,CS2,AT2,RT2,CT3,CS3,AT3,RT3,fs)

%%                                Introduction
%       This is a multiband compressor function ,which is 
%       used to split the input signal to 3 different bands.
%       And these 3 bands are separated by crossover frequencies,which could
%       be decided by users.(LFC,HFC)
%       After that , the separated signals will be processed by the
%       compressor,and combined together at last . 
%
%       Where:
%         x  = Input signal 
%       Mode =Mode selection 'chebyshev' or 'butterworth' .
%             Because the different characteristic of these two IIR filters,
%             the users could make their own choices for the
%             filter to process the input signal. 
%             Please note that : the quotation mark must be
%             concluded!
%       LFC = low crossover frequency
%       HFC = high crossover frequency 
%       CT1 = The threshold of compressor for low frequency band in dB
%       (30<CT<0)
%       CS1 = The slope(ratio) of compressor for low frequency band 
%       
%       AT1 = The attack time  of compressor for low frequency band 
%       RT1 = The release time  of compressor for low frequency band 
% 
%       CT2 = The threshold of compressor for mid frequency band band in dB
%       (30<CT<0)
%       CS2 = The slope(ratio) of compressor for low frequency band 
%       AT2 = The attack time  of compressor for mid frequency band 
%       RT2 = The release time  of compressor for mid frequency band 
% 
%       CT3 =The threshold of compressor for high frequency band band in dB
%       (30<CT<0)
%       CS3 = The slope(ratio) of compressor for low frequency band 
%       AT3 = The attack time  of compressor for high frequency band 
%       RT3 = The release time  of compressor for high frequency band 
% 
%       fs =  Sampling frequency 

%%                                        Filter parameter setting 
 
WnL=LFC./(fs./2);  % Because the use of matlab built-in function 'butter' 
                   % or 'cheby2' need normalize the frequency first 
WnH=HFC./(fs./2);  % So this step is to normalize the frequency of the input 
                   % signal. fs/2 is Nyquist frequency,which means 
                   % the highest frequency of the input signal .So WnL and WnH are the two
                   % normalized frequency variables
order =3;          % determine the order of the filter, in this function ,3-order filter are be used
                   
%%                                      chebyshev 2 filter design
if  strcmp(Mode, 'chebyshev')        % To determine which kind of filter will be used ,and this means the users choose chebyshev filter
Rp = 30;                              % stopband ripple which means the ripple in dB in stopband of the filter 
[bl,al] =  cheby2(order,Rp,WnL,'low');   %determine the numerator and denominator ,which are important parameters of the filter.
% And in this case ,the matlab built-in function "cheby2" are used ,which means the filter's type is "'chebyselv 2 filter",
% "order" is the number of order of the filter."Rp " is the stopband ripple ,
% and "WnL" is the stop frequency as well as the low crossover frequency 
% which is determined by the users.And 'low' is the built-in type selection option ,which means the filter is a lowpass filter

[bm,am] =  cheby2(order,Rp,[WnL WnH]) ; %determine the numerator and the denominator of the bandpass filter
[bh,ah] =  cheby2(order,Rp,WnH,'high') ; %determine the numerator and the denominator of the highpass filter

%%                        processing the signal with filters(mono)
[Paul,channels] = size(x); % We use the size function to determine the
                           % size of the input signal in all its dimensions. We are interested in the
                           % second dimension (columns) which is the number of channels (the first dimension
                           % (rows) is the number of samples).
if  channels == 1; % We use an if statement to allow our functions to 
                  % work with mono or stereo signals,and if channels = 1 means the signal
                  % is a mono one .
    signal_low_cheby = filtfilt(bl,al,x);  %in this step ,the matlab built-in function "filtfilt"are used.
                                           %which performs zero-phase
                                           %digital filtering by processing
                                           %the input data,x.And 'bl' and 
                                           %'al' are obtained above .so in this case the input signal will be filtered by 
                                           %the low pass filter of "chebyshev 2"  
                                         
    signal_mid_cheby = filtfilt(bm,am,x);  %filter the signal by bandpass filter
    signal_high_cheby = filtfilt(bh,ah,x); %filter the signal by highpass filter

compressed_signal_low =compexp(signal_low_cheby,CT1,CS1,AT1,RT1);
% call compression function for low band
compressed_signal_low = compressed_signal_low/max(abs(compressed_signal_low));
%normalize compressed signal

compressed_signal_mid =compexp(signal_mid_cheby,CT2,CS2,AT2,RT2);
% call compression function for mid band
compressed_signal_mid = compressed_signal_mid/max(abs(compressed_signal_mid));
%normalize compressed signal

compressed_signal_high =compexp(signal_high_cheby,CT3,CS3,AT3,RT3);
% call compression function for high band
compressed_signal_high = compressed_signal_high/max(abs(compressed_signal_high));
%normalize compressed signal
compressed_signal = compressed_signal_low + compressed_signal_mid + compressed_signal_high ; 
%combine the compressed signals back together to obtain the compressed
%signal
y = compressed_signal./max(abs(compressed_signal));
%normalize compressed signal 

z = abs(x)-abs(y);              % create gain reduction vector for display only
 


%%                      processing the signal with filters(stereo)
elseif channels == 2; % in this case, the input signal is a stereo signal 
    
inputL = x(:,1);        %Split input channels into separate vectors,inputL means the first column of the input signal
inputR = x(:,2);        %Split input channels into separate vectors,inputL means the second column of the input signal
input = [inputL inputR];  
for  i = 1:2                   % we use a for loop to process the signal column by column.
    signal_low_cheby = filtfilt(bl,al,input(:,i));    %filter the signal by lowpass filter
    signal_mid_cheby = filtfilt(bm,am,input(:,i));    %filter the signal by bandpass filter
    signal_high_cheby = filtfilt(bh,ah,input(:,i));    %filter the signal by highpass filter
compressed_signal_low =compexp(signal_low_cheby,CT1,CS1,AT1,RT1);
% call compression function for low band
compressed_signal_low = compressed_signal_low/max(abs(compressed_signal_low));
%normalize compressed signal

compressed_signal_mid =compexp(signal_mid_cheby,CT2,CS2,AT2,RT2);
% call compression function for High band
compressed_signal_mid = compressed_signal_mid/max(abs(compressed_signal_mid));
%normalize compressed signal

compressed_signal_high =compexp(signal_high_cheby,CT3,CS3,AT3,RT3);
% call compression function for High band
compressed_signal_high = compressed_signal_high/max(abs(compressed_signal_high));
%normalize compressed sound

output(:,i)=compressed_signal_low+compressed_signal_mid+compressed_signal_high;
%combine the compressed signals back together to obtain the compressed
%signal

end
    y = [output(:,1),output(:,2)];
%obtain the final output signal .
a=max(abs(y));
y(:,1)=y(:,1)/a(1);
y(:,2)=y(:,2)/a(2);
z = abs(x)-abs(y);   
% create gain reduction vector for display only
else 
    error('This function is only a mono or stereo compressor'); 
% We will add an error option that lets the user know there is an error  
% if they input a file with more than three channels.
end



%%                                Butterworth  filter design
elseif strcmp(Mode, 'butterworth')     % To determine which kind of filter will be used ,
                                        %and this means the users choose butterworth filter
[bl,al] = butter(order,WnL,'low');      %determine the numerator and denominator ,which are important parameters of the filter.
% And in this case ,the matlab built-in function "butter" are used ,which means the filter's type is "'butterworth fliter",
% "order" is the number of order of the filter,and "WnL" is the stop frequency as well as the low crossover frequency 
% which is determined by the users.And 'low' is the built-in type selection option ,which means the filter is a lowpass filter
[bm,am] = butter(order,[WnL WnH]);      %determine thenumerator and the denominator of the bandpass filter
[bh,ah] = butter(order,WnH,'high');     %determine the numerator and the denominator of the highpass filterer

%%                          processing the signal with filters(mono)

[Paul,channels] = size(x); % We use the size function to determine the
% size of the input signal in all its dimensions. We are interested in the
% second dimension (columns) which is the number of channels (the first dimension
% (rows) is the number of samples).
if channels == 1; %% We use an if statment to allow our functions to 
                  % work with mono or stereo signals,and if channels = 1 means the signal
                  % is a mono one .
    signal_low_butter = filtfilt(bl,al,x);%filter the signal by lowpass filter
    signal_mid_butter= filtfilt(bm,am,x); %filter the signal by bandpass filter
    signal_high_butter = filtfilt(bh,ah,x); %filter the signal by highpass filter

compressed_signal_low =compexp(signal_low_butter,CT1,CS1,AT1,RT1);
% call compression function for low band
compressed_signal_low = compressed_signal_low/max(abs(compressed_signal_low));
%normalize compressed signal

compressed_signal_mid =compexp(signal_mid_butter,CT2,CS2,AT2,RT2);
% call compression function for mid band
compressed_signal_mid = compressed_signal_mid/max(abs(compressed_signal_mid));
%normalize compressed signal

compressed_signal_high =compexp(signal_high_butter,CT3,CS3,AT3,RT3);
% call compression function for High band
compressed_signal_high = compressed_signal_high/max(abs(compressed_signal_high));
%normalize compressed signal
compressed_signal = compressed_signal_low + compressed_signal_mid + compressed_signal_high ;
%combine the compressed signals back together to obtain the compressed signal
y = compressed_signal.*max(abs(compressed_signal))/max(abs(compressed_signal));
%normalize compressed signal
z = abs(x)-abs(y);              % create gain reduction vector for display only
   

%%                             processing the signal with filters(stereo)
elseif channels == 2;   % in this case, the input signal is a stereo signal 
    
    
inputL = x(:,1);        
inputR = x(:,2);
input = [inputL inputR];
for i = 1:2
    signal_low_butter = filtfilt(bl,al,input(:,i));
    signal_mid_butter = filtfilt(bm,am,input(:,i));
    signal_high_butter = filtfilt(bh,ah,input(:,i));
compressed_signal_low =compexp(signal_low_butter,CT1,CS1,AT1,RT1);
% call compression function for low band
compressed_signal_low = compressed_signal_low/max(abs(compressed_signal_low));
%normalize compressed signal

compressed_signal_mid =compexp(signal_mid_butter,CT2,CS2,AT2,RT1);
% call compression function for mid band
compressed_signal_mid = compressed_signal_mid/max(abs(compressed_signal_mid));
%normalize compressed signal

compressed_signal_high =compexp(signal_high_butter,CT3,CS3,AT3,RT3);
% call compression function for High band
compressed_signal_high = compressed_signal_high/max(abs(compressed_signal_high));
%normalize compressed signal

output(:,i)=compressed_signal_low+compressed_signal_mid+compressed_signal_high;

end
y = [output(:,1),output(:,2)]; %obtain the final output signal .
a=max(abs(y));
y(:,1)=y(:,1)/a(1);
y(:,2)=y(:,2)/a(2);
     %normalize the output signal
z = abs(x)-abs(y);             % create gain reduction vector for display only
        

else 
    error('This function is only a mono or stereo compressor'); 
    % We will add an error option that let the users know there is an error 
    % if they input a file with more than three channels.
end
else 
    error('please select mode from butterworth or chebyshev'); % We will add an error option that 
                                                               % let the users know there is an error 
                                                               %  if they
                                                               %  do not choose from 'butterworth' or 'chebyshev'
                                                              
end       
        


%% obtain figures for the comparison of input and output signals in time  domain
figure(1)        
        subplot(3,1,1), plot(x),title('Input'),xlabel('Samples'), ylabel('Amplitude'); 
        subplot(3,1,2), plot(y), title('Output'),xlabel('Samples'), ylabel('Amplitude'); 
        subplot(3,1,3), plot(z), title('Gain Reduction'),xlabel('Samples'), ylabel('Reduction'); 
      
%% obtain the figures  of input and output signals in frequency domain
figure(2)


freq_x = fft(x(:,1));
fft_frequencies_vector = linspace(0, fs, length(x(:,1))); % vector
% of frequencies from 0 to the sampling frequency in steps quantity equal
% to the length of our signal
magnitude_spectrum = abs(freq_x); % obtain the magnitude
% response from complex spectrum
dbmag_spectrum = 20 .* log10(magnitude_spectrum); % obtain the dB
% magnitude from a magnitude response
Nyquist_index = length(freq_x) ./ 2; % obtain the index
% for the Nyquist frequency
subplot(2,1,1),
semilogx(fft_frequencies_vector(1:Nyquist_index),dbmag_spectrum(1:Nyquist_index)); % plot magnitude response using
% logarithmic scale up to the Nyquist frequency

freq_y = fft(y(:,1));
fft_frequencies_vector = linspace(0, fs, length(y(:,1))); % vector
% of frequencies from 0 to the sampling frequency in steps quantity equal
% to the length of our signal
magnitude_spectrum = abs(freq_y); % obtain the magnitude
% response from complex spectrum
dbmag_spectrum = 20 .* log10(magnitude_spectrum); % obtain the dB
% magnitude from a magnitude response
Nyquist_index = length(freq_y) ./ 2; % obtain the index
% for the Nyquist frequency
subplot(2,1,2),
semilogx(fft_frequencies_vector(1:Nyquist_index),dbmag_spectrum(1:Nyquist_index)); % plot magnitude response using
% logarithmic scale up to the Nyquist frequency


%% obtain the frequency response of the filter
figure(3)
 fvtool(bl,al); %use matlab built-in function "fvtool",to obtain the frequency response of the low-pass filter
 fvtool(bm,am); %use matlab built-in function "fvtool",to obtain the frequency response of the band-pass filter
fvtool(bh,ah); %use matlab built-in function "fvtool",to obtain the frequency response of the high-pass filter


end

