%%%%%%%%%%%%%%%%%%% Multiband Delay By Trenton Smith 2016%%%%%%%%%%%%%%%
 % This script is designed to take one mono or stereo file and with
 % controllable parameters be able to adjust the delay time for three
 % separate frequency bands. The delay design and filter options are
 % below. This script does take a longer time depending on the delay and
 % audio input length
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Load Audio File
 [x, fs] = audioread('clap.wav');

 %%%%%%%%%%%%%%%%%%%%%%%%%%%Delay Design by BPM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % The function 'tempo' takes BPM and coverts it into the corrisponding 
% miliseconds delay value for 1/4,1/8,1/16 or 1/4(triplet) note delay

% Choose the total BPM to sync 

% Chose a mode of delay
% Mode 1 = 1/4 delay
% Mode 2 = 1/8 delay
% Mode 3 = 1/16 delay
% Mode 4 = 1/4(triplet) delay

% If a delay time greater than 4 is selected that number now represents the
% delay in Miliseconds.

% If 0 is selected delay time will = 5ms (smallest delay acceptable) 
 

Delay_time1 = tempo(300,120); % Set mode and BPM for LP Filter (mode,bpm)
Delay_time2 = tempo(500,120); % Set mode and BPM for Bandpass Filter(mode,bpm)
Delay_time3 = tempo(1000,120); % Set mode and BPM for HP Filter(mode,bpm)

Fb1 = 0.5  ;    % Set feedback for LP Delay
Fb2  = 0.5 ;  % Set feedback for Bandpass Delay
Fb3 = 0.5;  % Set feedback for HP Delay

%%%%%%%%%%%%%%%%%%%%%%%%%%Filter Design%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xover1 = 200;    %Set crossover point between LP and BP
xover2 = 2000;   %Set crossover point between BP and HP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%Ready to run the script!!%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% WARNING!! LONGER THE DELAY - LONGER THE PROCESSING TIME %%%%%%%%%%


D1Fp = xover1-(xover1./10);     %Set passband end frequency for lowpass filter (LP)
D1Fst = xover1+(xover1./0.5);   %Set stopband end frequency for lowpass filter (LP)

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

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

Hd1 = design(D1,'IIR'); %Implement IIR version of LP filter
Hd2 = design(D2,'IIR'); %Implement IIR version of BP filter
Hd3 = design(D3,'FIR'); %Implement FIR version of HP filter

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);
    input=x;

     %% Process each channel separately through Multiband delay if signal is Stereo

    
band1L = filter(Hd1, inputL); %Run input throught LP filter
band2L = filter(Hd2, inputL); %Run input throught BP filter
band3L = filter(Hd3, inputL); %Run input throught HP filter

band1R = filter(Hd1, inputR); %Run input throught LP filter
band2R = filter(Hd2, inputR); %Run input throught BP filter
band3R = filter(Hd3, inputR); %Run input throught HP filter

band1delayL = my_delay(band1L, Fb1, Delay_time1, fs);   %Run LP signal through delay
band2delayL = my_delay(band2L, Fb2, Delay_time2, fs);     %Run BP signal through delay
band3delayL = my_delay(band3L, Fb3, Delay_time3, fs);     %Run HP signal through delay

band1delayR = my_delay(band1R, Fb1, Delay_time1, fs);   %Run LP signal through delay
band2delayR = my_delay(band2R, Fb2, Delay_time2, fs);     %Run MP signal through delay
band3delayR = my_delay(band3R, Fb3, Delay_time3, fs);     %Run HP signal through delay

% obtain the length of each delayed band 
banddelayL_L(1) = length(band1delayL);
banddelayL_L(2) = length(band2delayL);
banddelayL_L(3) = length(band3delayL);

% obtain the length of each delayed band 
banddelayR_L(1) = length(band1delayR);
banddelayR_L(2) = length(band2delayR);
banddelayR_L(3) = length(band3delayR);

% calculate the length difference between each delayed band
banddelay_maxL = max(banddelayL_L);
banddelayL_D(1) = banddelay_maxL-banddelayL_L(1);
banddelayL_D(2) = banddelay_maxL-banddelayL_L(2);
banddelayL_D(3) = banddelay_maxL-banddelayL_L(3);

% calculate the length difference between each delayed band
banddelay_maxR = max(banddelayR_L);
banddelayR_D(1) = banddelay_maxR-banddelayR_L(1);
banddelayR_D(2) = banddelay_maxR-banddelayR_L(2);
banddelayR_D(3) = banddelay_maxR-banddelayR_L(3);

