% Add ScatNet to path
addpath('../scatnet-v0.2.1');
addpath_scatnet;

% Set options
directory = '../SCATT_BIRD_50mostE_perfile_april2014_RB_HG/';
opt.averaging_length = [4096 4096];
nOctaves = log2(opt.averaging_length);
opt.nFilters_per_octave = 16;
opt.quality_factor = 16;
opt.compression_exponent = 1;
opt.nFreq_filters_per_octave = 4;
opt.freq_quality_factor = 4;
% full frequency transposition invariance
opt.frequency_transposition_invariance = ...
	opt.quality_factor * nOctaves(2);
opt.odf_hop_length = 512;
opt.fr_oversampling = 1;

% Explore directory
category_index = 1;
subdirectory = sprintf('%d/',category_index);
file_name_array = dir([directory subdirectory '*.wav']);

% Make file path
file_index = 5;
file_name = file_name_array(file_index).name;
file_path = [directory subdirectory file_name];

% Read waveform
[waveform sample_rate] = wavread(file_path);
waveform_length = length(waveform);

% Compute scalogram
scalogram = computeScalogram(waveform,opt);

% Build the 1d filter bank along the time axis
time_filter_options.Q = 1; % number of filters per octave
time_filter_options.B = 1; % reciprocal octave bandwidth
time_filter_options.J = T_to_J(opt.averaging_length(2),time_filter_options);
time_filters = filter_bank(waveform_length,time_filter_options);

% Build the 1d filter bank along the log-frequency axis
N_freq = 2^nextpow2(numel(scalogram.signal));
freq_filter_options.Q = opt.nFreq_filters_per_octave; % number of filters per octave
freq_filter_options.B = opt.freq_quality_factor; % reciprocal octave bandwidth
freq_filter_options.J = ...
	T_to_J(opt.frequency_transposition_invariance,freq_filter_options);
freq_filters = filter_bank(N_freq,freq_filter_options);
time_frequency_filter_cell = {freq_filters{1},time_filters{1}};

% Perform joint wavelet transform of the scalogram
joint_scattering_options.oversampling = log2(opt.averaging_length(2)/opt.odf_hop_length);
joint_scattering_options.fr_oversampling = opt.fr_oversampling;
[~,joint_wavelet_transform] = ...
        joint_tf_wavelet_layer_1d(scalogram,time_frequency_filter_cell,joint_scattering_options);

% Take element-wise modulus
jointgram = modulus_layer(joint_wavelet_transform);

% Average jointgram over acoustic frequencies
averaged_jointgram = ...
        joint_tf_wavelet_layer_1d(jointgram,time_frequency_filter_cell,joint_scattering_options);

% Get indices in the averaged jointgram
j2_array = averaged_jointgram.meta.j(2,:);
fr_j_array = averaged_jointgram.meta.fr_j;

