% Evan Shenkman
% MMI610
% Project 1
% Tone Burst OAEs

% --- Setup simulation params
fs = 48000;                                 % sampling rate
t = 0:1.0/fs:0.05;                          % time vec
f0 = [500 1000 1500 2000 3000 4000 6000];   % the freqs under test vector
minCycles = 5;                              % min number of periods in each toneburst
minSamples = 1/min(f0) * fs * minCycles;    % maintain signal length for equal energy
nsigs = 4*length(f0);                       % number of signals


% --- Model params
spl = zeros(1, nsigs);
irron = zeros(1, nsigs);
normalizeRms = ones(1, nsigs);
subjectNo = 1;
fc = 'all';


% --- Signal Generation
% The input to the model is Nx(4*length(f0)) (i.e. 28 N-point signals)
% There are four signals per frequency
% - Sig 1: a positive, 5 cycle toneburst
% - Sig 2: a positive, 5 cycle toneburst
% - Sig 3: a negative, 5 cycle toneburst
% - Sig 4: a positive, 5 cycle toneburst

% - allocate space for the signal matrix AND the rest of the simulation params
insig = zeros(length(t), nsigs);

% Do the Signal Generation
% - including spl setting for DNLR technique
% - including whether or not irregularities are on
i = 1; % signal number counter
for f = f0
    % Determine the number of cycles for each freq
    cycles = ceil(minSamples * f / fs);
    
    % Generate the 5 cycle toneburst for frequency f
    tone = sin(2*pi*f*t(1:ceil(cycles*fs/f)+1))';
    
    % Set the input signals accordingly
    insig(1:length(tone),i) = tone;
    insig(1:length(tone),i+1) = tone;
    insig(1:length(tone),i+2) = -1 * tone;
    insig(1:length(tone),i+3) = tone;
    
    % Set the SPLs and irregularity flags
    spl(i:i+3) = [50, 50, 50+mag2db(2), 50];
    irron(i:i+3) = [1, 1, 1, 0];
    
    % Move up the counter by 4
    i = i+4;
end


% --- Run the model for ALL signals
[~,~,E,~]=verhulst2012(insig,fs,fc,spl,normalizeRms,subjectNo,irron);


% --- Calculate the OAEs (two different ways)
% - #1: Model Reference OAE
%   Subtract output 4 from output 1 (for each freq)
%   This leaves us just with the OAE.
%   The stimulus artifact is subtracted out.
%   The columns or OAE1 correspond to the OAEs for the individual tonebursts.
OAE1 = E(:,1:4:end) - E(:,4:4:end);


% - #2: Implement the DNLR technique for calculating the OAEs
%   OAE = in + in + (-2*in) will leave us with the non-linear OAE
% or
%   OAE = output 1 + output 2 + output 3
%   The columns or OAE2 correspond to the OAEs for the individual tonebursts.
OAE2 = E(:,1:4:end) + E(:,2:4:end) + E(:,3:4:end);


% --- Take the fft of the OAEs
% FFT params
fftSize = 1024;
freqs = 0:fs/fftSize:fs/2;          % Frequency vector in Hz

% Take the ffts (of each column)
OAE1freq = mag2db(abs(fft(OAE1, fftSize, 1)));
OAE2freq = mag2db(abs(fft(OAE2, fftSize, 1)));

% Only need the positive spectrum
OAE1freq = OAE1freq(1:fftSize/2+1, :);
OAE2freq = OAE2freq(1:fftSize/2+1, :);


% --- Determine the latency per each toneburst frequency
latencies = zeros(1, length(f0));

% A threshold tuned by trial and error
threshold = 1.3e-5;

% Find where the instaneous derivative (slope) of the OAE is above the threshold
detect = 3e-3 * (abs(diff(OAE1)) > threshold);
for i = 1:length(f0)
    inds = find(detect(:,i));
    latencies(i) = inds(1) / fs;
end


% --- Fit an exponential curve to the latencies
%expFit = fit(f0, latencies, 'exp1');
expFit = load('evan.mat');
xx = linspace(0, 7e3, 500);
yy = expFit.a * exp(expFit.b * xx);



%% Plotting

% --- Plot the Model Reference OAE for each frequency
% - Impulse Response and Frequency Response pair
h = figure(1);
set(h, 'Position', [0, 0, 1400, 800]);
subplotNum = 1;
i = 1;
suptitle('Model Reference OAE');
while i <= length(f0)
    subplot(length(f0), 2, subplotNum);
    plot(1e3*t, OAE1(:,i));
    title(sprintf('%d Hz Toneburst  -  Impulse Reponse', f0(i)));
    xlabel('Time (ms)');
    ylabel('Amplitude');
    ylim([-1.65e-2 1.65e-2]);
    grid on;
    
    subplot(length(f0), 2, subplotNum+1);
    plot(freqs/1e3, OAE1freq(:,i));
    title(sprintf('%d Hz Toneburst  -  Frequency Reponse', f0(i)));
    xlabel('Frequency (kHz)');
    xlim([0 fs/2/1e3]);
    ylabel('Magnitude (dB)');
    ylim([-120 12]);
    grid on;
    
    subplotNum = subplotNum + 2;
    i = i + 1;
end
%saveas(h, 'tboae1.png');

% --- Plot each DNLR OAE
% - Impulse Response and Frequency Response pair
h = figure(2);
set(h, 'Position', [0, 0, 1400, 800]);
subplotNum = 1;
i = 1;
suptitle('Derived Non-Linear OAE');
while i <= length(f0)
    subplot(length(f0), 2, subplotNum);
    plot(1e3*t, OAE2(:,i));
    title(sprintf('%d Hz Toneburst  -  Impulse Reponse', f0(i)));
    xlabel('Time (ms)');
    ylabel('Amplitude');
    ylim([-1.65e-2 1.65e-2]);
    grid on;
    
    subplot(length(f0), 2, subplotNum+1);
    plot(freqs/1e3, OAE2freq(:,i));
    title(sprintf('%d Hz Toneburst  -  Frequency Reponse', f0(i)));
    xlabel('Frequency (kHz)');
    xlim([0 fs/2/1e3]);
    ylabel('Magnitude (dB)');
    ylim([-120 12]);
    grid on;
    
    subplotNum = subplotNum + 2;
    i = i + 1;
end
%saveas(h, 'tboae2.png');

% --- Plot and overlay the starts of the OAEs, one at a time to view the
% latency
h = figure(3);
for i = 1:length(f0)
    plot(1e3*t, OAE1(:,i));
    title('Overlaid OAEs');
    xlabel('Time (ms)');
    xlim([0 7]);
    ylabel('Amplitude');
    ylim([-1.65e-2 1.65e-2]);
    grid on;
    hold on;
    
    %saveas(h, sprintf('tboae3_%d', i), 'png');
    %waitforbuttonpress();
end

% --- Plot the latency as a function of frequency
h = figure(4);
stem(f0/1e3, 1e3*latencies);
hold on;
plot(xx/1e3, 1e3*yy);
title('Model Reference OAE Latency vs Toneburst Frequency');
xlabel('Frequency (kHz)');
xlim([0 7]);
ylabel('Latency (ms)');
ylim([0 2.75]);
%saveas(h, 'tboae4.png');