function [out_wave] = hammondorganism(input, vibrato, tremolo, freq, depth, output)
% Single Celled Hammond Organism Digital Audio Effect
% Developed by: Hayden P. Tweedley
% Final Revision Date: 08/05/2012
%
% Code created for educational purposes by Hayden P. Tweedley
% extended from code written by S. Disch published in DAFX book 2nd edition, 
% pp.75 - 76, copyright Wiley & Sons 2011, available at http://www.dafx.de.
% It may be used for educational purposes and not for commercial 
% applications without further permission.
%
% Description:
%
%  hammondorganism creates a Leslie-rotary-loudspeaker-like effect by
%  simulating some of the fundamental constituent aspects of the Leslie
%  sound.  These currently include sine modulated vibrato and tremolo.
%
% Input Arguments:
%
% * input: specifies input audio filename of the form 'input.wav'.  Inverted
%   commas must be included and file must be of wav filetype.  File must be
%   included in the same directory from which the hammondorganism function is
%   called.
% 
% * vibrato: specifies depth of the vibrato constituent effect.  Argument is
%   of the range 0 to 1 where 0 indicates no vibrato effect and 1 indicates
%   maximum vibrato effect.  For best results, a value of 0.75 or less is
%   recommended.  0.4 is close to realistic.
%
% * tremolo: specifies depth of the tremolo constituent effect.  Argument is
%   of the range 0 to 1 where 0 indicates no tremolo effect and 1 indicates
%   maximum tremolo effect.  For best results, a value of 0.5 or less is
%   recommended.  0.2 is close to realistic.
% 
% * freq: specifies the modulation frequency governing both vibrato and
%   tremolo in Hz.  For best results, a value of between 1 and 10 is 
%   recommended.
%   This spans the modulation frequency between the slow and fast speeds of a
%   Leslie rotary loudspeaker.  A value of 0 is equivalent to no vibrato 
%   effect and no tremolo effect.  6.5 Hz is close to realistic for fast
%   ("tremolo") Leslie, 3Hz is close to realistic for slow ("chorale") Leslie.
% 
% * depth: specifies delay in seconds governing pitch depth of vibrato.  
%   This delay is responsible for the range of pitch variation heard in the 
%   vibrato effect. For best results, a value between 0.0005 and 0.001 is 
%   recommended.  0.0009 is close to realistic.
%
% * output: specifies output audio filename of the form 'output.wav'.  
%   Inverted commas must be included and file must be of wav filetype 
%   IE. filename must include .wav extension.
%
% Output Arguments:
%
% * Function returns vector out_wave containing the effected audio data.  
%   This is duplicated in the output file as specified by the input 
%   argument: "output".  Call function by:
%   "x = hammondorganism(input, vibrato, tremolo, freq, flanging, output)"
%   (without quotation marks) to write returned out_wave to new vector x,
%   else out_wave is written to default vector ans.
%
% Other Outputs:
%
% * An output file is generated, the filename of which is specified by the 
%   input argument: "output".
%
% * A graph of the input waveform is generated in figure 1
%
% * A graph of the output waveform is generated in figure 2
%
%
% END HELP PAGE

[in_wave, Fs, NBITS ] = wavread(input); % read input waveform, 
% and write: data to variable in_wave, samplingfrequency to variable Fs, 
% number of bits to variable NBITS.

% variable definitions for ease of debug
vib = vibrato; % percent influence of vibrato - try 0.5 or less
trem = tremolo; % percent influence of tremolo - try 0.5 or less
Modfreq = freq; % modulation frequency specified in Hz - try 10Hz or less
Delay = depth;  % basic delay of input sample in sec - try 0.0015 or less

Width = Delay;  % Set width equal to delay (cannot be greater than delay)
% for this application, else width just scales delay, IE. same as depth 
% input argument except that a minimum basic delay can be hard-coded and the
% width variable can then increase upon this.

% DELAY = Basic delay in # of samples
DELAY=round(Delay*Fs); % = basic delay of input sample in seconds, 
% mutliplied by sample rate, rounded to nearest integer

% Modulation width in # of samples (same as DELAY)
WIDTH=round(Width*Fs);

% Error Checking:
if WIDTH>DELAY 
  error('delay greater than basic delay!');
  return;
end

% Modulation frequency in expressed as sample frequency (1/samples) 
MODFREQ=Modfreq/Fs;  %%IE. frequency is 1/time  and FS is sample/time, 
% therefore modfreq/fs equals time/sample * 1/time = 1/sample

% Number of samples in WAV-file
LEN=length(in_wave); % Total number of samples in input signal

% Length of the entire delay
L=2+DELAY+WIDTH*2;



% Memory allocation for delay
Delayline=zeros(L,1); % generate a vector of length equal to length of 
% entire delay line

% Memory allocation for output vector
out_wave = zeros(length(in_wave),1); % generate a vector of length equal to
% length of entire input signal 
% (since output signal will be same length as input signal)  



for n=1:(LEN-1)  % for n = 1 to the length of in_wave by 1 sample increments
    

   M=MODFREQ;  % Modulation frequency expressed as sample frequency
   MOD=sin(M*2*pi*n);  % Sinusoidal modulation function of frequency M weighted by n
   ZEIGER=1+DELAY+WIDTH*(1+vib*MOD); %begin a count 'zeiger' of total delay plus 
% total delay modulated by the function MOD weighted by the vibrato intensity coefficient. 
   i=floor(ZEIGER); % floor rounds to nearest integer towards negative infinity
   frac=ZEIGER-i;  % frac is difference between zeiger and zeiger rounded down to nearest integer.
         
   Delayline=[in_wave(n);Delayline(1:L-1)];  % construct a column vector delayline.  
% Add nth sample of in_wave to the start of the colum vector delayline 
% while shifting all existing values down by one row.

%%% Linear Interpolation (not used in this implementation):
% out_wave(n,1)=((1+trem*MOD)*(Delayline(i+1)*frac+Delayline(i)*(1-frac))); % Linear Interpolation.   
% write output audio data into vector sample by sample, indexed by nth sample, IE. n by 1 vector.
% output audio data is interpolated linearly between samples, and multiplied
% in amplitude by the modulation function as weighted by the tremolo
% intensity coefficient.
%%%

%%% 3rd order Spline interpolation (superior to linear interpolation):
out_wave(n,1)=(1+trem*MOD)*(Delayline(i+1)*frac^3/6+Delayline(i)*...
((1+frac)^3-4*frac^3)/6+Delayline(i-1)*((2-frac)^3-4*(1-frac)^3)/6+Delayline(i-2)*(1-frac)^3/6);% Linear Interpolation.   
% write output audio data into vector sample by sample, indexed by nth sample, IE. n by 1 vector.
% output audio data is interpolated by 3rd order spline between samples, and multiplied
% in amplitude by the modulation function as weighted by the tremolo
% intensity coefficient.
%%%
     
end 

out_wave = out_wave./max(abs(out_wave)); %normalising output audio data     

wavwrite(out_wave, Fs, NBITS,output); % write audio output data to output file
% as specified by input parameter.

%plot input and outputs waveforms respectively
figure(1)
plot(in_wave,'r');
title('Input Waveform');
figure(2)
plot(out_wave,'b');
title('Output Waveform');
%%% END
