function y=Spatial_chorus(x,SAMPLERATE, WidthL, WidthR, DelayL, DelayR, ModRate, ModCorr, ModType, ModFilterF, ModFilterQ, Interpolation, mix, fx_shift, output_level)

% STEREO SPATIAL CHORUS with Fractal modulation. 

% Rob Taylor Student# 311123392
% Lab report 1
% DAS - Semester 1 2012


% Parameters:
% Width.. value between (less width) 0 and (more width) 1
% Delay.. in milliseconds
% ModRate.. > 1 (noise resolution)
% ModCorr..  (less correlated) 0 and (more correlated) 1
% ModType... 1 = fractal, 0 = white noise
% NoiseFilter.. filter frequency
% ModFilterQ...
% Interpolation... 1 = Spline, 2 = All-pass
% mix... mix level, 0 for equal amplitude, -1 for dry, 1 for wet
% fx_shift...  shifts the modulated channel path in milliseconds
% output_level... amplitude of output between 0 and 1

%%
% Default parameters:

if nargin < 15, output_level = 0.999; end
if nargin < 14, fx_shift = 0; end
if nargin < 13, mix = 0; end
if nargin < 12, Interpolation = 1; end
if nargin < 11, ModFilterQ = 2; end
if nargin < 10, ModFilterF = 8; end
if nargin < 9, ModType = 1; end
if nargin < 8, ModCorr = 0; end
if nargin < 7, ModRate = 10; end
if nargin < 6, DelayR = 12; end
if nargin < 5, DelayL = 8; end
if nargin < 4, WidthR = 0.09; end
if nargin < 3, WidthL = 0.05; end
if nargin < 2, SAMPLERATE = 44100; end

%%
% Correlation warning for incorrect paramter

if ModCorr > 1 || ModCorr < 0;
    disp('Oops! ModCorr must be less than 1 or greater than zero');
    ModCorr = 0.5;
end

%%
% Chorus delay and width

DELAYL=round((DelayL./1000)*SAMPLERATE); % in ms
DELAYR=round((DelayR./1000)*SAMPLERATE);

DELAYMAX = max([DELAYL, DELAYR]);

WidthMax = max([WidthL, WidthR]);

LEN=length(x); % # of samples in WAV-file

noisesamples = round(ModRate.*(LEN/SAMPLERATE));

%%
% Modulation type: Case 1 = fractal modulation, Case 2 = white noise modulation

switch ModType
    case 1
        FracCycles = ceil(noisesamples./64);
        noiseL = frand64(FracCycles);
        noiseL = noiseL';
        noiseL = noiseL./max(abs(noiseL));
        noiseR = frand64(FracCycles);
        noiseR = noiseR';
        noiseR = noiseR./max(abs(noiseR));
    case 0
        noiseL = randsig(noisesamples);
        noiseR = randsig(noisesamples);

end

%%
% Modulation filter

pole_coefsL = M_fq2coef(ModFilterF,ModFilterQ,ModRate);
pole_coefsR = M_fq2coef(ModFilterF,ModFilterQ,ModRate);
pole_sumL = sum(pole_coefsL);
pole_sumR = sum(pole_coefsR);
z_coefsL = pole_sumL;
z_coefsR = pole_sumR;


noisefilteredL = filter(z_coefsL, pole_coefsL, noiseL);
noisefilteredR = filter(z_coefsR, pole_coefsR, noiseR);

denominator = round(LEN./noisesamples);

resampled_noiseL = resample(noisefilteredL, denominator, 1);
resampled_noiseR = resample(noisefilteredR, denominator, 1);


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

%%
% Modulation signal normalisation

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

%%
% Memory allocation for delay

L = ceil(6 + DELAYMAX + (DELAYMAX*WidthMax)); % length of the entire delay  
DelaylineL=zeros(L,1); % memory allocation for delay
DelaylineR=zeros(L,1);
yL=zeros(length(x),1); % memory allocation for output vector
yR=zeros(length(x),1);

[resampled_noiseL_mixed, resampled_noiseR_mixed ] = pcor_sigs( resampled_noiseL, resampled_noiseR, ModCorr );


%%
% Input channel number select and process

channels = size(x);
channels = channels(2);

% Mono input path -----------------------

