% find freqs in big wave files
reload=1;
if reload==1
    %Y=wavread('../wav/betweenMLandInfinite_1-2.wav');
    Y=wavread('../wav/outsideSmallAC_1-2.wav');
    %Y=wavread('../wav/nitrogenTanks_1-2.wav');
    %Y=wavread('../wav/underDeskML1_1-2.wav');
    %Y=wavread('../wav/doorWhack_1-2.wav');
    %Y=wavread('../wav/avtroom.wav');
    %Y=wavread('../wav/referenceAnote');
    %Y=wavread('../wav/treeLabQuad');
    %Y=wavread('../wav/elisazTelephonePole.wav');
    %Y=wavread('../wav/powerPlant.wav');
    
    yf=fft(Y(:,1));
end

yfm=abs(yf);

N = size(yfm,1);

yfm2=yfm(1:N/2)';

plot(yfm2);

topFreq=44100;

binInc = topFreq/N;

fKHz = (1:N/2)*binInc/1000;  % 1000 cuz of KHz

% apply outer and inner ear transfer function
% (from Terhardt), Calculating vertual pitch.  Hearing Research 1:155-18x,
% via traistan jehan
% AdB(fKHz) =  -3.64 * f^(-0.8) + 6.5 exp (-0.6 * f-3.3)^2) - 10^-3 * f^4
AdB =  -3.64 * fKHz.^(-0.8) + 6.5 .* exp (-0.6 * (fKHz -3.3).^2) - 10^-3 .* fKHz.^4;
scale = 10.^(AdB./20);

plot(fKHz(800:N/2),AdB(800:N/2))

% scaled magnitude 
sMag = yfm2 .* scale;
%sMag = yfm2;

plot(fKHz(800:N/2),sMag(800:N/2))


%%[v,i] = localmax(yfm);
%[v,i] = sort(yfm2);
%v2 = flipdim(v,1);
%i2 = flipdim(i,1);

%hist(i2(1:10000)*binInc,1000)

% bin them into semitones (no harmonics)
energy1=zeros(1,120);
lowC=8.125;
for semi=0:120
    fun1=lowC*2^((semi-.5)/12);
    fun2=lowC*2^((semi+.5)/12);
    bin1 = round(fun1/binInc);
    bin2 = round(fun2/binInc);
    %nrg(1+semi) = sum(sMag(bin1:bin2));    
    nrg(1+semi) = sum(sMag(bin1:bin2))/(bin2-bin1+1);    
    stone(1+semi) = semi;
end

plot(stone,nrg,'.');

[v,i] = sort(nrg);
v2 = flipdim(v,2);
i2 = flipdim(i,2);
i3= i2 -1;    % stone(1) contains data for semitone 1
i3(1:12)


% bin them into semitones (including harmonics)
energy1=zeros(1,120);
lowC=8.125;
for semi=0:120
    for harm=1:20
        fun1=lowC*2^((semi-.5)/12)*harm;
        fun2=lowC*2^((semi+.5)/12)*harm;
        bin1 = round(fun1/binInc);
        bin2 = round(fun2/binInc);
        if bin1 < N/2
            if bin2 < N/2
                %nrg(1+semi) = sum(sMag(bin1:bin2))/(bin2-bin1+1);    
                nrg(1+semi) = sum(sMag(bin1:bin2)) * harm^-1;    
            end
        end
    end
    stone(1+semi) = semi;
    semi;
end

plot(stone,nrg,'.'); 
mnrg=mean(nrg);

[v,i] = sort(nrg);
v2 = flipdim(v,2);
i2 = flipdim(i,2);
i3= i2 -1;    % stone(1) contains data for semitone 1
i3(1:12)