%%This function creates a spectrogram with variable window sizing based off
%%frequency information obtained from an initial low-resolution
%%spectrogram, in order to increase the precision of pitch-detection.
%%hpsonoff is a boolean that controls whether or not the spectrogram is a
%%Harmonic Product Spectrum. isNoisy is a boolean that specifies if the user 
%%wants the signal to have its DC component removed.

function [S_tot F T_tot Y] = varWinSpectrogram(wave, fs, hpsonoff, isNoisy)

%%For a simple synthetic sinusoid, finding the peak of each bin is
%%sufficient enough to find the approximate pitch during that time-window.
%%However, if there is a large DC component, the window sizing algorithm
%%will be ineffectual, so the signal is put through a simple high-pass filter with
%%a low cutoff frequency.
if(isNoisy ~= 0)
    cutoff = 100;
	d = fdesign.highpass('N,F3dB',2,cutoff/(fs/2));
    hd = design(d,'butter');
    wave = filter(hd, wave);
end

%The average lowest frequency the human ear can detect is about 20Hz. Corresponding
%period is 1/20 = 50ms. The length of the Hanning window in the initial
%spectrogram will be the equivalent of 50ms in sample length,
%and then taken to the nearest power of 2 for efficiency's sake.
maxwin=2^nextpow2(ceil(0.05*fs));

%%Initial low-resolution spectrogram. Windows are non-overlapping to make
%%later calculations substantially easier.
[S F T]=spectrogram(wave, hanning(maxwin), 0, maxwin, fs);
Y = displaydBSpectrogram(S,F,T);

bins = length(Y(1,:));
S_cell = cell(bins);
T_cell = cell(bins);
S_tot = [];
X_tot = [];

%%Pre-allocates memory for arrays
maxpeaks = zeros(1, bins);
indices = zeros(1, bins);
approxpitch = zeros(1, bins);
approxperiod = zeros(1, bins);
winsize = zeros(1, bins);


%%Calculates the approximates pitches of each time-window in the STFT
for i = 1:bins
    %%If specified, convert Y into a HPS
    if (hpsonoff ~= 0)
        Y(:,i) = hps(Y(:,i));
    end
    maxpeaks(i) = max(Y(:,i));
    indices(i) = find(Y(:,i) == maxpeaks(i),1,'first');
    approxpitch(i) = F(indices(i));
    approxperiod(i) = 1/approxpitch(i);
end

%%scaling is a constant factor found from experimentation
scaling = 20;
winsize = floor(approxperiod * fs * scaling);

for i = 1:bins
    if(winsize(i) > maxwin) %Window can not be larger than the length of the signal
        winsize(i) = maxwin;
    end
    if(winsize(i) <= 256) %Minimum window size is 256
        winsize(i) = 256;
    end
end

for i = 1:bins
   %%Sample space for each sub-STFT. The last bin may be slightly longer than
   %%the rest to make sure the whole signal is transformed.
   start = ((i-1)*maxwin+1);
   over = (maxwin*i);
   if(i==bins)
         over=length(wave);
   end
   %%nfft must remain constant at maxwin for the arrays to
   %%concatenate properly
   [S F T] = spectrogram(wave(start:over), hanning(winsize(i)), [], maxwin, fs);
   %%S_cell hold all of the spectrograms
   S_cell{i} = S;
   %For some reason the time indices were wrong (some were even negative), so I
   %manually calculated them
   T_cell{i} = linspace(start/fs,over/fs,length(T)+1);
   T_cell{i} = T_cell{i}(2:length(T_cell{i})); %%Prevents repeats
end

S_tot = S_cell{1};
T_tot = T_cell{1};

%%This for-loop is separate from the last in order to properly initialize
%%the size of arrays S_tot and T_tot, which is necessary for
%%concatenation
for i = 2:bins
  S_tot = [S_tot, S_cell{i}];
  T_tot = [T_tot, T_cell{i}]; 
end

figure;

newbins = length(S(1,:));

%%If specified, convert Y into a HPS
if (hpsonoff ~= 0)
    Y = abs(S_tot);
    for i = 1:newbins;
        Y(:,i) = hps(Y(:,i));
    end
    displaySpectrogram(Y, F, T);
else
Y = displaydBSpectrogram(S_tot, F, T_tot);
end

end