% Data initial from 3-10

 S  = 22050;					 % sample rate
 b=120;
 fMin= 27.5;         
 fMax = 4186; 
 fb  = 2^(1/b);       
 winnam = 'hamming';
 Q = 1/(fb - 1);
windsizmax = 2205;	
 hopsiz = floor(0.15*S);  
% hopsiz=512;
[X,Fs,Bits]=wavread('Fall.wav');
% plot(X);
%bin=fix(b*log2(fMax/fMin));     % The  maxmum of k calculate the bINS
bin= 1 + fix( log(fMax/fMin)/log(fb) ); 
% fk = fMin * (fb .^ [(0:1:bin-1)]);
 Q = 1/(fb - 1);
fk = fMin * (fb .^ [(0:1:bin-1)]);

% fprintf('Fk is %.0f',fk(1));
% bin=length(fk);
Between =  fk * 2*pi./S;
% Nk=fix(S*Q./fk);  % windSizeOk

    Nk = fix (2*pi*Q ./Between);  
    dir= find(Nk > windsizmax);
 
    Nk(dir) = windsizmax;
    
    kerncos = zeros(bin, Nk(1));
    kernsin = zeros(bin, Nk(1));
    numzeros = Nk(1) - Nk;
    numzerosO2 = round(numzeros/2);
    



           
% Set the transform of kernl
      for k = 1:bin
          sz = Nk(k);
              winstr = [ winnam '(' num2str(sz) ')'];  % have the window amounts
     %Wn=0.54+0.46*cos((2*pi*k)/Nk);
     %X(k)=Wn*sin((2*pi*Q*k)/Nk);
    

     wind = eval(winstr); 
           wind = wind';
          
          numz = 1;
          if numzerosO2(k) ~= 0, numz =  numzerosO2(k); end;
          kerncos(k, numz: numz + Nk(k)-1) = (1/Nk(k)) * ...
               cos((Between(k))*( -sz/2 : sz/2 - 1 )).* wind;
        
          kernsin(k, numz: numz + Nk(k)-1) = (1/Nk(k)) * ...
                    sin((Between(k))*( -sz/2 : sz/2 - 1 )).* wind;
                
       % kerncos(k,1: Nk(k)) = (1/Nk(k))* cos((fk(k)*2*pi/S)*Nk(k)).* wind;
       % instead of Between(k);
       % kernsin(k,1:  Nk(k)) = (1/Nk(k)) *sin((fk(k)*2*pi/S)*Nk(k)).* wind; 

      end
  % Over

  
  % calculate the Xcq

[nrows ncols] = size(X);

[Crows Ccols]= size(kerncos);
nhops = 1 + floor((length(X)-windsizmax)/hopsiz) ;
Xcq = zeros(bin, nhops);

% Fill in the Xcq Matrix
    for jj = 0:nhops-1
	 	ninit = jj*hopsiz+1;
 		nfnl  =  jj*hopsiz +1 + windsizmax - 1;
		Y= X(ninit:nfnl);
	 	reOinfil = kerncos * Y;
	 	imOinfil = kernsin * Y;
 		    Xcq(:,jj+1) = ((reOinfil .^ 2) + (imOinfil .^ 2)) .^(.5);      % how happen:
    end
    Xcq = Xcq';
    [rowcq collcq] = size(Xcq);
     time = (hopsiz / (S/3) )  * [1 : rowcq] ;
    %time=Nk/S;

% Over 

    figure
	 % waterfall(fk, time, abs(Xcq));  % view(0, 70)
     mymode=gray;
     mymode=mymode(64:-1:1,:);
    imagesc(time,fk,abs(Xcq));
    %  specgram(abs(Xcq),512,S,512,384);
    %plot(bin,time);
   % set( gca, 'XScale', 'log'); axis tight; grid off
	%set(gca,'XTick', [100 300 1000 3000 ])