function [lower_freq, upper_freq, center_freq, filterbank] = make_filterbank(num_chan, filt_ord, low_freq, high_freq, nyquist_freq)

%% AAB and FG 10/18
% Function to create log-spaced filterbank
%
% Inputs]
% num_chan - number of vocoder channels
% filt_ord - filter order
% low_freq - lowest channel frequency lower cutoff (Hz)
% high_freq - highest channel frequency upper cutoff (Hz)
% nyquist_freq - highest frequency that can be represented at a given
% sample rate, sample_rate/2
%
% Outputs
% upper_freq(ii) - upper frequency cutoff of channel ii (Hz)
% lower_freq(ii) - lower frequency cutoff of channel ii (Hz)
% center_freq(ii) - centre frequency of channel ii (Hz)
% filter_a(jj,:) - analysis/synthesis filter coefficients 'a'
% filter_b(jj,:) - analysis/synthesis filter coefficients 'a' 
%
%e.g. [upper_freq, lower_freq, center_freq, filter_a, filter_b] = make_filterbank(16, 6, 250, 8000, 22050);

log_freq_range = log10(high_freq / low_freq); % Frequency range in logarithmic units
freq_interval = log_freq_range / num_chan; % Logarithmic frequency interval per channel

lower_freq = zeros (1, num_chan); % Make a zero vector for the center frequencies of each channel
upper_freq = lower_freq; % Make a zero vector for the upper cutoff frequencies of each channel
center_freq = lower_freq; % Make a zero vector for the lower cutoff frequencies of each channel 

for ii = 1:num_chan
    
    upper_freq(ii) = low_freq*10^(freq_interval*ii); % Upper frequency cutoff of channel ii (Hz)
    lower_freq(ii) = low_freq*10^(freq_interval*(ii-1)); % Lower frequency cutoff of channel ii (Hz)
    center_freq(ii) = 0.5*(lower_freq(ii)+upper_freq(ii)); % Centre frequency of channel ii (Hz)
    
end

% Make zero vectors for analysis and synthesis filterbank filter coefficients
filterbank = zeros(num_chan, 2, filt_ord+1); 

for jj=1:num_chan
    
    w = [lower_freq(jj)/nyquist_freq, upper_freq(jj)/nyquist_freq];
    [b,a] = butter(filt_ord/2, w); %  Obtain Butterworth filter coefficients. Not using a highpass only filter for top channel. Check order is correct for bandpass
    
    filterbank(jj, 1, 1:filt_ord+1) = a;   % Save the analysis/synthesis filter coefficients 'a'
    filterbank(jj, 2, 1:filt_ord+1) = b;   % Save the analysis/synthesis filter coefficients 'b'
    
end