function output = adaptivetremolo(input,fs,ascending,depth,figs)
%ADAPTIVETREMOLO Applies a time-varying tremolo or amplitude
%   modulation to an input mono audio signal, where the modulation
%   frequency at the output changes with the detected fundamental
%   frequency of the input.
%   
%   Input arguments:
%   input - the input mono audio signal
%   fs - the sample rate of the input audio
%   ascending - specifies how modulation rate varies with
%   fundamental frequency. Set to 1 for ascending or 0 for descending
%   modulation rate with increasing fundamental frequency.
%   depth - depth of the tremolo effect. Set to an integer between
%   0 and 100, where 0 is minimum (no effect) and 100 is maximum depth.
%   figs - set to 1 to generate plots of the detected fundamental
%   frequency and the time-varying amplitude modulator
%        
% Disclaimer: The yinDAFX pitch tracking code by Adrian v.d. Knesebeck
% is sourced from Chapter 9 "Adaptive digital audio effects" of DAFX:
% Digital Audio Effects 2nd Edition Ed. U. Zölzer (2011) p352-353


%% PITCH TRACKING

% initialisation 

% The first three variables were originally input arguments
% for the yinDAFX function, but for this function are
% initialised here instead. The execution timing variables
% tic and toc were also added for testing purposes


x = input;
f0min = 50;  % minimum detectable pitch
hop = 1024;  % hop size
yinTolerance = 0.22;
taumax = round(1/f0min*fs);
yinLen = 1024;
k = 0;

tic % begin recording execution time of the YIN algorithm

% frame processing
for i = 1:hop:(length(x)-(yinLen+taumax))
    k=k+1;
    xframe = x(i:i+(yinLen+taumax));
    yinTemp = zeros(1,taumax);
    % calculate the square differences
        for tau=1:taumax
            for j=1:yinLen
                yinTemp(tau) = yinTemp(tau) + (xframe(j) - xframe(j+tau))^2;
            end
        end
    
    % calculate cumulated normalization
        tmp = 0;
        yinTemp(1) = 1;
        for tau=2:taumax
            tmp = tmp + yinTemp(tau);
            yinTemp(tau) = yinTemp(tau) *(tau/tmp);
        end
    
    % determine lowest pitch
      tau=1; 
      while(tau<taumax)
          if(yinTemp(tau) < yinTolerance)  
              % search turning point
              while (yinTemp(tau+1) < yinTemp(tau)) 
                  tau = tau+1;
              end
              pitch(k) = fs/tau;
              break
          else
              tau = tau+1;
          end
      % if no pitch detected
      pitch(k) = 0;
      end     
end
toc % end recording execution time

%% POST-PROCESSING
% Sample and hold freq values to create a
% frequency vector equal in length to the input signal

hold = repmat(pitch,hop,1); % duplicate the pitch row vector 
% into a matrix with a number - equal to the hop size - of copies 
% of the same row.
f0 = hold(:)';  % merge all the columns into a single column and
% transpose into a row
padsize = length(input) - length(f0); % calculate remaining samples to pad
% if necessary
f0 = [f0 zeros(1,padsize)]; % applies zero padding at the end

%% TREMOLO MAPPING
% Compresses the range of detected fundamental frequencies
% into a range of modulation rates between 1 Hz and 14 Hz,
% either ascending or descending with increasing fundamental
% frequency, depending on the ascending input argument

f_min = min(f0(f0~=0)); % find the min frequency value not equal to zero
f_max = max(f0); % find the max frequency value

switch ascending
    case 1 % modulation rate ascending with fundamental freq
        f_mod = 1 + 13*((f0-f_min)/(f_max-f_min));
    case 0 % modulation rate descending with fundamental freq
        f_mod = 1 + 13*((f_max-f0)/(f_max-f_min));
end

%% TREMOLO OUTPUT
% Create the time-varying amplitude modulation and apply
% it to the output

n = 0:(length(x) - 1); % create samples vector

amp_mod = 10.^(depth*(sin(2*pi*n.*f_mod/fs)-1)/40); % amplitude modulation
c = amp_mod';

output = input .* (1 + c); % tremolo output


%% PLOTTING
% Generates plots of the detected fundamental
% frequency and the time-varying amplitude modulator
% depending on the value of the figs input argument

if figs == 1
    figure(1)
    plot(f0)
    title('Detected Fundamental Frequencies')
    xlabel('Samples')
    ylabel('Frequency (Hz)')
    figure(2)
    plot(c)
    title('Time-Varying Amplitude Modulator')
    xlabel('Samples')
    ylabel('Normalised Amplitude')
end
end

