% Evan Shenkman
% MMI610
% Project 1
% Distortion Product OAEs

% --- Setup simulation params
fs = 48000;                                 % sampling rate
t = 0:1.0/fs:0.04;                          % time vec
f1 = [500 1000 1500 2000 3000 4000 6000];   % f1s
f2 = 1.22 * f1;                             % f2s
minCycles = 5;                              % number of periods in each toneburst
minSamples = 1/min(f1) * fs * minCycles;    % maintain signal length for equal energy
f_dp = 2*f1 - f2;                           % lower cubic distortion product frequencies
nsigs = 3*length(f1);

% --- Model Params
spl = 70 * ones(1,nsigs);
spl(3:3:end) = 70 + mag2db(2);
irron = ones(1,nsigs);
normalizeRms = ones(1,nsigs);
subjectNo = 1;
fc = 'all';

% --- Signal Generation
% The input to the model is Nx(3*length(f1)) (i.e. 21 N-point signals)
% There are three signals per frequency
% - Sig 1: a positive, 5 cycle, scaled summation of f1 and f2
% - Sig 2: a positive, 5 cycle, scaled summation of f1 and f2
% - Sig 3: a negative, 5 cycle, scaled summation of f1 and f2
insig = zeros(length(t), 3*length(f1));
i = 1;
for f_i = 1:length(f1)
    % Determine the number of cycles for each freq
    cycles = ceil(minSamples * f1(f_i) / fs);
    
    % Generate the f1 tone
    s1 = sin(2*pi*t(1:ceil(cycles*fs/f1(f_i))+1)'*f1(f_i));
    
    % Generate the f2 tone and make it the same length as s1
    s2 = sin(2*pi*t(1:ceil(cycles*fs/f2(f_i)))'*f2(f_i));
    s2 = [s2; zeros(length(s1) - length(s2), 1)]; %#ok<AGROW>
    
    % Compute the weighted sum of the signals
    s3 = s1 + s2 * db2mag(-10);
    
    % Normalize the signal
    s3 = s3 / max(abs(s3));
    
    % Assign the signals to the columns of insig
    insig(1:length(s3), i) = s3;
    insig(1:length(s3), i+1) = s3;
    insig(1:length(s3), i+2) = -1 * s3;
    
    % Increment index counter
    i = i + 3;
end


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


% --- Calculate the DPOAE (two different ways)
% - #1: DNLR technique
OAEdnlr = E(:,1:3:end) + E(:,2:3:end) + E(:,3:3:end);

% - #2: High-order BPF
OAE_preFilt = E(:,1:3:end);
OAE = zeros(size(OAE_preFilt));
filtOrder = 1024;
filters = zeros(filtOrder+1, length(f1));
for i = 1:length(f_dp)
    % Design the BPF
    hiF = (2*f_dp(i)/3 + f1(i)/3);
    loF = f_dp(i) - (hiF - f_dp(i));
    B = fir1(filtOrder, (2/fs)*[loF hiF], 'bandpass');
    
    % Save the filter IR
    filters(:,i) = B';
    
    % Apply the filter
    OAE(:,i) = filter(B, 1, OAE_preFilt(:,i));
end


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

% Take the ffts (of each column)
OAEdnlrfreq = mag2db(abs(fft(OAEdnlr, fftSize, 1)));
OAEfreq = mag2db(abs(fft(OAE, fftSize, 1)));
OAEfreq_preFilt = mag2db(abs(fft(E(:,1:3:end), fftSize, 1)));

% Only need the positive spectrum
OAEdnlrfreq = OAEdnlrfreq(1:fftSize/2+1, :);
OAEfreq = OAEfreq(1:fftSize/2+1, :);
OAEfreq_preFilt = OAEfreq_preFilt(1:fftSize/2+1, :);


% --- Compute the energy of a 20Hz band around the cubic distortion frequency
OAEenergy = zeros(1, length(f1));
for i = 1:size(OAEfreq,2)
    % Make the spectrum linear again and square it
    linFFT = (db2mag(OAEfreq(:,i))).^2;
    
    % Get the index of the DP frequency (freq resolution of 2Hz / bin)
    index = floor(f_dp(i)/10);
    
    % Integrate over the band centered around the DP frequency
    OAEenergy(i) = sum(linFFT(index-1:index+1))/3;
end

% Smooth the data
xx = linspace(0, 7e3, 32);
yy = spline(f_dp, OAEenergy, xx);


%% --- Plotting
% Plot the response of each Filter
for i = 1:size(filters,2)
    h = figure(1);
    set(h, 'Position', [100 100 900 700]);
    
    % Plot the filter response
    freqz(filters(:,i));
    subplot(2,1,1);
    ylim([-120 12]);
    subplot(2,1,2)
    ylim([-5000 1000]);
    
    suptitle(sprintf('%dHz BPF Response', f_dp(i)));
    saveas(h, sprintf('dpoae1_%d.png', i));
    
    close(h);
end

% Plot the OAEs (before and after filtering)
h = figure(2);
set(h, 'Position', [100 100 900 700]);
suptitle('DPOAE Impulse Responses');
for i = 1:7
    if i == 7
        subplot(4,2,7.5);
    else
        subplot(4,2,i);
    end
    plot(1e3*t, OAE_preFilt(:,i));
    hold on;
    plot(1e3*t, OAE(:,i));
    %plot(1e3*t, OAEdnlr(:,i));
    title(sprintf('DPOAE - %.0fHz', f_dp(i)));
    xlabel('Time (ms)');
    ylabel('Amplitude');
    legend('Pre Filt', 'BPF');
end
saveas(h, 'dpoae2.png');


% Plot the FFT of the above OAEs
h = figure(3);
set(h, 'Position', [100 100 900 700]);
suptitle('DPOAE Frequency Response');
for i = 1:7
    if i == 7
        subplot(4,2,7.5);
    else
        subplot(4,2,i);
    end
    
    plot(freqs/1e3, OAEfreq_preFilt(:,i));
    hold on;
    plot(freqs/1e3, OAEfreq(:,i));
    %plot(freqs/1e3, OAEdnlrfreq(:,i));
    hold on;
    hstem = stem([f_dp(i) f1(i) f2(i)]/1e3, -150*[1 1 1], 'ShowBaseLine', 'off');
    xpos = get(hstem, 'XData');
    ypos = get(hstem, 'YData');
    labels = {'f_{dp}' 'f1' 'f2'};
    for j = 1:1
        text(xpos(j), ypos(j) - 25, labels{j}, 'HorizontalAlignment', 'center');
    end
    title(sprintf('DPOAE using DNLR - %.0fHz', f_dp(i)));
    xlabel('Frequency (kHz)');
    xlim([0 fs/2/1e3]);
    ylabel('Magnitude');
    ylim([-200 100]);
    legend('Pre Filt', 'BPF');
end
saveas(h, 'dpoae3.png');
    

% Plot the "energy" about the DP frequency
h = figure(4);
semilogx(f_dp, 10*log10(OAEenergy), 'o')
set(gca,'xtick',f_dp);
set(gca, 'XTickLabel', strread(num2str(f_dp),'%s')); %#ok<DSTRRD>
hold on;
plot(xx, 10*log10(yy))
title('DPOAE Energy vs F1 Frequency');
xlabel('Frequency (Hz)');
xlim([0 7e3]);
ylabel('Energy (dB)');
ylim([0 35]);
saveas(h, 'dpoae4.png');