view Code/Descriptors/Matlab/MPEG7/AudioSpectralFlatness.m @ 4:92ca03a8fa99 tip

Update to ICASSP 2013 benchmark

function [ASF] = AudioSpectralFlatness(raw, fs, noOfFrames, frameLength)

check this works!

hopSize = frameLength;
windowsize = frameLength;
FFTsize = frameLength;
windowType = hamming(frameLength);

% bin = fs/512;
% bandsNumber = 24;
% loEdge = 250;
% hiEdege = 16000;


loEdge = 250;
hiEdge = 16000;  % Setting default hiedge

if (hiEdge >= fs/2) % Check if value for hiedge is valid
    % sr/2-1 : Skipping extraction up to sr/2; Not possible due to 5 percent band overlap
  hiEdge =  min(hiEdge,fs/2-1); 
end
hiEdge = 2^(floor(log2(hiEdge/1000)*4)/4)*1000 ; % Setting exact value for hiEdge, rounding to next lower frequency if necessary
if (hiEdge*1.05 >= fs/2) 
    hiEdge = hiEdge/2^0.25 ; 
end; %Now it's possible to check if hiEdge is valid
	
N = frameLength;

numbands = floor(4*log2(hiEdge/loEdge));
firstband = round(log2(loEdge/1000)*4);
overlap = 0.05;
grpsize = 1;

[fftout,phase] = mpeg7getspec( raw, fs, hopSize, windowsize, windowType, FFTsize );
band=[];
for k = 1:numbands
  f_lo = loEdge * (2^((k-1)/4)) * (1-overlap);
  f_hi = loEdge * (2^((k  )/4)) * (1+overlap);
  
  i_lo = round( f_lo/(fs/N) ) + 1;
  i_hi = round( f_hi/(fs/N) ) + 1;
  
  fband(k,1) = f_lo;
  fband(k,2) = f_hi;
  
  iband(k,1) = i_lo;
  iband(k,2) = i_hi;

  % Rounding of upper index according due to coefficient grouping
  if (k+firstband-1 >= 0)                   %Start grouping at 1kHz 
    grpsize = 2^ceil( (k+firstband )/4);
    i_hi = round((i_hi-i_lo+1)/grpsize)*grpsize + i_lo-1 ;
  else
    grpsize = 1;
  end
  tmp = fftout(i_lo:i_hi,:) .^ 2;         % PSD coefficients
  ncoeffs = i_hi - i_lo + 1;
 
  if (k+firstband-1 >= 0)                   % Coefficient grouping
    tmp2 = tmp(1:grpsize:ncoeffs,:);
    for g=2:grpsize
      tmp2 = tmp2 + tmp(g:grpsize:ncoeffs,:) ;
    end
    tmp = tmp2;
  end
  % Actual calculation
  ncoeffs = ncoeffs/grpsize ;
  tmp = tmp + 1e-50;       % avoid underflow for zero signals
  gm(k,:) = exp( sum(log(tmp))/ncoeffs ); % log processing avoids overflow
  am(k,:) = sum(tmp) / ncoeffs;
end
ASF = (gm./am)';

% myColor = 'rx';
% for i=1:(length(raw)/frameLength)
%     subplot(411); plot( raw( (i-1)*frameLength+1 : (i*frameLength) ));
%     subplot(412); plot( fftout(:,i) );
%     subplot(425); plot( i, ASF(i,1), myColor ); hold on;
%     subplot(426); plot( i, ASF(i,2), myColor ); hold on;
%     subplot(427); plot( i, ASF(i,3), myColor ); hold on;
%     subplot(428); plot( i, ASF(i,4), myColor ); hold on;
% end
% lo_edge = loEdge;
% hi_edge = hiEdge;

% for n = 1:noOfFrames
%     j = (n-1) * 512 + 1;
%     k = n * 512;
%     if(k<=length(raw))
%         temp = raw(j:k);
%     end
%     if(k>length(raw))
%         
%         temp = [raw(j:end)' zeros(512-length(raw(j:end)),1)'];
%     end
%     
%         FFTX=fft(temp,hopSize);
%         subplot(311); plot(abs(FFTX(1:length(FFTX)/2)));
%         BW=1.5*fs*1024;
%         PSD=FFTX/BW;
%      
%         for b = 1:24                                
%             floKb = 0.95*loEdge * 2^(0.25*b-0.25);
%             fhiKb = 1.05*loEdge * 2^(0.25*b);
%             lowKbp = round( floKb/(fs/512) ) + 1;
%             hiKbp = round( fhiKb/(fs/512) ) + 1;
%   
%             if (b-9 >= 0)            
%                 groupSize = 2^ceil( (b-8)/4);
%                 hiKbp = round((hiKbp-lowKbp+1)/groupSize)*groupSize + lowKbp -1 ;
%             else
%                 groupSize = 1;
%             end
%         
%             newtemp = PSD(lowKbp:hiKbp).*conj(PSD(lowKbp:hiKbp));
%             
%             if (b-9 >= 0)                
%                 temp1 = newtemp(1:groupSize:(hiKbp-lowKbp+1)) ;
%                 for i=2:groupSize
%                 temp1 = temp1 + newtemp(i:groupSize:(hiKbp-lowKbp+1)) ;
%                 end
%                 newtemp = temp1;
%             end        
%             GM(b,n) = exp( sum(log(newtemp))/ (hiKbp-lowKbp+1));
%             AM(b,n) = sum(newtemp) / (hiKbp-lowKbp+1);
%         end
%             ASF = (GM./AM)';
%             subplot(312); plot(temp);subplot(313); plot(ASF(n,1:6));
% end