function out = comp (X, fs, AT, RT, Comp, Cth, Lth, subp)
% Syntax : comp (X, Fs , AT, RT, Comp, Cth, Lth)
% X: Input sound
% fs : Sample frequency
% AT: Attack time (0.00016 - 2.6 sec)
% RT: Release time (0.001 - 5 sec)
% Comp: Compression factor (0 - 1)
% Cth: Compression threshold (in dB)
% Lth: Limiter threshold (in dB)
% subp: Plot graphs, 1 for yes, 0 for no.
%
% Based on the limiter.m example p99 in DAFX.

if nargin < 8
    subp = 0;
elseif nargin < 7
    error ('Incorrect number of input arguments. See help comp for more info.');
end;

%Check so that all parameter values are valid.
if 0.00016 > AT || AT > 2.6
    error ('Bad attack time');
end;
if 0.001 > RT || RT > 5
    error ('Bad release time');
end;
if 0 > Comp || Comp > 1
    error ('Not a good compression factor');
end;
if 20 < Cth || Cth < -120
    error ('Not a good compression threashold value');
end;
if 20 < Lth || Lth < -120 || Lth < Cth
    error ('Not a good limiter threashold value');
end;


%Convert from dB.
Cth = 10^(Cth/20);
Lth = 10^(Lth/20);

%Calculate at what level the limiter kicks in.
Limit = (Lth - Cth)*Comp + Cth;

%Attack and release times in number of samples.
AT = 1-exp(-2.2* fs/AT);
RT = 1-exp(-2.2* fs/RT);

if size(X, 2) > size (X,1)
    X = X';
end;

%Initialize peak detector vector and gain factor vector.
xd = zeros (length(X),1);
f = zeros (length(X),1);

%Go through the sound vector and compare the signal values with the
%threshold values. Set gain factor vector and multiply.
for n = 2:length(X)
    
    a = abs (X(n,1)) - xd (n - 1);
    if a < 0
        a = 0;
    end;
    xd(n) = xd(n - 1)*(1-RT)+AT*a;
    
    if Cth < xd(n) && xd(n) < Lth
            f(n) = 10^(-Comp*(log10(xd(n))-log10(Cth)));
        elseif xd(n) >= Lth
            f(n) = 10^(-1*(log10(xd(n))-log10(Limit)));
    else
            f(n) = 1;
    end;
end;

%Multiply with the gain factor vector.
out (:,1) = X(:,1).*f;
%Maximize
out(:,1) = out(:,1)/max(abs(out(:,1)));

%Plot results
if subp == 1
    subplot(2,2,1), plot(X(:,1)), title ('Input channel 1'), xlabel('Samples'), ylabel('Amplitude');
    subplot(2,2,2), plot(f), title('Gain '), xlabel('Samples'), ylabel('Gain factor');
    subplot(2,2,3), plot(xd), title('Peak detector'), xlabel('Samples'), ylabel('Peak amplitude');
    subplot(2,2,4), plot(out(:,1)), title('Out'), xlabel('Samples'), ylabel('Amplitude');
end;

%If the sound is in stereo, repeat the same as above for the second
%channel.
if size(X,2) > 1
    xd = zeros(length(X),1);
    f = zeros(length(X),1);
    for n = 2:length(X)
        a = abs(X(n,2)) - xd(n-1);
        if a<0
            a = 0;
        end;
        xd(n) = xd(n-1)*(1-RT)+AT*a;
    
        if Cth < xd(n) && xd(n) < Lth
            f(n) = 10^(-Comp*(log10(xd(n))-log10(Cth)));
        elseif xd(n) >= Lth
            f(n) = 10^(-1*(log10(xd(n))-log10(Limit)));
        else
            f(n) = 1;
        end;
    end;
    out(:,2) = X(:,2 ).*f;
    out (:,2) = out(:,2)/max(abs(out(:,2)));

    if subp == 1
        figure;
        subplot(2,2,1), plot (X(:,2)), title('Input channel 2'), xlabel('Samples'), ylabel('Amplitude');
        subplot(2,2,2), plot(f), title('Gain'), xlabel('Samples'), ylabel('Gain factor');
        subplot(2,2,3), plot(xd), title('Peak detector'), xlabel('Samples'), ylabel('Peak amplitude');
        subplot(2,2,4), plot(out(:,2)), title('Out'), xlabel('Samples'), ylabel('Amplitude');
    end;
end;