function [y] = distortion(x,gain,tone,fs)
%distortion simulates a guitar distortion
%x = input signal
%fs = sampling rate
%drive = level of distortion, best results I found around 5-15dB
%tone = linearly mixed bandpass filter, 0-1 where 0 is more bass and 1 is more tontreble 

% drive - normalize input signal and apply gain
s=x;
s=s/max(abs(s));
s=drive*s;


% Nonlinearity - ends asymptotically at 1 and is mostly linear close to values at +/- 0
s=erf(s);


% Tone - 2 overlapping 1st order butterworth filters
% From several test trials, I chose to lowpass at 1.3kHz and highpass at 260Hz
[ZL,PL,KL]=butter(1,1300/(fs/2),'low');
[ZH,PH,KH]=butter(1,260/(fs/2),'high');
[BL,AL]=zp2tf(ZL,PL,KL);
[BH,AH]=zp2tf(ZH,PH,KH);


% Output - filter based on position/value of tone control
y=(1-tone)*filter(BL,AL,s)+tone*filter(BH,AH,s);


end