function [output] = SpatialChorus(input, fs, WidthL, WidthR, DelayL, DelayR, ModRate, ModCorr, ModType, ModFilterF, ModFilterQ, Mix, WetShift, DryShift, OutputLevel)
%
% SPATIAL CHORUS WITH FRACTAL MODULATION AND EXTENDED USER CONTROLS
%
% Robert W Taylor - 311123392(SID)
% Digital Audio Systems, DESC9115, Semester 1 2012
% Graduate Program in Audio and Acoustics
% Faculty of Architecture, Design and Planning, The University of Sydney
%   
% This function represents a new design for a spatial chorus effect that 
% incorporates a modulation signal derived from a fractal algorithm with 
% low-pass filter control, and the implementation of extensive user controls 
% and parameters that center around the de-correlation of left/right signal 
% paths resulting in a spatially enhanced outcome.
%
% SpatialChorus(input, fs, WidthL, WidthR, DelayL, DelayR, ModRate, ModCorr, 
% ModType, ModFilterF, ModFilterQ, Mix, WetShift, DryShift, OutputLevel)
%
% ***Parameters and suggested working ranges***
% 
% fs = Sample rate of sound source in Hz.
% 
% DepthL, DepthR = Peak deviation from base delay measured as a ratio of 
%   the maximum amount of delay imposed. 
%   Value between 0(less width) and 1(more width). 
%   Independent controls for left/right output signals.
% 
% DelayL, DelayR = Delay offset of modulating signal in milliseconds. 
%   Independent controls for left/right output signals.
% 
% ModRate = Rate of modulation in Hz. (Suggested: 1 - 30 Hz)
% 
% ModCorr = Correlation of left and right modulation signals. 
%   Value between 0 (no correlation) and 1 (correlation). 
% 
% ModType = Select the type of modulation signal 
%   0 = Fractal, 1 = White noise, 2 = Brown noise, 3 = Pink noise. 
%   (Suggested: 0 = Fractal)
% 
% ModFilterF = Modulation signal filter frequency selection in Hz. 
%   (Suggested: 1 - 30 Hz)
% 
% ModFilterQ = Modulation signal filter bandwidth (Q) in Hz. 
%   (Suggested: 1 - 5 Hz)
% 
% Mix = Mix level between dry/wet signals. 
%   Values between -1 (100% dry), 0 (50/50), 1 (100% wet).  
%   (Suggested: 0 (50/50))
% 
% WetShift = Delay of the modulated channel path in milliseconds.
% 
% DryShift = Delay of the un-modulated channel path in milliseconds.
% 
% OutputLevel = Overall amplitude of output signal.  
%   Range between 0 and 1.
%
% ***Diagnostics***
%
%       The primary aim of the function is to create spatial effects from 
%       a mono source and therefore this function only accepts mono inputs.
%
%       Four additional functions (Spatial Chorus toolkit) provide key 
%       components to the working effect and are required to be present 
%       in the MatLab environment for it to function.
%
%               - Fractal modulator [frand64.m] : This function implements 
%                 a 64-sample fractal algorithm to create a stochastic, 
%                 statistically self-similar pseudo-random noise signal. 
%                 (Created by Assoc. Professor William L. Martens) 
%
%               - Generation of filter coefficients [M_fq2coef.m] :  
%                 This function converts the frequency and bandwidth parameters 
%                 to coefficients. It returns a standard filter vector (a). 
%                 (Created by Luis Miranda)
%
%               - Independent left/right modulation signal mixing [pcor_sigs.m]: 
%                 This function controls the degree of correlation between 
%                 the left and right independent modulator signals. 
%                 (Created by Assoc. Professor William L. Martens and Luis Miranda)
%
%               - Noise generator [makenoise.m] : This general purpose function 
%                 creates various types of noise signals as alternatives 
%                 to the fractal version for testing and evaluation.  
%                 (Created by Dr Densil Cabrera)
%
% ***Tips***
%
%       Although the parameters supply a large range and quite severe sounding 
%       results are possible, it is advisable to stay within specific value 
%       ranges to achieve more standard results. Generally a modulation rate 
%       of 5-10ms, a depth of between 0.02 and .1, and a delay offset 
%       of between 10-15ms will provide a pleasing result for most applications.
%
%       To obtain spatial [width] it is advised to have the modulators 
%       de-correlated (0), and slight variations between left and right 
%       depth and delay parameters
%
%       To obtain spatial [movement] it is advised to delay the un-modulated 
%       (dry) signal at the same rate as the delay offset (5-10ms).
%
%       Combining the two spatial attributes together will achieve a 
%       significant increase in spatial perception as a whole.
%
%
%       For more detailed information seek 'Help on selection'


%--Delay Width Length------------------------

DelayL=round((DelayL./1000)*fs); % in ms
DelayR=round((DelayR./1000)*fs);

DelayMax = max([DelayL, DelayR]);

WidthMax = max([WidthL, WidthR]);

Length=length(input); % # of samples in WAV-file

noisesamples = round(ModRate.*(Length/fs));

microsModL = ((DelayL * WidthL) + DelayL ) ./ 10;
microsModR = ((DelayR * WidthR) + DelayR ) ./ 10;


%--Modulation type switch---------------------

