function joint_matrix = computeJointMatrix(file_path,opt)
% Read waveform and truncate it to desired length
[original_waveform,original_sample_rate] = wavread(file_path);
if isfield(opt,'waveform_length')
	waveform = original_waveform(1:opt.waveform_length);
else
	waveform = original_waveform;
end

% 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.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 jointgram corresponding to desired j2
index_array = find(averaged_jointgram.meta.j(2,:)==opt.j2);
nCells = length(index_array);

% Compute correspondent arrays of indices j1 and fr_j
j1_array = averaged_jointgram.meta.j(1,index_array);
fr_j_array = averaged_jointgram.meta.fr_j(index_array);
max_j1 = max(j1_array);
max_fr_j = max(fr_j_array)-1; % do not show maximum scale across log-frequencies

% Initialize joint matrix
joint_matrix = zeros(1+max_j1,1+2*max_fr_j);

% Compute expected scattering by taking the mean through time
expected_jointgram_cell = cellfun(@mean,averaged_jointgram.signal,'UniformOutput',false);
expected_jointgram = cell2mat(expected_jointgram_cell);

% Fill in joint matrix
for fr_j = -max_fr_j:max_fr_j
        subindex_array = find(fr_j_array==fr_j);
        column = expected_jointgram(index_array(subindex_array));
        joint_matrix(:,1+fr_j+max_fr_j) = interpft(column,1+max_j1);
end
end
