% Evan Shenkman
% MMI610
% Project 1
% Click Evoked OAEs

% --- Setup the simulation params
fs = 48000;                                     % sampling rate
t = 0:1.0/fs:0.05;                              % time vec


% -- Signal Generation
% The input to the model is Nx4 where N is the number of samples
% - Sig 1: a positive 100us pulse
% - Sig 2: a positive 100us pulse
% - Sig 3: a negative 100us pulse
% - Sig 4: a positive 100us pulse

% Start with all zeros
insig = zeros(length(t),4);

% Determine the pulse width
pulseWidth = ceil(1e-4 * fs)+1;

% Add the pulse to the zeros
insig(1:pulseWidth, :) = ones(pulseWidth, 4);

% Invert the "phase" of the third signal
% - used for the DNLR technique
insig(:,3) = -1 * insig(:,3);


% --- Model params
spl = [50,50,50+mag2db(2),50];  % make signal 3, the negative pulse, twice as large for DNLR
irron = [1,1,1,0];              % use irregularities for all but the last signal (reference)
normalizeRms = [1,1,1,1];
subjectNo = 1;
fc = 'all';


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


% --- Calculate the OAE (three different ways)
% - #1: Model Reference OAE
%   Subtract output 4 from output 1
%   This leaves us just with the OAE.
%   The stimulus artifact is subtracted out.
OAE1 = E(:,1) - E(:,4);


% - #2: Window the single output -- 5ms window
%   Manually remove the signal artifact using a rectangular window
% Make the window
window = ones(length(t), 1);
windowLen = 5e-3;
windowWidth = ceil(windowLen * fs);
% Apply the window
window(1:windowWidth) = zeros(windowWidth, 1);
OAE2 = window .* E(:,1);


% - #3: Implement the DNLR technique for calculating the OAE
%   OAE = in + in + (-2*in) will leave us with the non-linear OAE
%   or
%   OAE = output 1 + output 2 + output 3
OAE3 = E(:,1) + E(:,2) + E(:,3);


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

% Take the three ffts
OAE1freq = mag2db(abs(fft(OAE1, fftSize)));
OAE2freq = mag2db(abs(fft(OAE2, fftSize)));
OAE3freq = mag2db(abs(fft(OAE3, fftSize)));

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



%% Plotting

% --- Plot the three OAEs
% - #1: The "reference" OAE
h = figure(1);
set(h, 'Position', [100, 100, 800, 500]);
subplot(3,1,1)
plot(1e3*t, OAE1);
title('Model Reference OAE  -  Impulse Reponse');
xlabel('Time (ms)');
ylabel('Amplitude');
ylim([-4e-3 4e-3]);
subplot(3,1,2)
plot(freqs/1e3, OAE1freq);
title('Model Reference OAE  -  Frequency Reponse');
xlabel('Frequency (kHz)');
xlim([0 fs/2/1e3]);
ylabel('Magnitude (dB)');
ylim([-120 6]);
subplot(3,1,3)
sgram(OAE1, fs, 'nocolorbar');
title('Model Reference OAE  -  Spectrogram');
%saveas(h, 'ceoae1.png');


% - #2: The windowed OAE
h = figure(2);
set(h, 'Position', [100, 100, 800, 500]);
subplot(3,1,1)
plot(1e3*t, OAE2);
title('Windowed OAE  -  Impulse Reponse');
xlabel('Time (ms)');
ylabel('Amplitude');
ylim([-4e-3 4e-3]);
subplot(3,1,2)
plot(freqs/1e3, OAE2freq);
title('Windowed OAE  -  Frequency Reponse');
xlabel('Frequency (kHz)');
xlim([0 fs/2/1e3]);
ylabel('Magnitude');
ylim([-72 6]);
subplot(3,1,3)
sgram(OAE2, fs, 'nocolorbar');
title('Windowed OAE  -  Spectrogram');
%saveas(h, 'ceoae2.png');


% - #3: The DNLR OAE
h = figure(3);
set(h, 'Position', [100, 100, 800, 500]);
subplot(3,1,1)
plot(1e3*t, OAE3);
title('Derived Non-Linear OAE  -  Impulse Reponse');
xlabel('Time (ms)');
ylabel('Magnitude');
ylim([-4e-3 4e-3]);
subplot(3,1,2)
plot(freqs/1e3, OAE3freq);
title('Derived Non-Linear OAE  -  Frequency Reponse');
xlabel('Frequency (Hz)');
xlim([0 fs/2/1e3]);
ylabel('Magnitude (dB)');
ylim([-120 6]);
subplot(3,1,3)
sgram(OAE3, fs, 'nocolorbar');
title('Derived Non-Linear OAE  -  Spectrogram');
%saveas(h, 'ceoae3.png');