function output=VIBRATOOL(Input,FS,Width,frequency,q,type,fade)
% function output=VIBRATOOL(Input,FS,Width,frequency,q,type,fade)
%        
%       This function allows you to apply a vibrato (delay modulation)
%       effect to an input wave. The function gives you the option of 5
%       different types of modulation specified by 'type'. This function
%       was created with violin in mind and as such includes the functions
%       'proclpc' and 'synlpc'. These functions seperate the source
%       (excitation sound) from the filter (resonant structure) of the
%       input wave so that vibrato is only applied to the source sound
%       whilst the resonance of the input is maintained. The default values
%       for when no input arguments are given are also set with the violin in mind.
%       The use of filtered fractal noise as the modulation wave is chosen
%       as default as it is quasi-periodic and hopefully better matches a more 'human'
%       vibrato than the other modulation types. The other types of
%       modulation are provided to allow the user to experiment.
%
% Syntax:
% output=VIBRATOOL(Input,FS,Width,frequency,q,type)
%           
%       Input       = Input wave to apply vibrato effect to (.wav file)
%       FS          = Sampling Frequency (generally 44100 or 48000)
%       Width       = Modulation Width in seconds (0.002-0.006)
%       frequency   = Modulation frequency, AKA rate (2-8)
%       q           = Width of filter used in 'type' cases 0 and 1 (15-50)
%       type        = Type of modulation (0=Fractal noise, 1=White noise,
%                     2=sine wave, 3=square wave, 4=sawtooth wave
%       fade        = Vibrato fade in (0=no fade in, 1=1 second fade in,
%                     2=2 second fade in, 3=3second fade in)
%
% Signal Flow:
%
% Input --> Source/Filter Separation --> Modulation wave created -->
% --> Vibrato applied to source using modulation wave --> Source/Filter resynthesis --> Output as
% variable 'output'.
%
%                   Example of default values:
%           output=VIBRATOOL(Input,44100,0.003,5,30,0,0)
%
%
%       This function was adapted from 'vibrato.m' by Udo Zoelzer, ed. "DAFX: Digital
%       Audio Effects". Wiley, 2002, pp. 68-69. Additional code was adapted
%       from 'get_LP_filtered_mod.m' by William Martens and
%       'generatesounddemo.m' by Densil Cabrera.

if nargin <7, fade = 0; end         % default no fade in of vibrato
if nargin <6, type = 0; end         % default modulation type is fractal noise filtered to desired frequency
if nargin <5, q = 30; end           % default q for filter is 30 - higher values create a narrower filter, only applies to type 0 and 1
if nargin <4, frequency = 5; end    % default modulation frequency is 5Hz
if nargin <3, Width = 0.003; end    % default Width is 0.004
if nargin <2, FS = 44100; end     % default sampling frequency is 44100


[aCoeff,resid,pitch,G,parcor,stream] = proclpc(Input,FS); % Separate source and filter of input wave
wavwrite(stream,FS,'stream.wav'); % Write source to disk as a wave file named 'stream.wav'
[stream,FS]=wavread('stream.wav'); % Import wave file 'stream.wav', testing shows smoother operation by doing these two steps

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

duration=round(LEN/FS*1.1); % length in seconds of 'stream' + 10% (to ensure modulating wave is definitely long enough for input)


% generate a list of times (in seconds) for each sample
samples = 1:round(FS*duration); % sample numbers (1-FS)
times = (samples-1) ./ FS; % times in seconds (time in seconds at each sampling point)

switch type % allows you to switch between different code for for 'type', dependant on the number of cases below
    case 0
        mod = get_LP_filtered_mod(frequency,q,duration,FS); %call function 'get_LP_filtered_mod' to create filtered fractal noise
    case 1
        % generate white noise and then filter it to specified frequency
        mod = rand(1,FS*duration)-0.5; %create variable mod using rand function
        p_coefs=M_fq2coef(frequency,q,FS); % Filter mod to desired frequency and q using function 'M_fq2coef'
        pole_sum = sum(p_coefs); %sum output of function 'M_fq2coef'
        z_coefs=pole_sum; 
        mod = filter(z_coefs, p_coefs, mod); % Use 'filter' function to filter mod with p_coefs and z_coefs
        mod = normsig(mod); %normalise 'mod'
    case 2
        % generate a sine wave
        mod = sin(2*pi*frequency.*times); %create sine wave of specified frequency using built-in 'sin' function
    case 3
        % generate a square wave
        mod = square(2*pi*frequency.*times); %create square wave of specified frequency using built-in 'sin' function
    case 4
        % generate a sawtooth wave
        mod = sawtooth(2*pi*frequency.*times); %create sawtooth wave of specified frequency using built-in 'sin' function
end

switch fade
    case 0
        mod=mod; %no fade in of vibrato
    case 1
        scalingsig1 = linspace(0,1,FS); % create vector of length FS rising from 0 to 1
        scalingsig2 = zeros(1,FS*duration); %create vector of zeros equal to length of stream-FS
        scalingsig2 = scalingsig2+1; % make those zeros into ones
        scalingsignal = [scalingsig1,scalingsig2]; %create vector with fade in in of 1 seconds and then ones
        scalingsignal = .99*normsig(scalingsignal); %normalise scaling signal
        scalingsignal = scalingsignal(1:(length(mod))); %Ensure scalingsignal is equal to mod length
        mod = mod.* scalingsignal; %make modulating sine wave with fade in equal to 1 seconds (FS)
    case 2
        scalingsig1 = linspace(0,1,2*FS); % create vector of length 2*FS rising from 0 to 1
        scalingsig2 = zeros(1,FS*duration); %create vector of zeros equal to length of stream-2*FS
        scalingsig2 = scalingsig2+1; % make those zeros into ones
        scalingsignal = [scalingsig1,scalingsig2]; %create vector with fade in in of 2 seconds and then ones
        scalingsignal = .99*normsig(scalingsignal); %normalise scaling signal
        scalingsignal = scalingsignal(1:(length(mod))); %Ensure scalingsignal is equal to mod length
        mod = mod.* scalingsignal; %make modulating sine wave with fade in equal to 2 seconds (FS)
    case 3
        scalingsig1 = linspace(0,1,3*FS); % create vector of length 3*FS rising from 0 to 1
        scalingsig2 = zeros(1,FS*duration); %create vector of zeros equal to length of stream-3*FS
        scalingsig2 = scalingsig2+1; % make those zeros into ones
        scalingsignal = [scalingsig1,scalingsig2]; %create vector with fade in in of 3 seconds and then ones
        scalingsignal = .99*normsig(scalingsignal); %normalise scaling signal
        scalingsignal = scalingsignal(1:(length(mod))); %Ensure scalingsignal is equal to mod length
        mod = mod.* scalingsignal; %make modulating sine wave with fade in equal to 3 seconds (FS)
end


mod=mod(1:LEN); %makes mod variable same length as stream

ya_alt=0;
% Note that Delay=Width, the basic delay of input sample in sec
DELAY=round(Width*FS); % basic delay in # samples
WIDTH=round(Width*FS); % modulation width in # samples
L=2+DELAY+WIDTH*2;    % length of the entire delay  
Delayline=zeros(L,1); % memory allocation for delay
y=zeros(size(stream));     % memory allocation for output vector

for n=1:(LEN-1) %create vector 'n' with a point for every sampling point of input
   ZEIGER=1+DELAY+WIDTH*mod(n); %use the width input agument and the modulating wave 'mod' to calculate delay at each point 'n'
   i=floor(ZEIGER); % get red of fractional result usinf 'floor' function
   frac=ZEIGER-i; % work out fractional portion we just got rid of
   Delayline=[stream(n);Delayline(1:L-1)]; %delay 'stream' at each point 'n'
   %---Allpass Interpolation------------------------------
   y(n,1)=(Delayline(i+1)+(1-frac)*Delayline(i)-(1-frac)*ya_alt);  %interpolate delay values to sampling frequency
   ya_alt=y(n,1);
end


synWave = synlpc(aCoeff,y,FS,G); %reconstruct vibrato output with Input filter (resonance)
output=synWave; %name 'output'
output = .99*normsig(output); %normalise 'output'



