Dawn@4: function [ASF] = AudioSpectralFlatness(raw, fs, noOfFrames, frameLength)
Dawn@4: 
Dawn@4: check this works!
Dawn@4: 
Dawn@4: hopSize = frameLength;
Dawn@4: windowsize = frameLength;
Dawn@4: FFTsize = frameLength;
Dawn@4: windowType = hamming(frameLength);
Dawn@4: 
Dawn@4: % bin = fs/512;
Dawn@4: % bandsNumber = 24;
Dawn@4: % loEdge = 250;
Dawn@4: % hiEdege = 16000;
Dawn@4: 
Dawn@4: 
Dawn@4: loEdge = 250;
Dawn@4: hiEdge = 16000;  % Setting default hiedge
Dawn@4: 
Dawn@4: if (hiEdge >= fs/2) % Check if value for hiedge is valid
Dawn@4:     % sr/2-1 : Skipping extraction up to sr/2; Not possible due to 5 percent band overlap
Dawn@4:   hiEdge =  min(hiEdge,fs/2-1); 
Dawn@4: end
Dawn@4: hiEdge = 2^(floor(log2(hiEdge/1000)*4)/4)*1000 ; % Setting exact value for hiEdge, rounding to next lower frequency if necessary
Dawn@4: if (hiEdge*1.05 >= fs/2) 
Dawn@4:     hiEdge = hiEdge/2^0.25 ; 
Dawn@4: end; %Now it's possible to check if hiEdge is valid
Dawn@4: 	
Dawn@4: N = frameLength;
Dawn@4: 
Dawn@4: numbands = floor(4*log2(hiEdge/loEdge));
Dawn@4: firstband = round(log2(loEdge/1000)*4);
Dawn@4: overlap = 0.05;
Dawn@4: grpsize = 1;
Dawn@4: 
Dawn@4: [fftout,phase] = mpeg7getspec( raw, fs, hopSize, windowsize, windowType, FFTsize );
Dawn@4: band=[];
Dawn@4: for k = 1:numbands
Dawn@4:   f_lo = loEdge * (2^((k-1)/4)) * (1-overlap);
Dawn@4:   f_hi = loEdge * (2^((k  )/4)) * (1+overlap);
Dawn@4:   
Dawn@4:   i_lo = round( f_lo/(fs/N) ) + 1;
Dawn@4:   i_hi = round( f_hi/(fs/N) ) + 1;
Dawn@4:   
Dawn@4:   fband(k,1) = f_lo;
Dawn@4:   fband(k,2) = f_hi;
Dawn@4:   
Dawn@4:   iband(k,1) = i_lo;
Dawn@4:   iband(k,2) = i_hi;
Dawn@4: 
Dawn@4:   % Rounding of upper index according due to coefficient grouping
Dawn@4:   if (k+firstband-1 >= 0)                   %Start grouping at 1kHz 
Dawn@4:     grpsize = 2^ceil( (k+firstband )/4);
Dawn@4:     i_hi = round((i_hi-i_lo+1)/grpsize)*grpsize + i_lo-1 ;
Dawn@4:   else
Dawn@4:     grpsize = 1;
Dawn@4:   end
Dawn@4:   tmp = fftout(i_lo:i_hi,:) .^ 2;         % PSD coefficients
Dawn@4:   ncoeffs = i_hi - i_lo + 1;
Dawn@4:  
Dawn@4:   if (k+firstband-1 >= 0)                   % Coefficient grouping
Dawn@4:     tmp2 = tmp(1:grpsize:ncoeffs,:);
Dawn@4:     for g=2:grpsize
Dawn@4:       tmp2 = tmp2 + tmp(g:grpsize:ncoeffs,:) ;
Dawn@4:     end
Dawn@4:     tmp = tmp2;
Dawn@4:   end
Dawn@4:   % Actual calculation
Dawn@4:   ncoeffs = ncoeffs/grpsize ;
Dawn@4:   tmp = tmp + 1e-50;       % avoid underflow for zero signals
Dawn@4:   gm(k,:) = exp( sum(log(tmp))/ncoeffs ); % log processing avoids overflow
Dawn@4:   am(k,:) = sum(tmp) / ncoeffs;
Dawn@4: end
Dawn@4: ASF = (gm./am)';
Dawn@4: 
Dawn@4: % myColor = 'rx';
Dawn@4: % for i=1:(length(raw)/frameLength)
Dawn@4: %     subplot(411); plot( raw( (i-1)*frameLength+1 : (i*frameLength) ));
Dawn@4: %     subplot(412); plot( fftout(:,i) );
Dawn@4: %     subplot(425); plot( i, ASF(i,1), myColor ); hold on;
Dawn@4: %     subplot(426); plot( i, ASF(i,2), myColor ); hold on;
Dawn@4: %     subplot(427); plot( i, ASF(i,3), myColor ); hold on;
Dawn@4: %     subplot(428); plot( i, ASF(i,4), myColor ); hold on;
Dawn@4: % end
Dawn@4: % lo_edge = loEdge;
Dawn@4: % hi_edge = hiEdge;
Dawn@4: 
Dawn@4: % for n = 1:noOfFrames
Dawn@4: %     j = (n-1) * 512 + 1;
Dawn@4: %     k = n * 512;
Dawn@4: %     if(k<=length(raw))
Dawn@4: %         temp = raw(j:k);
Dawn@4: %     end
Dawn@4: %     if(k>length(raw))
Dawn@4: %         
Dawn@4: %         temp = [raw(j:end)' zeros(512-length(raw(j:end)),1)'];
Dawn@4: %     end
Dawn@4: %     
Dawn@4: %         FFTX=fft(temp,hopSize);
Dawn@4: %         subplot(311); plot(abs(FFTX(1:length(FFTX)/2)));
Dawn@4: %         BW=1.5*fs*1024;
Dawn@4: %         PSD=FFTX/BW;
Dawn@4: %      
Dawn@4: %         for b = 1:24                                
Dawn@4: %             floKb = 0.95*loEdge * 2^(0.25*b-0.25);
Dawn@4: %             fhiKb = 1.05*loEdge * 2^(0.25*b);
Dawn@4: %             lowKbp = round( floKb/(fs/512) ) + 1;
Dawn@4: %             hiKbp = round( fhiKb/(fs/512) ) + 1;
Dawn@4: %   
Dawn@4: %             if (b-9 >= 0)            
Dawn@4: %                 groupSize = 2^ceil( (b-8)/4);
Dawn@4: %                 hiKbp = round((hiKbp-lowKbp+1)/groupSize)*groupSize + lowKbp -1 ;
Dawn@4: %             else
Dawn@4: %                 groupSize = 1;
Dawn@4: %             end
Dawn@4: %         
Dawn@4: %             newtemp = PSD(lowKbp:hiKbp).*conj(PSD(lowKbp:hiKbp));
Dawn@4: %             
Dawn@4: %             if (b-9 >= 0)                
Dawn@4: %                 temp1 = newtemp(1:groupSize:(hiKbp-lowKbp+1)) ;
Dawn@4: %                 for i=2:groupSize
Dawn@4: %                 temp1 = temp1 + newtemp(i:groupSize:(hiKbp-lowKbp+1)) ;
Dawn@4: %                 end
Dawn@4: %                 newtemp = temp1;
Dawn@4: %             end        
Dawn@4: %             GM(b,n) = exp( sum(log(newtemp))/ (hiKbp-lowKbp+1));
Dawn@4: %             AM(b,n) = sum(newtemp) / (hiKbp-lowKbp+1);
Dawn@4: %         end
Dawn@4: %             ASF = (GM./AM)';
Dawn@4: %             subplot(312); plot(temp);subplot(313); plot(ASF(n,1:6));
Dawn@4: % end
Dawn@4: 
Dawn@4: 
Dawn@4: