comparison 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
+
+
+