% This script is intended to speed up the work of a student who has been
% given many hours of sound files and told to listen for bird songs.
% Answer the INITIAL QUESTIONS below and press Run when instructed.
% You will still need to know what your bird sounds like.  Good luck!
% - Dr. Benjamin Shank, Hope College Physics Dept. (shank@hope.edu)

% P.S. This script has to run from the same folder as MakeSpectrogram.m,
% MakeTimeStamp.m and SuggestDynamicDB.m.  Please move these files together.
% This script requires the DSP Toolbox and the Curve Fitting Toolbox.

clear variables  % Not a question, just making space

% TO START, WHERE ARE YOUR SONG FILES?
readpath = 'D:\RCCB zoom recordings\anatomy birds';

% WHICH FILE DO YOU WANT PARSED INTO SONGS?
filename = 'M_gr997_b7_12-37_211115-T003_Tr1.WAV';

% WHERE WOULD YOU LIKE THE SONGS SAVED? (UNCOMMENT YOUR FAVORITE SUGGESTION)
%writepath = readpath;    %Save songs to same folder
writepath = [readpath '\' filename(1:end-4)]; %New local sub-folder with this file name
%writepath = 'C:\Somewhere\Else\Entirely';
    %Careful! Matlab will create whatever mis-typed file tree you enter here!

% HOW MANY SECONDS BETWEEN NOTES BEFORE THE BIRD HAS "STOPPED SINGING"?
songGap = 1.0;       %Ask your advisor about literature standards

% WHAT IS THE MINIMUM LENGTH OF A SONG (IN SECONDS)?
shortestSong = 0.1;  %Use this to cut out short flutters, etc.

% WHAT IS THE MAXIMUM LENGTH OF A SONG (IN SECONDS)?
longestSong = 60;  %Set high (This just protects against certain noise lockups)

% HOW MUCH SILENCE (IN SECONDS) WOULD YOU LIKE BEFORE AND AFTER EACH SAVED SONG?
bufferTime = 0.4;  % Must be half of songGap or less

% WHAT FREQUENCY RANGE (IN HERTZ) DOES YOUR BIRD SING IN?
Fmin = 900;    %Minimum  (Typically >500 Hz, most room/traffic noise is <1500 Hz)
Fmax = 20000;   %Maximum  (Typically about 10000 Hz, but start high.)
% DO YOU WANT TO BANDPASS FILTER THE SAVED SONGS? (TRUE/FALSE)
bFilterBandPass = true;

% HOW LONG (IN SAMPLES) WOULD YOU LIKE YOUR SPECTROGRAM BUFFERS?
frameN = 512;  %Commonly use 512 points -> 11-12 msec in modern 44-48 kHz formats

% WHAT IS THE EXPECTED DYNAMIC RANGE OF YOUR SOUND FILE (IN DECIBELS)?
dynamicDB = 35;     %Expected dB Variation Between Loudest And Quietest Note
% DO YOU WANT HELP CHOOSING dynamicDB? (TRUE/FALSE)
bSuggestDynamicDB = true;

% DO YOU WANT TO SEE THE SPECTROGRAM FOR EACH SONG? (TRUE/FALSE)
bShowSpect = true;   %Reading spectrograms is a key skill. Set this to true.

% DO YOU WANT TO HEAR EACH SONG? (TRUE/FALSE)
bSound = true;  %Requires speakers

% PLEASE CHOOSE A FEW SPECTROGRAM DISPLAY SETTINGS
logDepth = 5;    %PSD Orders of Magnitude On Spectrogram (Typically 4-6)
yScale = 'linear';       %'log' or 'linear' Frequency Axis
% SHOULD A SPECTROGRAM BE SAVED FOR EACH SONG? (TRUE/FALSE)
bSaveSpect = false;     %Spectrogram saving works even if bShowSpect = false

% *****  ADVANCED FEATURES  *****
% SHOULD PARAMETERS FROM THIS SEARCH BE SAVED WITH THE RESULTS? (TRUE/FALSE)
% (This sets you up for forthcoming syllable classification scripts)
bSaveParameters = true;

% PLEASE SELECT PEAK TRACKER SETTINGS. (10,4,0.18 works for most small birds)
noteDB = 10;         %Expected dB Variation Between 'On' and 'Off'
instDB = 4;         %Instant Drop In Tracked Peak When Note Ends
decayTime = .18;    %Time Constant Of Peak Tracker In Seconds

% WHAT FREQUENCY RANGE (IN HERTZ) CONTAINS ANNOYING IN-BAND NOISE?
% (If bFilterDetect = false and bFilterSave = false, this range is unused.)
FminNoise = 2400;    % This was added to deal with a shrill fan
FmaxNoise = 2600; 
% DO YOU WANT TO USE A NOTCH FILTER WHILE DETECTING SONGS? (TRUE/FALSE)
bFilterDetect = false;
% DO YOU WANT TO USE A NOTCH FILTER ON SAVED SONGS? (TRUE/FALSE)
bFilterSave = false;      %Set to false unless the fan is deafening

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OKAY, MATLAB CAN TAKE IT FROM HERE.  PRESS THE RUN BUTTON OR HIT F5.
% WE'LL MEET YOU IN THE COMMAND WINDOW AT THE BOTTOM OF THIS SCREEN.

% By the way, for every (y/n) question you can also:
% Type s to hear the latest sound again
% Type c to exit the program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%Get Basic File Info
info = audioinfo([readpath '\' filename]);
fileLength = info.TotalSamples;
Fs = info.SampleRate;

if Fmax >= Fs/2
    Fmax = Fs/2-1;
    disp(['Max frequency cannot exceed half the sample rate. We changed it to ' num2str(Fmax,'%i') ' Hz.'])
end
if FmaxNoise >= Fs/2
    FmaxNoise = Fs/2-1;
    disp(['Max frequency for noise removal cannot exceed half the sample rate. We changed it to ' num2str(Fmax,'%i') ' Hz.'])
end

if ~isfolder(writepath)
    mkdir(writepath);
end
if bSaveParameters
    save([writepath '\SongFinderParameters.mat']); %Save search variables for future analysis
end

%Set Up Filter(s) if Desired
if bFilterBandPass
    [bPass,aPass] = butter(4,[Fmin Fmax]/(Fs/2), 'bandpass');
end
if bFilterDetect || bFilterSave
    [bNotch,aNotch] = butter(4,[FminNoise FmaxNoise]/(Fs/2), 'stop');
end

frameTime = frameN/Fs;
window = ones(frameN,1);   %Blackman Windowing Function
for n = 1:frameN
    window(n) = (0.42 - 0.5*cos(2*pi*n/(frameN-1)) + 0.08*cos(4*pi*n/(frameN-1)));
end
windRMS = sqrt(mean(window.^2));  %Make sure the window does not change the RMS level
window = window/windRMS;

if bShowSpect || bSaveSpect
    close all   %Get rid of any open graphs (We could be making a lot of new ones)
    Acutoff = 10^(-logDepth);   %This Is Only Used To Scale The Spectrogram
    F = (1:frameN/2)/frameN*Fs;   %Master List Of Frequencies To Analyze
    cmap = hot(64);    %Custom Colormap Converts to B&W Well
    cmap = cmap(end:-1:1,:);
end

%%
% We're just going to go through the file once to get RMS volume data
t = frameTime/2 * (1:floor(fileLength*2/frameN));     %Place To Store Time Variable
RMS = zeros(size(t));  %Place To Store Single-Frame RMS Values
frameEnd = frameN;  %End Of First Frame
i = 1; %Storage Index
jumpAhead = 0;  %File Index Offset for latest read (Look, these files are big!)
y = [];   %Dummy read to get things started
bStopEarly = false;   %This catches an error in audioread that can abort the script
while frameEnd < fileLength - frameN/2 && ~bStopEarly
    if frameEnd > length(y)+jumpAhead  %We don't have enough data for the next frame
        jumpAhead = frameEnd - frameN; %Start at the next unprocessed window
        readSize = min(1e6,fileLength-jumpAhead); %Read 1 million points or to end of file
        y = audioread([readpath '\' filename],[jumpAhead+1, jumpAhead+readSize]);
        y = y(:,1);        %Grab Only The First Channel
        if bFilterBandPass
            y = filtfilt(bPass,aPass,y);  %filtfilt Prevents Phase Delays, Which Limits Peak Distortion
        end
        if bFilterDetect
            y = filtfilt(bNotch,aNotch,y);
        end
    end  %Now we should have plenty of data. Go back to making "spectrogram" buffers.

    frame = y(frameEnd-frameN+1-jumpAhead : frameEnd-jumpAhead);  %Grab The Next Data Set...
    frame = frame.*window;                    %...And Window It
    RMS(i) = std(frame);         %Store The Total Volume Level In This Frame
    i = i+1;
    frameEnd = frameEnd + frameN/2; %Standard 50% Overlap
    if length(y) < readSize && frameEnd > length(y)+jumpAhead
        bStopEarly = true;
    end
end   %Now we have volume vs time data


%%
% Find the loudest sound and make sure it is a bird
displayIndex = 1;  %Figure number for display
decayFrames = floor(2*decayTime/frameTime); %# of windows in one decay time
response = 'a';  %#ok<*NASGU> %Non-sense initialization 
bContinue = bShowSpect || bSound;
while (bShowSpect || bSound) && bContinue %Keep asking until loudest sound is a bird
    z1 = 10^(noteDB/20);    %RMS Ratio Between 'On' and 'Off'
    maxIdx = find(RMS == max(RMS),1);
    startIdx = find(RMS(1:maxIdx-decayFrames) < max(RMS)/z1,1,"last");
    stopIdx = maxIdx+decayFrames + find(RMS(maxIdx+decayFrames+1:end)<max(RMS)/z1,1,"first");
    readStart = max(startIdx*frameN/2, 1);
    readStop = min(stopIdx*frameN/2, fileLength);
    y = audioread([readpath '\' filename],[readStart, readStop]);
    y = y(:,1);        %Grab Only The First Channel
    if bFilterBandPass
        y = filtfilt(bPass,aPass,y); 
    end
    if bFilterDetect
        y = filtfilt(bNotch,aNotch,y); 
    end
    if bShowSpect
        MakeTimestamp   %Secondary script in this folder
        MakeSpectrogram   %Secondary script in this folder
    end
    if bSound
        sound(y,Fs);
    end
    response = input("Here's the loudest sound. Is this a bird? (y/n) ","s");
    if isempty(response)  %So an accidental Enter doesn't give an error code
        response = 's';
    end
    if bShowSpect
        pause(0.1)
        figure(displayIndex);  %Just bringing this to the front again
    end
    if response == 'y' || response == 'Y' || response == "yes"
        bContinue = false;
    elseif response == 's' || response == 'S' || response == "sound"
        while response == 's' || response == 'S' || response == "sound"
            sound(y,Fs);
            response = input("Here's the sound again. Is this a bird? (y/n) ","s");
            if isempty(response)  %So an accidental Enter doesn't give an error code
                response = 's';
            end
            if response == 'y' || response == 'Y' || response == "yes"
                bContinue = false;
            elseif response == 'c' || response == 'C' || response == "close"
                disp("Song finder terminated by user")
                return
            end
        end
    elseif response == 'c' || response == 'C' || response == "close"
        disp("Song finder terminated by user")
        return
    else
        RMS(startIdx:stopIdx) = 0; %Ignore the offending loud noise from now on
    end
end  %End of 'loudest sound' check

% Provide dynamic range help if asked
if bSuggestDynamicDB
    SuggestDynamicDB
end

%%
% Hopefully we've now removed all the loud extraneous noises from this
% file. Next, go through and identify all the "notes".

z1 = 10^(noteDB/20);    %RMS Ratio Between 'On' and 'Off'
z2 = 10^(instDB/20);    %Instantaneous 'On' Threshold Shift At Each 'Off' (Optional)
z3 = 10^(dynamicDB/20); %Ratio Between Noisest And Quietest Notes
maxRMS = max(RMS);
minLvl = maxRMS / z3; %How Low Can The 'On' Threshold Decay To?
mult = 1 - frameTime/(2*decayTime);  %Convert Time Constant To Multiplier In Each Frame
startIdx = [];        %(Empty) List of Note Parameters  
stopIdx = [];
pk = minLvl;
bNote = false;  %Are We In The Middle Of A Note?
i = 1;          %Block Index
while i <= length(t)
    if ~bNote && RMS(i) > pk
        bNote = true;
        startIdx = [startIdx i]; %#ok<*AGROW>
    end
    if bNote && RMS(i) < pk/z1
        bNote = false;
        stopIdx = [stopIdx i];
        pk = pk/z2; %Implement quikDB
    end
    decayedPk = min([mult*pk RMS(i)*z3]);  %Drop The Peak Rapidly If The Sound Really Cuts Out
    pk = max([decayedPk RMS(i) minLvl]);   %Can't Drop Below Current Volume Or Dynamic Threshold
    i = i+1;  %Move To Next Frame
end

%%
% Now use the songGap to group "notes" into "songs"
gapFrames = ceil(2*songGap/frameTime);
i = 1;
while i < length(startIdx)  %String together notes with short gaps in between
    if startIdx(i+1) < stopIdx(i)+gapFrames
        startIdx(i+1) = [];
        stopIdx(i) = [];
    else
        i = i+1;
    end
end
minFrames = ceil(2*shortestSong/frameTime);
maxFrames = floor(2*longestSong/frameTime);
i = 1;
while i < length(startIdx)  %Cut out songs that are too short
    if stopIdx(i) - startIdx(i) < minFrames
        startIdx(i) = []; %#ok<*SAGROW> 
        stopIdx(i) = [];
    end
    if stopIdx(i) - startIdx(i) > maxFrames
        stopIdx = [stopIdx(1:i-1) startIdx(i)+maxFrames-gapFrames stopIdx(i:end)];
        startIdx = [startIdx(1:i) startIdx(i)+maxFrames startIdx(i+1:end)];
    end
    i = i+1;
end

%%
% Best to show them what's going on before we write a bunch of memory.
displayIndex = 1;  %Figure number for display
if bufferTime > songGap/2
    bufferTime = songGap/2;
end
bufferN = floor(bufferTime*Fs);  %How many points before and after each song?
filestub = filename(1:end-4); %Remove the extension
response = 'y';  % Start with a positive question
while displayIndex <= length(stopIdx)
    readStart = max(startIdx(displayIndex)*frameN/2 - bufferN, 1);
    readStop = min(stopIdx(displayIndex)*frameN/2 + bufferN, fileLength);
    y = audioread([readpath '\' filename],[readStart, readStop]);
    y = y(:,1);        %Grab Only The First Channel
    if bFilterBandPass
        y = filtfilt(bPass,aPass,y);
    end
    if bFilterSave
        y = filtfilt(bNotch,aNotch,y);
    end
    MakeTimestamp   %Secondary script in this folder (Used in saved filename)
    if bShowSpect || bSaveSpect
        MakeSpectrogram   %Secondary script in this folder
    end
    if bSound
        sound(y,Fs);
    end
    if bShowSpect || bSound
        if bShowSpect
            pause(0.1)
            figure(displayIndex);  %Just bringing this to the front again
        end
        if response == 'y' || response == 'Y' || response == "yes"
            response = input(['Here is song #' num2str(displayIndex,'%i') '. Is this a bird? (y/n) '],"s");

        else
            response = input(['Okay. How about this for song #' num2str(displayIndex,'%i') '? (y/n) '],"s");
        end
        if isempty(response)  %So an accidental Enter doesn't give an error code
            response = 's';
        end
        if response == 'y' || response == 'Y' || response == "yes"
            audiowrite([writepath '\' filestub '_Song' num2str(displayIndex,'%.3i')...
                '_' timestamp '.wav'],y,Fs);
            if bSaveSpect
                saveas(fig,[writepath '\' filestub '_Song' num2str(displayIndex,'%.3i')...
                    '_' timestamp '.png'])
            end
            displayIndex = displayIndex + 1;  % Go to next "song"
        elseif response == 's' || response == 'S' || response == "sound"
            while response == 's' || response == 'S' || response == "sound"
                sound(y,Fs);
                response = input(['Here is song #' num2str(displayIndex,'%i') ' again. Is this a bird? (y/n) '],"s");
                if isempty(response)  %So an accidental Enter doesn't give an error code
                    response = 's';
                end
                if response == 'y' || response == 'Y' || response == "yes"
                    audiowrite([writepath '\' filestub '_Song' num2str(displayIndex,'%.3i')...
                        '_' timestamp '.wav'],y,Fs);
                    if bSaveSpect
                        saveas(fig,[writepath '\' filestub '_Song' num2str(displayIndex,'%.3i')...
                            '_' timestamp '.png'])
                    end
                    displayIndex = displayIndex + 1;  % Go to next "song"
                elseif response == 'n' || response == 'N' || response == "no"
                    %Delete this "song" to preserve the numbering scheme and move on
                    startIdx(displayIndex) = [];
                    stopIdx(displayIndex) = [];
                elseif response == 'c' || response == 'C' || response == "close"
                    if displayIndex > 1
                        startIdx = startIdx(1:displayIndex-1);
                        stopIdx = stopIdx(1:displayIndex-1);
                        if bSaveParameters
                            save([writepath '\SongFinderSummary.mat'],'maxRMS','startIdx','stopIdx'); %Save result variables for future analysis
                        end
                    end
                    disp("Song finder terminated by user (Songs to this point were saved)")
                    return
                end
            end
        elseif response == 'c' || response == 'C' || response == "close"
            if displayIndex > 1
                startIdx = startIdx(1:displayIndex-1);
                stopIdx = stopIdx(1:displayIndex-1);
                if bSaveParameters
                    save([writepath '\SongFinderSummary.mat'],'maxRMS','startIdx','stopIdx'); %Save result variables for future analysis
                end
            end
            disp("Song finder terminated by user (Songs to this point were saved)")
            return
        else
            %Delete this "song" to preserve the numbering scheme and move on
            startIdx(displayIndex) = [];
            stopIdx(displayIndex) = [];
        end
    else  % They don't want to check *anything*.  Write all the files.
        audiowrite([writepath '\' filestub '_Song' num2str(displayIndex,'%.3i')...
            '_' timestamp '.wav'],y,Fs);
        if bSaveSpect
            saveas(fig,[writepath '\' filestub '_Song' num2str(displayIndex,'%.3i')...
                '_' timestamp '.png'])
        end
        displayIndex = displayIndex + 1;  % Go to next "song"
    end
end

if bSaveParameters
    save([writepath '\SongFinderSummary.mat'],'maxRMS','startIdx','stopIdx'); %Save result variables for future analysis
end
disp(['Great, that should be it! We identified ' num2str(length(stopIdx),'%i') ' songs.'])
disp(['Songs were saved to ' writepath])