%noisetest
%generate white noise
whiten = rand(60000,1)-0.5;  % from -0.5 to 0.5

cuts=[0.25 2 5 10 18 25]; % filterbank cutoff frequencies
fs = 50.0;
framehz=0.5;  % frequency in Hz we are taking frames in filterbank
[enoisewh,time]=auditoryfbank(whiten,fs,framehz,cuts);  % filterbanks
colorwh13 = enoisewh(:,1) - enoisewh(:,3);
dbwhite= median(colorwh13)  % answer is -4.9db! for white noise!


% generate pink or 1/f noise
% http://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html
Nx = 2^16;  % number of samples to synthesize
B = [0.049922035 -0.095993537 0.050612699 -0.004408786];
A = [1 -2.494956002   2.017265875  -0.522189400];
nT60 = round(log(1000)/(1-max(abs(roots(A))))); % T60 est.
v = randn(1,Nx+nT60); % Gaussian white noise: N(0,1)
pinkn = filter(B,A,v);        % Apply 1/F roll-off to PSD
pinkn = pinkn(nT60+1:end);    % Skip transient response

[enoisepink,time]=auditoryfbank(pinkn,fs,framehz,cuts);  % filterbanks
colorpink13 = enoisepink(:,1) - enoisepink(:,3);
dbpink=median(color13)  % answer is +4.47db! for pink noise!
diff = dbwhite-dbpink
% therefor our events have power propto f^alpha where 
% alpha = (mycolor - dbwhite)/diff
% our 10 db corresponds to alpha about -1.6
clim = 10.0; % db for our clim
alpha = (clim - dbwhite)/diff
