function y = multi_band_limiter(x, lt1, lt2, lt3, fs);

%% A stereo 3-band multiband limiter
%
%       This function takes the input 'x' and splits it into three separate
%       band limited signals. The crossover points being 150 Hz and 2000 Hz.
%       These signals are independantly processed through a compression 
%       algorithm and then combined back together.
%       A plot is then generated to show the response of the filters and
%       the transfer function of the system. The sound is then played.
%
%       The function is able to detect whether the input is stereo or mono.
%       If the input is stereo it processes the Left and Right channels
%       separately through then writes then back as a stereo file. If the
%       file is stereo the plot only shows the transfer function of the
%       left channel.
%
%       x =         Input signal (Column Vector, 1 column for mono, 2
%                   columns for stereo)
%
%       lt1 =       treshold value for band pass 1
%       lt2 =       treshold value for band pass 2
%       lt3 =       treshold value for band pass 3
%       fs =        Sampling Frequency (Hz)

        %% Design Filters

xover1 = 200;    %Set crossover point between LF and MF
xover2 = 2000;   %Set crossover point between MF and HF
        
D1Fp = xover1-(xover1./10);     %Set passband end frequency for lowpass filter (LF)
D1Fst = xover1+(xover1./0.5);   %Set stopband end frequency for lowpass filter (LF)

D2Nb = 50;                      %Set IIR numerator order for bandpassfilter (MF) 
D2Na = 10;                      %Set IIR denominator order for bandpassfilter (MF)
D2Fst1 = xover1-(xover1./5);    %Set low stopband end frequency for bandpass filter (MF)
D2Fp1 = xover1+(xover1./5);     %Set low passband end frequency for bandpass filter (MF)
D2Fst2 = xover2+(xover2./1.5);  %Set high stopband end frequency for bandpass filter (MF)
D2Fp2 = xover2-(xover2./5);     %Set high passband end frequency for bandpass filter (MF)

D3Fp = xover2+(xover2./3);      %Set passband end frequency for highpass filter (HF)
D3Fst = xover2-(xover2./1.5);   %Set stopband end frequency for highpass filter (HF)
     
D1 = fdesign.lowpass('Fp,Fst,Ap,Ast', D1Fp, D1Fst, 1, 60, fs);                                  %Create lowpass filter (LF)
D2 = fdesign.bandpass('Nb,Na,Fst1,Fp1,Fp2,Fst2',D2Nb, D2Na, D2Fst1,D2Fp1,D2Fp2,D2Fst2, fs);     %Create bandpass filter (MF)
D3 = fdesign.highpass('Fst,Fp,Ast,Ap', D3Fst, D3Fp, 60, 1, fs);                                 %Create highpassfilter (HF)

Hd1 = design(D1,'IIR'); %Implement IIR version of LF filter
Hd2 = design(D2,'IIR'); %Implement IIR version of MF filter
Hd3 = design(D3,'FIR'); %Implement FIR version of HF filter


    %% Set limiter time settings
    
at = 9;    %Set limiter attack time
rt = 500;   %Set limiter release time

        %% Test for number of input channels
        
stereo_input_test = size(x); %Returns dimensions of input array to determine if stereo or mono

if stereo_input_test(:,2)==2 %If number of columns in stereo_input_test = 2 then input is stereo
    
    inputL = x(:,1);        %Split input channels into separate vectors
    inputR = x(:,2);
    
     %% Process each channel separately through Multiband Limiter if signal is Stereo
for input = [inputL inputR];
    
    band1 = filter(Hd1, input); %Run input throught LF filter
    band2 = filter(Hd2, input); %Run input throught MF filter
    band3 = filter(Hd3, input); %Run input throught HF filter

    band1limit = newlimit(band1, lt1, at, rt, fs);    %Run LF signal through compressor 
    band2limit = newlimit(band2, lt2, at, rt, fs);    %Run MF signal through compressor
    band3limit = newlimit(band3, lt3, at, rt, fs);    %Run HF signal through compressor
    
    if input==inputL
    output(:,1) = band1limit+band2limit+band3limit;    %Recombine compressed Left signals
    end
    if input==inputR
    output(:,2) = band1limit+band2limit+band3limit;    %Recombine compressed Right signals
    end