% padd with zeros to match the length of all delayed bands
band1delayL = [band1delayL';zeros(banddelayL_D(1),1)];
band2delayL = [band2delayL';zeros(banddelayL_D(2),1)];
band3delayL = [band3delayL';zeros(banddelayL_D(3),1)];

% padd with zeros to match the length of all delayed bands
band1delayR = [band1delayR';zeros(banddelayR_D(1),1)];
band2delayR = [band2delayR';zeros(banddelayR_D(2),1)];
band3delayR = [band3delayR';zeros(banddelayR_D(3),1)];
    

 outputL = band1delayL+band2delayL+band3delayL;
 outputR = band1delayR+band2delayR+band3delayR;
 %%%%%%%%%%%%%%%%%%%%%%% Create a fade out %%%%%%%%%%%%%%%%%%%%%%
    
decayConstant = 2;

decayLength = length(outputL);
 
for t = 1:decayLength
   
  
    outputL(t) = outputL(t) * exp(-decayConstant * (t/decayLength));
   
end
 
decayLength = length(outputR);
 
for t = 1:decayLength
   
  
    outputR(t) = outputR(t) * exp(-decayConstant * (t/decayLength));
   
end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%% Sum L and R %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 output = [outputL,outputR];

 output(:,1) = output(:,1)./max(abs(output(:,1)));  %Normalize output
output(:,2) = output(:,2)./max(abs(output(:,2)));

elseif stereo_input_test(:,2)==1

    %% Process Mono channel through Multiband delay
    input = x;
    
    band1 = filter(Hd1, input); %Run input throught LP filter
    band2 = filter(Hd2, input); %Run input throught BP filter
    band3 = filter(Hd3, input); %Run input throught HP filter

    band1delay = my_delay(band1, Fb1, Delay_time1, fs);     %Run LP signal through delay
    band2delay = my_delay(band2, Fb2, Delay_time2, fs);      %Run BP signal through delay
    band3delay = my_delay(band3, Fb3, Delay_time3, fs);      %Run HP signal through delay
    
    % obtain the length of each delayed band 
    banddelay_L(1) = length(band1delay);
    banddelay_L(2) = length(band2delay);
    banddelay_L(3) = length(band3delay);
    
    % calculate the length difference between each delayed band
    banddelay_max = max(banddelay_L);
    banddelay_D(1) = banddelay_max-banddelay_L(1);
    banddelay_D(2) = banddelay_max-banddelay_L(2);
    banddelay_D(3) = banddelay_max-banddelay_L(3);
    
    % padd with zeros tomatch the length of all delayed bands
    band1delay = [band1delay';zeros(banddelay_D(1),1)];
    band2delay = [band2delay';zeros(banddelay_D(2),1)];
    band3delay = [band3delay';zeros(banddelay_D(3),1)];
    
    output = band1delay+band2delay+band3delay; %Recombine delayed signals
   
    %%%%%%%%%%%%%%%%%%%%%%% Create a fade out %%%%%%%%%%%%%%%%%%%%%%
    
decayConstant = 2;

decayLength = length(output);
 
for t = 1:decayLength
   
  
    output(t) = output(t) * exp(-decayConstant * (t/decayLength));
   
end
                 %%%%%%%%%%%%%%%%%%Normalise%%%%%%%%%%%%%%%%%%%%%%%
    output = output./max(abs(output)); % normalisation for processing one-channel (mono) signal
    
end

        %% Generate Plot
    
H1=freqz(Hd1);          %Returns magnitude data from LP filter
H2=freqz(Hd2);          %Returns magnitude data from BP filter
H3=freqz(Hd3);          %Returns magnitude data from HP filter
H1mag=abs(H1);          %Convert magnitude data from LP filter from complex form to absolute values 
H2mag=abs(H2);          %Convert magnitude data from BP filter from complex form to absolute values 
H3mag=abs(H3);          %Convert magnitude data from HP filter from complex form to absolute values 
H1dB=20.*log10(H1mag);  %Convert absolute magnitude data to dB for LP filter
H2dB=20.*log10(H2mag);  %Convert absolute magnitude data to dB for BP filter
H3dB=20.*log10(H3mag);  %Convert absolute magnitude data to dB for HP 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 LP 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('LP','BP','HP','Transfer Function');
set(legendPlot, 'Location', 'NorthEastOutside');
set(outfig, 'Position', [0 0 1200 800]);
    
        %% Play Sound
sound(output,fs)        %Play output sound

audiowrite('outputdelay.wav',output,fs)