
function [out]=tapedelay(x,drygain,delaygain,modgain,delay,feedback,cutofflowpass,cutoffhighpass,filterorderlowpass,filterorderhighpass,saturation,Modfreq,width,fs)
%Tape Delay function, with adjustable parameters pertaining to: lowpass
%filtering, highpass filtering, saturation and delay modulation of delayed signal. The sonic
%goal of this function is to emulate an analogue tape delay.
%
%   example:[out]=tapedelay(x,0.8,0.4,0.4,0.3,0.4,20000,1000,100,100,2.2,6,0.01,fs);
%           
%           where: x = signal in
%            drygain = the level of the dry(uneffected signal) 
%            wetgain = the level of the delayed signal                   
%              delay = the amount of delay time (in seconds)
%           feedback = feedback gain, amount of delayed signal to be fedback.  
%     cuttofflowpass = the cutt off frequency for the lowpass filter of the delayed signal
%        filterorder = the order of the lowpass filter: 1 < dB/octave, 10 > dB/octave (a steeper slope)
%         saturation = the amount of saturation(distortion) added to
%         delayed signal (between 1-3) 3 being quite distorted.
%            Modfreq = the frequency that the delay will be modulated (low
%            frequency recommended (2-10 Hz)
%              width = modulation width
%                 fs = sampling frequency (44100)



dbstop if error



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MODULATING DELAY from DAFX (vibrato)



%width cannot be greater than delay

ya_alt=0;
bdelay=round(delay*fs); % basic delay in # samples
bwidth=round(width*fs); % basic width in # samples
if bwidth>bdelay %setting limits to parameters 
    error('delay greater than basic delay')
    return;
end
MODFREQ=Modfreq/fs; % modulation frequency in # samples
LEN=length(x) ; % # of samples in WAV-file
L=2+bdelay+bwidth*2; % length of the entire delay
Delayline=zeros(L,1); % memory allocation for delay
for d=1:(LEN-1); %setting vector dimensions of 'd' 
    MOD=sin(MODFREQ*2*pi*d) ; %sine wave used for modulation
    ZEIGER=1+bdelay+bwidth*MOD; % incomperating delay and width with the modulating sine wave
    i = floor (ZEIGER); %using floor function to round towards negative infinity
    frac=ZEIGER-i; %rounding down
    Delayline=[x(d);Delayline(1:L-1)]; 
    %allpassinterpoaltion
    y(d,1)=(Delayline(i+1)+(1-frac)*Delayline(i)-(1-frac)*ya_alt); %vibrato output(modulating delay)
    ya_alt=y(d,1);
end
        
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Delay
sampleDelay =ceil(delay*fs); % delay in sec * samp/sec = samples
delayedSignal = [zeros(sampleDelay,1);x(1:end);zeros((6*fs)-sampleDelay,1)];
% first delay by padding beginning with proper # zeros
% then zero pad the end of the delayed signal for the rest of the delays
% that come from feedback

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_extra_zeros = zeros(length(delayedSignal)-length(y),1);%matching matrix lengths of y and delayed signal
ypadded = [y;y_extra_zeros];%adding zeros to match length of delayedSignal

ypadded = ypadded./max(abs(ypadded)); %normalizing to match signals
delayedSignal = delayedSignal./max(abs(delayedSignal));%normalizing to match signals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Summed Signal

summedsignal = (ypadded*modgain) + (delayedSignal*delaygain); %summed signal of both the 
%modulated delay and the delayed signal to be processed together through the next stage(filter/loop delay)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Filter summed Signal put through feedback loop
%
%Lowpass/highpass filter
%h=fdesign.lowpass/(or 'highpass')('N,Fc',filterorder,frequency,samplingfrequency)
%hd = design(h);
%inputsig = fftfilt(hd.Numerator,inputsig);
%
%Saturate Filtered 
%saturation= gainfactor(0>3) * atan(k * datain/gainfactor);
%k=constant used to increase or decrease distortion;
%delayedSignal=saturation * atan (k * delayedSignal/saturation);
%
H= fdesign.highpass('N,Fc',filterorderhighpass,cutoffhighpass,fs);%calling matlab function highpass filter
HD= design(H); %design function creates the high pass filter to be convolved with signal.
h = fdesign.lowpass('N,Fc',filterorderlowpass,cutofflowpass,fs);%calling matlab function lowpass filter
hd = design(h); %design function creates the low pass filter to be convolved with signal.
n=1; %iterator variable for # delays
while feedback >= 0.0001 %going to continually decrease feedback inside loop
    summedsignal=summedsignal+feedback*[zeros(sampleDelay*n,1); 
    summedsignal(1:end-sampleDelay*n)]; 
    summedsignal = fftfilt(hd.Numerator,summedsignal);%filtered delay signal low pass
    summedsignal = fftfilt(HD.Numerator,summedsignal);%filtered delay signal high pass
    summedsignal=saturation * atan (1 * summedsignal/saturation);%saturated (distortion added) to delay signal
    %summedsignal=(summedsignal./(max(abs(summedsignal)))); %normalized
    %summedsignal= summedsignal * (drygain*0.99); %leveling of normalized signal(multiplying
    %by drygian so that the ratio between wet and dry is the same, also max(abs(x))=0.99 hence using that value)
    feedback=feedback^2;%feedback exponentially decaying
    n=n+1;
end


x_extra_zeros = zeros(length(summedsignal)-length(x),1);%matching matrix lengths of x and summedsignal
xpadded = [x;x_extra_zeros]; %vector of zeros added to end of input signal to match summedsignal matrix dimensions
xpadded = xpadded./max(abs(xpadded)); %dry singal normalized 
xpadded = xpadded * drygain; %normalized signal multiplied by a gain factor of preference.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%output signal direct combined with delayed.
%gain factors for wet and dry


out = xpadded + summedsignal;%final output of feedforward and summedsignal signal both with gain factors

out = (out./max(abs(out))) * 0.8; % normalizing final signal out with a reducing gain factor to give it a reasonable output level.

% removing excess zeros from the end
content = find(out,1,'last'); %find function used to locate the last non zero value
out=out(1:content); % the output vector now has no excess zeros at the end of the signal having used the above line.

end