function [ audio_out ] = autowah4( audio_in, fs )
% By Bella
% Applies an enveloped autowah function with 2 formants to a given audio file.
% Script inputs:
% - audio_in: audio file for autowah to be applied.
% - fs: Sample rate, given by audio file.
% In dialogue box:
% - min_f1: Minimum cutoff frequency of bandpass for formant 1 (recommend 300).
% - max_f1: Maximum cutoff frequency of bandpass for formant 1 (recommend 1000).
% - min_f2: Minimum cutoff frequency of bandpass for formant 2 (recommend 1500).
% - max_f2: Maximum cutoff frequency of bandpass for formant 2 (recommend 2200).
% - mix: Wet/dry mix 0 to 1; 0=dry, 1=wet.

% Dialogue box for inputs
prompt = {'Oscillation rate (Hz)','Minimum cutoff frequency of bandpass for formant 1 (Hz)','Maximum cutoff frequency of bandpass for formant 1 (Hz)','Minimum cutoff frequency of bandpass for formant 2 (Hz). If second formant is not required, enter 0','Maximum cutoff frequency of bandpass for formant 2 (Hz). If second formant is not required, enter 0','Wet/dry mix from 0 to 1 where 0=dry, 1=wet'};
dlg_title = 'Autowah input';
num_lines = 1;
defaultans = {'1000','300','1000','1500','2200','0.5'};
G = inputdlg(prompt,dlg_title,num_lines,defaultans);

% Read inputs
rate = str2num(G{1});
min_freq1 = str2num(G{2});
max_freq1 = str2num(G{3});
min_freq2 = str2num(G{4});
max_freq2 = str2num(G{5});
mix = str2num(G{6});

% No. channels variable
[~, channels] = size(audio_in);

% Create envelope
env = abs(hilbert(audio_in)); 

% Define fc
centre_freq = rate/fs;

% Damping factor
damping = 0.05;

% FORMANT 1

% Sweep fc between min and max frequency
cutoff_freq1=min_freq1:centre_freq:max_freq1;
while(length(cutoff_freq1) < length(audio_in) )
cutoff_freq1 = [ cutoff_freq1 (max_freq1:-centre_freq:min_freq1) ];
cutoff_freq1 = [ cutoff_freq1 (min_freq1:centre_freq:max_freq1) ];
end

% Set length of modulation signal to match audio input
cutoff_freq1 = cutoff_freq1(1:length(env));

% State variable filters
F1 = 2*sin((pi*cutoff_freq1(1))/fs);
Q1 = 2*damping;

% Create and zero vectors that match length of audio input file
highpass1=zeros(size(env)); 
bandpass1=zeros(size(env));
lowpass1=zeros(size(env));

% Begin processing from first sample
highpass1(1) = env(1);
bandpass1(1) = F1*highpass1(1);
lowpass1(1) = F1*bandpass1(1);

% Processing using the differential equation
for i =1:channels;
    for n=2:length(env);
    highpass1(n,i) = env(n,i) - lowpass1(n-1,i) - Q1*bandpass1(n-1,i);
    bandpass1(n,i) = F1*highpass1(n,i) + bandpass1(n-1,i);
    lowpass1(n,i) = F1*bandpass1(n,i) + lowpass1(n-1, i);
    F1 = 2*sin((pi*cutoff_freq1(n))/fs);
    end
end

% FORMANT 2

% Sweep fc between min and max frequency
cutoff_freq2=min_freq2:centre_freq:max_freq2;
while(length(cutoff_freq2) < length(audio_in) )
cutoff_freq2 = [ cutoff_freq2 (max_freq2:-centre_freq:min_freq2) ];
cutoff_freq2 = [ cutoff_freq2 (min_freq2:centre_freq:max_freq2) ];
end

% Set length of modulation signal to match audio input
cutoff_freq2 = cutoff_freq2(1:length(env));

% State variable filters
F1 = 2*sin((pi*cutoff_freq2(1))/fs);
Q1 = 2*damping;

% Create and zero vectors that match length of audio input file
highpass2=zeros(size(env)); 
bandpass2=zeros(size(env));
lowpass2=zeros(size(env));

% Begin processing from first sample
highpass2(1) = env(1);
bandpass2(1) = F1*highpass1(1);
lowpass2(1) = F1*bandpass1(1);

% Processing using the differential equation
for j =1:channels;
    for n=2:length(env);
    highpass2(n,j) = env(n,j) - lowpass2(n-1,j) - Q1*bandpass2(n-1,j);
    bandpass2(n,j) = F1*highpass2(n,j) + bandpass2(n-1,j);
    lowpass2(n,j) = F1*bandpass2(n,j) + lowpass2(n-1, j);
    F1 = 2*sin((pi*cutoff_freq2(n))/fs);
    end
end

% Combine input with envelope and multiply by wet/dry mix
wah = audio_in + bandpass1.*mix + bandpass2.*mix;

% Normalised output
audio_out = wah./max(max(abs(wah))); 

end