end

output(:,1) = 0.99.*output(:,1)./max(abs(output(:,1)));  %Normalize output multiplied by 0.99 to avoid clipping
output(:,2) = 0.99.*output(:,2)./max(abs(output(:,2)));

wavwrite(output, fs, 'outputstereolimit.wav'); %Write wave file of output

end

    %% Process Mono channel through Multiband Limiter

if stereo_input_test(:,2)==1
    input = x;
    
    band1 = filter(Hd1, input); %Run input throught LF filter
    band2 = filter(Hd2, input); %Run input throught MF filter
    band3 = filter(Hd3, input); %Run input throught HF filter

    band1limit = newlimit(band1, lt1, at, rt, fs);    %Run LF signal through compressor newlimit(x,lt,at,rt,fs)
    band2limit = newlimit(band2, lt2, at, rt, fs);    %Run MF signal through compressor
    band3limit = newlimit(band3, lt3, at, rt, fs);    %Run HF signal through compressor
    
    output = band1limit+band2limit+band3limit; %Recombine compressed signals
    
    wavwrite(output, fs, 'outputmonolimit.wav'); %Write wave file of output
end

        %% Generate Plot
    
H1=freqz(Hd1);          %Returns magnitude data from LF filter
H2=freqz(Hd2);          %Returns magnitude data from MF filter
H3=freqz(Hd3);          %Returns magnitude data from HF filter
H1mag=abs(H1);          %Convert magnitude data from LF filter from complex form to absolute values 
H2mag=abs(H2);          %Convert magnitude data from MF filter from complex form to absolute values 
H3mag=abs(H3);          %Convert magnitude data from HF filter from complex form to absolute values 
H1dB=20.*log10(H1mag);  %Convert absolute magnitude data to dB for LF filter
H2dB=20.*log10(H2mag);  %Convert absolute magnitude data to dB for MF filter
H3dB=20.*log10(H3mag);  %Convert absolute magnitude data to dB for HF filter

in= freqz(x(:,1));           %Returns magnitude data for input
out= freqz(output(:,1));     %Returns magnitude data for output
inmag= abs(in);              %Convert magnitude data from input from complex form to absolute values
outmag= abs(out);            %Convert magnitude data from output from complex form to absolute values

resp= outmag./inmag;        %Get the magnitude transfer function of the system

respdB= 20.*log10(resp);    %Convert transfer function to dB form
len_resp = length(respdB);  %Get length of transfer function
len_filt = length(H1dB);    %Get length of LF filter

frequencies_filt = fs .* ((1:len_filt)-1) ./ len_filt;  %Generate X-Axis values for filter data

frequencies_resp = fs .* ((1:len_resp)-1) ./ len_resp;  %Generate X-Axis values for transfer function

figure (1)
outfig = figure (1);
plot(frequencies_filt, H1dB)        % Plot all filters and transfer function on the same plot
hold on
plot (frequencies_filt, H2dB, 'r')
plot (frequencies_filt, H3dB,'g')
plot(frequencies_resp, respdB, 'k')

title('Filter Response & Transfer Functon');    %Set graphical parameters for the plot
ylabel('Magnitude (dB)');
xlabel('Frequency (Hz)');
set(gca, 'XScale', 'log');
set(gca, 'YLim', [-80 10]);
set(gca, 'XLim', [20 20000]);
legendPlot = legend('LF','MF','HF','Transfer Function');
set(legendPlot, 'Location', 'NorthEastOutside');
set(outfig, 'Position', [0 0 1200 800]);
    
        %% Play Sound
sound(output,fs)        %Play output sound

end
    