%% Input must be a spectrogram

function s2 = whiten(s,fmax)

%Values output by spectrogram are already in [0, 1] range

%Square root of spectrogram (of the absolute value?)
s1 = sqrt(abs(s)); 
%s1 = sqrt(s);

%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));
    s2(f,:) = s1(f,:) ./P(f);
end