if channels == 1;

    for n=1:(LEN-1)
   
       % Delay line (mono path)   
       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=[x(n);DelaylineL(1:L-1)]; 
       DelaylineR=[x(n);DelaylineR(1:L-1)]; 

       % Interpolation (mono path)
       switch Interpolation
    
        case 1 %---Spline Interpolation-------------------------------
        if iL - 2 <= 0;
            iL = 3;
        end
        
        if iR - 2 <= 0;
            iR = 3;
        end
        
        yL(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; 
        yR(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;


        case 2 %---Allpass Interpolation------------------------------
        yL(n,1)=(DelaylineL(iL+1)+(1-fracL)*DelaylineL(iL)-(1-fracL)*ya_alt);
        yR(n,1)=(DelaylineR(iR+1)+(1-fracR)*DelaylineR(iR)-(1-fracR)*ya_alt);


       end         
           
    end
    
    % Wet/dry mix (mono path)
    mixnorm = mix + 1;
    mixnorm = mixnorm./2;

    drylevel = 1 - mixnorm;
    wetlevel = mixnorm;
    
     x_gain_corrected = x .* drylevel;
     yL = yL .* wetlevel;
     yR = yR .* wetlevel;
    
    % Modulated channel path shift in milliseconds(mono path)
    if fx_shift == 0;
        yL = yL;
        yR = yR;
        x_gain_corrected = x_gain_corrected;
    elseif fx_shift > 0;
        delay_samples = zeros(round(((fx_shift./1000).*SAMPLERATE)),1);
        yL = [delay_samples; yL];
        yR = [delay_samples; yR];
        x_gain_corrected = [x_gain_corrected; delay_samples];
    elseif fx_shift <0;
        delay_samples = zeros(round(((fx_shift./1000).*SAMPLERATE)),1);
        yL = [yL; delay_samples];
        yR = [yR; delay_samples];
        x_gain_corrected = [delay_samples; x_gain_corrected];
    end
        
    % Wet/dry signal mixer
    yL = yL + x_gain_corrected;
    yR = yR + x_gain_corrected;

% Stereo input path ------------------   

elseif channels == 2;
    
       for n=1:(LEN-1)

       % Delay line (stereo path)    
       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=[x(n,1);DelaylineL(1:L-1)]; 
       DelaylineR=[x(n,2);DelaylineR(1:L-1)]; 

       % Interpolation (stereo path)
       switch Interpolation

        case 1 %---Spline Interpolation-------------------------------
        if iL - 2 <= 0;
            iL = 3;
        end
        
        if iR - 2 <= 0;
            iR = 3;
        end
        yL(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; 
        yR(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;


        case 2 %---Allpass Interpolation------------------------------
        yL(n,1)=(DelaylineL(iL+1)+(1-fracL)*DelaylineL(iL)-(1-fracL)*ya_alt);
        yR(n,1)=(DelaylineR(iR+1)+(1-fracR)*DelaylineR(iR)-(1-fracR)*ya_alt);


       end
           
        end
        
    % Wet/dry mix (stereo path) 
    mixnorm = mix + 1;
    mixnorm = mixnorm./2;

    drylevel = 1 - mixnorm;
    wetlevel = mixnorm;
        
    x_gain_corrected = x .* drylevel;
    yL = yL .* wetlevel;
    yR = yR .* wetlevel;
     
    % Modulated channel path shift in millisecond (stereo path)
     if fx_shift == 0;
        yL = yL;
        yR = yR;
        x_gain_corrected = x_gain_corrected;
    elseif fx_shift > 0;
        delay_samples_mono = zeros(round(((fx_shift./1000).*SAMPLERATE)),1);
        delay_samples_stereo = zeros(round(((fx_shift./1000).*SAMPLERATE)),2);
        yL = [delay_samples_mono; yL];
        yR = [delay_samples_mono; yR];
        x_gain_corrected = [x_gain_corrected; delay_samples_stereo];
    elseif fx_shift <0;
        delay_samples_mono = zeros(round(((fx_shift./1000).*SAMPLERATE)),1);
        delay_samples_stereo = zeros(round(((fx_shift./1000).*SAMPLERATE)),2);
        yL = [yL; delay_samples_mono];
        yR = [yR; delay_samples_mono];
        x_gain_corrected = [delay_samples; x_gain_corrected];
    end
        
    % Wet/dry signal mixer (stereo path)
    yL = yL + x_gain_corrected(:,1);
    yR = yR + x_gain_corrected(:,2);     
    
end  

%%
% Two-channel matrix

y = [yL yR];

%%
% Output normalisation

y = (y./max(max(abs(y)))).*output_level;


