%% Input must be a spectrogram

function s2 = whiten(s,fmax)

%Not sure about this, but the paper seems to assume this
s = abs(s); 

%Values output by spectrogram are already in [0, 1] range?



%Square root of spectrogram (of the absolute value?)

%s1 = sqrt(abs(s)); %Version 1
s1 = s;             %Version for segmentation

%For each frame, find a quantity similar to the energy"

[Nf Nt] = size(s1);
energy = nan(Nt,1);

for t=1:Nt
   energy(t) = 1/fmax * sum(s1(:,t).^2);
end

threshold = quantile(energy, 0.2);

lowestix = find(energy<threshold);
%lowestenergy = energy(lowestix);


P = nan(Nf,1);
s2 = s1;
for f=1:Nf
    %P(f) = sqrt(eps + sum(s1(f,lowestix).^2));
    
    P(f) = sum(s1(f,lowestix)) / length(lowestix);  %Version for segmentation
    s2(f,:) = s1(f,:) ./P(f);
end