switch ModType
    case 0
        FracCycles = ceil(noisesamples./64);
        noiseL = frand64(FracCycles);
        noiseL = noiseL';
        noiseL = noiseL./max(abs(noiseL));
        noiseR = frand64(FracCycles);
        noiseR = noiseR';
        noiseR = noiseR./max(abs(noiseR));
        %'frand64' function
    case 1
        noiseL = randsig(noisesamples);
        noiseR = randsig(noisesamples);
        %'randsig' function
    case 2
        noise_duration = noisesamples./ModRate;
        noiseL = makenoise(-2, noise_duration, ModRate, 0, ModRate, 1, 0, 0);
        noiseR = makenoise(-2, noise_duration, ModRate, 0, ModRate, 1, 0, 0);
    case 3
        noise_duration = noisesamples./ModRate;
        noiseL = makenoise(-1, noise_duration, ModRate, 0, ModRate, 1, 0, 0);
        noiseR = makenoise(-1, noise_duration, ModRate, 0, ModRate, 1, 0, 0);
        %'makenoise' function

end

%--Resampling---------------------------------

denominator = round(Length./noisesamples);

resampled_noiseL = resample(noiseL, denominator, 1);
resampled_noiseR = resample(noiseR, denominator, 1);


if  length(resampled_noiseL) < Length;
    endzeros = zeros(Length-length(resampled_noiseL),1);
    resampled_noiseL = cat(1, resampled_noiseL, endzeros);  
    resampled_noiseR = cat(1, resampled_noiseR, endzeros);
end


%--Low pass filter----------------------------

pole_coefsL = M_fq2coef(ModFilterF,ModFilterQ,fs);
pole_coefsR = M_fq2coef(ModFilterF,ModFilterQ,fs);
% 'M_fq2coef' function

pole_sumL = sum(pole_coefsL);
pole_sumR = sum(pole_coefsR);

z_coefsL = pole_sumL;
z_coefsR = pole_sumR;

noisefilteredL = filter(z_coefsL, pole_coefsL, resampled_noiseL);
noisefilteredR = filter(z_coefsR, pole_coefsR, resampled_noiseR);



%--Normalise modulation signal----------------

resampled_noiseL = noisefilteredL./max(abs(noisefilteredL));
resampled_noiseR = noisefilteredR./max(abs(noisefilteredR));


%--Memory allocation--------------------------

L = ceil(6 + DelayMax + (DelayMax*WidthMax)); 
% Length of the entire delay  

DelayLineL=zeros(L,1); 
DelayLineR=zeros(L,1);
% memory allocation for delay

outputL=zeros(length(input),1); 
outputR=zeros(length(input),1);
% memory allocation for output vector


%--Modulation signal correlation mixing-------

[resampled_noiseL_mixed, resampled_noiseR_mixed ]...
= pcor_sigs( resampled_noiseL, resampled_noiseR, ModCorr );
% 'pcor_sigs' function


% Delay line---------------------------------- 

    for n=1:(Length-1)
    %    M=MODFREQ;
    %    MOD=sin(M*2*pi*n);

       MODL = resampled_noiseL_mixed(n,1); 
       MODR = resampled_noiseR_mixed(n,1);


       TAPL=1+DelayL+(DelayL*WidthL*MODL);
       TAPR=1+DelayR+(DelayR*WidthR*MODR);


       iL=floor(TAPL);
       iR=floor(TAPR);

       fracL=TAPL-iL;
       fracR=TAPR-iR;


       DelayLineL=[input(n);DelayLineL(1:L-1)]; 
       DelayLineR=[input(n);DelayLineR(1:L-1)]; 
 
%--Spline Interpolation-----------------------

    if iL - 2 <= 0;
       iL = 3;
    end
        
    if iR - 2 <= 0;
       iR = 3;
    end
        
outputL(n,1)=DelayLineL(iL+1)*fracL^3/6+DelayLineL(iL)*((1+fracL)^3-4*fracL^3)...
    /6+DelayLineL(iL-1)*((2-fracL)^3-4*(1-fracL)^3)/6+DelayLineL(iL-2)*(1-fracL)^3/6; 
outputR(n,1)=DelayLineR(iR+1)*fracR^3/6+DelayLineR(iR)*((1+fracR)^3-4*fracR^3)...
    /6+DelayLineR(iR-1)*((2-fracR)^3-4*(1-fracR)^3)/6+DelayLineR(iR-2)*(1-fracR)^3/6;

    end


%--Dry signal shift--------------------------       
           
    if DryShift == 0;
       outputL = outputL;
       outputR = outputR;
       input = input;
    elseif DryShift > 0;
       delay_samples = zeros(round(((DryShift./1000).*fs)),1);
       outputL = [outputL; delay_samples];
       outputR = [outputR; delay_samples];
       input = [delay_samples; input];
    elseif DryShift <0;
       delay_samples = zeros(round(((DryShift./1000).*fs)),1);
       outputL = [delay_samples; outputL];
       outputR = [delay_samples; outputR];
       input = [input; delay_samples];
    end
    
%--Wet signal shift--------------------------       
           
    if WetShift == 0;
       outputL = outputL;
       outputR = outputR;
       input = input;
    elseif WetShift > 0;
       delay_samples = zeros(round(((WetShift./1000).*fs)),1);
       outputL = [delay_samples; outputL];
       outputR = [delay_samples; outputR];
       input = [input; delay_samples];
    elseif WetShift <0;
       delay_samples = zeros(round(((WetShift./1000).*fs)),1);
       outputL = [outputL; delay_samples];
       outputR = [outputR; delay_samples];
       input = [delay_samples; input];
    end
 

%-- Dry/wet mix------------------------------  
    
mixnorm = Mix + 1;
mixnorm = mixnorm./2;

drylevel = 1 - mixnorm;
wetlevel = mixnorm;  
     
input_gain_corrected = input .* drylevel;
outputL = outputL .* wetlevel;
outputR = outputR .* wetlevel;       
    
outputL = outputL + input_gain_corrected;
outputR = outputR + input_gain_corrected;


%--Two channel matrix------------------------

output = [outputL outputR];

%--Normalise and final gain to output--------

output = (output./max(max(abs(output)))).*OutputLevel;

