comparison Code/Descriptors/Matlab/MPEG7/FromWeb/VoiceSauce/shrp.m @ 4:92ca03a8fa99 tip

Update to ICASSP 2013 benchmark

+function [f0_time,f0_value,SHR,f0_candidates]=shrp(Y,Fs,F0MinMax,frame_length,timestep,SHR_Threshold,ceiling,med_smooth,CHECK_VOICING)
+% SHRP - a pitch determination algorithm based on Subharmonic-to-Harmonic Ratio (SHR)
+% [f0_time,f0_value,SHR,f0_candidates]=shrp(Y,Fs[,F0MinMax,frame_length,TimeStep,SHR_Threshold,Ceiling,med_smooth,CHECK_VOICING]) 
+%
+%   Input parameters (There are 9):
+%
+%       Y:              Input data 
+%       Fs:             Sampling frequency (e.g., 16000 Hz)
+%       F0MinMax:       2-d array specifies the F0 range. [minf0 maxf0], default: [50 550]
+%                       Quick solutions:
+%                       For male speech: [50 250]
+%                       For female speech: [120 400]
+%       frame_length:   length of each frame in millisecond (default: 40 ms)
+%       TimeStep:       Interval for updating short-term analysis in millisecond (default: 10 ms)
+%       SHR_Threshold:  Subharmonic-to-harmonic ratio threshold in the range of [0,1] (default: 0.4). 
+%                       If the estimated SHR is greater than the threshold, the subharmonic is regarded as F0 candidate,
+%                       Otherwise, the harmonic is favored.
+%       Ceiling:        Upper bound of the frequencies that are used for estimating pitch. (default: 1250 Hz)       
+%       med_smooth:     the order of the median smoothing (default: 0 - no smoothing);                       
+%       CHECK_VOICING:  check voicing. Current voicing determination algorithm is kind of crude.
+%                       0: no voicing checking (default)
+%                       1: voicing checking
+%   Output parameters:
+%       
+%       f0_time:        an array stores the times for the F0 points
+%       f0_value:       an array stores F0 values
+%       SHR:            an array stores subharmonic-to-harmonic ratio for each frame
+%		f0_candidates:  a matrix stores the f0 candidates for each frames, currently two f0 values generated for each frame.
+%						Each row (a frame) contains two values in increasing order, i.e., [low_f0 higher_f0].
+%						For SHR=0, the first f0 is 0. The purpose of this is that when you want to test different SHR
+%						thresholds, you don't need to re-run the whole algorithm. You can choose to select the lower or higher
+%						value based on the shr value of this frame. 	
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%	Permission to use, copy, modify, and distribute this software without fee is hereby granted
+%	FOR RESEARCH PURPOSES only, provided that this copyright notice appears in all copies 
+% 	and in all supporting documentation.
+%
+%   This program is distributed in the hope that it will be useful,
+%   but WITHOUT ANY WARRANTY; without even the implied warranty of
+%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+%   
+%   For details of the algorithm, please see 
+%   Sun, X.,"Pitch determination and voice quality analysis using subharmonic-to-harmonic ratio" To appear in the Proc. of ICASSP2002, Orlando, Florida, May 13 -17, 2002.
+%   For update information, please check http://mel.speech.nwu.edu/sunxj/pda.htm.
+%
+%   Copyright (c) 2001 Xuejing Sun
+%   Department of Communication Sciences and Disorders
+%   Northwestern University, USA
+%   [e-mail]
+%
+%   Update history: 
+%	Added "f0_candidates" as a return value, Dec. 21, 2001
+%	Changed default median smoothing order from 5 to 0, Jan. 9, 2002
+%   Modified the GetLogSpectrum function, bug fixed due to Herbert Griebel. Jan. 15, 2002
+%   Several minor changes. Jan. 15,2002.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%t0 = clock;      
+%------------------ Get input arguments values and set default values -------------------------
+if nargin<9
+    CHECK_VOICING=0;
+end
+if nargin<8
+    med_smooth=0;
+end
+if nargin<7
+    ceiling=1250;
+end
+if nargin<6
+    SHR_Threshold=0.4; % subharmonic to harmonic ratio threshold
+end
+if nargin<5
+    timestep=10;
+   %timestep=6.4;
+end
+if nargin<4
+    frame_length=40; % default 40 ms
+end
+if nargin<3
+    minf0=50;
+    maxf0=500;
+else
+    minf0=F0MinMax(1);
+    maxf0=F0MinMax(2);
+end
+if nargin<2
+    error('Sampling rate must be supplied!')
+end
+segmentduration=frame_length; 
+
+%------------------- pre-processing input signal -------------------------
+Y=Y-mean(Y); % remove DC component
+Y=Y/max(abs(Y));  %normalization
+total_len=length(Y);
+%------------------ specify some algorithm-specific thresholds -------------------------
+interpolation_depth=0.5;  % for FFT length
+%--------------- derived thresholds specific to the algorithm -------------------------------
+maxlogf=log2(maxf0/2);
+minlogf=log2(minf0/2); % the search region to compute SHR is as low as 0.5 minf0.
+N=floor(ceiling/minf0); % maximum number harmonics
+m=mod(N,2);
+N=N-m;
+N=N*4; %In fact, in most cases we don't need to multiply N by 4 and get equally good results yet much faster.
+% derive how many frames we have based on segment length and timestep.
+segmentlen=round(segmentduration*(Fs/1000));
+inc=round(timestep*(Fs/1000));
+nf = fix((total_len-segmentlen+inc)/inc);
+n=(1:nf);
+f0_time=((n-1)*timestep+segmentduration/2)'; % anchor time for each frame, the middle point
+%f0_time=((n-1)*timestep)'; % anchor time for each frame, starting from zero
+%------------------ determine FFT length ---------------------
+fftlen=1;
+while (fftlen < segmentlen * (1 +interpolation_depth)) 
+    fftlen =fftlen* 2;
+end
+%----------------- derive linear and log frequency scale ----------------
+frequency=Fs*(1:fftlen/2)/fftlen; % we ignore frequency 0 here since we need to do log transformation later and won't use it anyway.
+limit=find(frequency>=ceiling);
+limit=limit(1); % only the first is useful
+frequency=frequency(1:limit);
+logf=log2(frequency);
+%% clear some variables to save memory
+clear frequency;
+min_bin=logf(end)-logf(end-1); % the minimum distance between two points after interpolation
+shift=log2(N); % shift distance
+shift_units=round(shift/min_bin); %the number of unit on the log x-axis
+i=(2:N);
+% ------------- the followings are universal for all the frames ---------------%%
+startpos=shift_units+1-round(log2(i)/min_bin);  % find out all the start position of each shift
+index=find(startpos<1); % find out those positions that are less than 1
+startpos(index)=1; % set them to 1 since the array index starts from 1 in matlab
+interp_logf=logf(1):min_bin:logf(end);
+interp_len=length(interp_logf);% new length of the amplitude spectrum after interpolation
+totallen=shift_units+interp_len;
+endpos=startpos+interp_len-1; %% note that : totallen=shift_units+interp_len;
+index=find(endpos>totallen);
+endpos(index)=totallen; % make sure all the end positions not greater than the totoal length of the shift spectrum
+
+newfre=2.^(interp_logf); % the linear Hz scale derived from the interpolated log scale
+upperbound=find(interp_logf>=maxlogf); % find out the index of upper bound of search region on the log frequency scale.
+upperbound=upperbound(1);% only the first element is useful
+lowerbound=find(interp_logf>=minlogf); % find out the index of lower bound of search region on the log frequency scale.
+lowerbound=lowerbound(1);
+
+%----------------- segmentation of speech ------------------------------
+curpos=round(f0_time/1000*Fs);   % position for each frame in terms of index, not time
+frames=toframes(Y,curpos,segmentlen,'hamm');
+[nf framelen]=size(frames);
+clear Y;
+%----------------- initialize vectors for f0 time, f0 values, and SHR
+f0_value=zeros(nf,1);
+SHR=zeros(nf,1);
+f0_time=f0_time(1:nf);
+f0_candidates=zeros(nf,2);
+%----------------- voicing determination ----------------------------
+if (CHECK_VOICING)
+    NoiseFloor=sum(frames(1,:).^2);   
+    voicing=vda(frames,segmentduration/1000,NoiseFloor);
+else
+    voicing=ones(nf,1);
+end
+%------------------- the main loop -----------------------
+curf0=0;
+cur_SHR=0;
+cur_cand1=0;
+cur_cand2=0;
+for n=1:nf
+    segment=frames(n,:);
+    curtime=f0_time(n);
+    if voicing(n)==0
+        curf0=0;
+        cur_SHR=0;
+    else
+        [log_spectrum]=GetLogSpectrum(segment,fftlen,limit,logf,interp_logf);
+        [peak_index,cur_SHR,shshift,all_peak_indices]=ComputeSHR(log_spectrum,min_bin,startpos,endpos,lowerbound,upperbound,N,shift_units,SHR_Threshold);
+        if (peak_index==-1) % -1 indicates a possibly unvoiced frame, if CHECK_VOICING, set f0 to 0, otherwise uses previous value
+            if (CHECK_VOICING)
+                curf0=0;
+                cur_cand1=0;
+                cur_cand2=0;
+            end
+            
+        else
+            curf0=newfre(peak_index)*2;  
+            if (curf0>maxf0)
+                curf0=curf0/2;
+            end
+            if (length(all_peak_indices)==1)
+            	cur_cand1=0;
+            	cur_cand2=newfre(all_peak_indices(1))*2;
+            else
+            	cur_cand1=newfre(all_peak_indices(1))*2;
+            	cur_cand2=newfre(all_peak_indices(2))*2;
+            end	
+            if (cur_cand1>maxf0)
+                cur_cand1=cur_cand1/2;
+            end
+            if (cur_cand2>maxf0)
+                cur_cand2=cur_cand2/2;
+            end
+            if (CHECK_VOICING)
+                voicing(n)=postvda(segment,curf0,Fs);
+                if (voicing(n)==0)
+                    curf0=0;
+                end
+            end
+        end
+    end
+    f0_value(n)=curf0;
+    SHR(n)=cur_SHR;
+    f0_candidates(n,1)=cur_cand1;
+    f0_candidates(n,2)=cur_cand2;
+    DEBUG=0;
+    if DEBUG
+        figure(9)
+        %subplot(5,1,1),plot(segment,'*') 
+        %title('windowed waveform segment')
+        subplot(2,2,1),plot(interp_logf,log_spectrum,'k*')
+        title('(a)')
+        grid
+        %('spectrum on log frequency scale')
+        %grid
+        shsodd=sum(shshift(1:2:N-1,:),1); 
+        shseven=sum(shshift(2:2:N,:),1);
+        difference=shseven-shsodd;
+        subplot(2,2,2),plot(interp_logf,shseven,'k*')
+        title('(b)')
+        %title('even')
+        grid
+        subplot(2,2,3),plot(interp_logf,shsodd,'k*')
+        title('(c)')
+        %title('odd')
+        grid
+        subplot(2,2,4), plot(interp_logf,difference,'k*')
+        title('(d)')
+        %title('difference (even-odd)')   
+        grid
+        curtime
+        curf0
+        cur_SHR
+        pause
+    end
+end
+%-------------- post-processing -------------------------------
+if (med_smooth > 0)
+    f0_value=medsmooth(f0_value,med_smooth);
+end
+%f0=linsmooth(f0,5); % this is really optional.
+
+%*****************************************************************************************
+%-------------- do FFT and get log spectrum  ---------------------------------
+%*****************************************************************************************
+function [interp_amplitude]=GetLogSpectrum(segment,fftlen,limit,logf,interp_logf)
+Spectra=fft(segment,fftlen);
+amplitude = abs(Spectra(1:fftlen/2+1)); % fftlen is always even here. Note: change fftlen/2 to fftlen/2+1. bug fixed due to Herbert Griebel
+amplitude=amplitude(2:limit+1); % ignore the zero frequency component
+%amplitude=log10(amplitude+1);
+interp_amplitude=interp1(logf,amplitude,interp_logf,'linear');
+interp_amplitude=interp_amplitude-min(interp_amplitude);
+%*****************************************************************************************
+%-------------- compute subharmonic-to-harmonic ratio  ---------------------------------
+%*****************************************************************************************
+function [peak_index,SHR,shshift,index]=ComputeSHR(log_spectrum,min_bin,startpos,endpos,lowerbound,upperbound,N,shift_units,SHR_Threshold)
+% computeshr: compute subharmonic-to-harmonic ratio for a short-term signal
+len_spectrum=length(log_spectrum);
+totallen=shift_units+len_spectrum;
+shshift=zeros(N,totallen); %initialize the subharmonic shift matrix; each row corresponds to a shift version
+shshift(1,(totallen-len_spectrum+1):totallen)=log_spectrum; % place the spectrum at the right end of the first row
+% note that here startpos and endpos has N-1 rows, so we start from 2
+% the first row in shshift is the original log spectrum
+for i=2:N
+    shshift(i,startpos(i-1):endpos(i-1))=log_spectrum(1:endpos(i-1)-startpos(i-1)+1); % store each shifted sequence
+end
+shshift=shshift(:,shift_units+1:totallen); % we don't need the stuff smaller than shift_units
+shsodd=sum(shshift(1:2:N-1,:),1); 
+shseven=sum(shshift(2:2:N,:),1);
+difference=shseven-shsodd;
+% peak picking process
+SHR=0;
+[mag,index]=twomax(difference,lowerbound,upperbound,min_bin); % only find two maxima
+% first mag is always the maximum, the second, if there is, is the second max
+NumPitchCandidates=length(mag);
+if (NumPitchCandidates == 1) % this is possible, mainly due to we put a constraint on search region, i.e., f0 range
+    if (mag <=0) % this must be an unvoiced frame
+        peak_index=-1;
+        return
+    end
+    peak_index=index;
+    SHR=0;
+else
+    SHR=(mag(1)-mag(2))/(mag(1)+mag(2));
+    if (SHR<=SHR_Threshold) 
+        peak_index=index(2);  % subharmonic is weak, so favor the harmonic
+    else
+        peak_index=index(1); % subharmonic is strong, so favor the subharmonic as F0
+    end
+end
+%%*****************************************************************************************
+%******************    this function only finds two maximum peaks   ************************
+function [mag,index]=twomax(x,lowerbound,upperbound,unitlen)
+%In descending order, the magnitude and index are returned in [mag,index], respectively
+lenx=length(x);
+halfoct=round(1/unitlen/2); % compute the number of units of half octave. log2(2)=1; 1/unitlen
+[mag,index]=max(x(lowerbound:upperbound));%find the maximum value
+if (mag<=0)
+    %    error('max is smaller than zero!') % return it!
+    return
+end
+index=index+lowerbound-1;
+harmonics=2;
+LIMIT=0.0625; % 1/8 octave
+startpos=index+round(log2(harmonics-LIMIT)/unitlen);
+if (startpos<=min(lenx,upperbound))
+    endpos=index+round(log2(harmonics+LIMIT)/unitlen); % for example, 100hz-200hz is one octave, 200hz-250hz is 1/4octave
+    if (endpos> min(lenx,upperbound))
+        endpos=min(lenx,upperbound);
+    end
+    [mag1,index1]=max(x(startpos:endpos));%find the maximum value at right side of last maximum 
+    if (mag1>0)
+        index1=index1+startpos-1;
+        mag=[mag;mag1];
+        index=[index;index1];
+    end
+end
+%*****************************************************************************************
+%%----------------------------------------------------------------------------------------
+%%-----------------------------------voicing determination -------------------------------
+function voice=vda(x,segmentdur,noisefloor,minzcr)
+%voice=vda(x) determine whether the segment is voiced, unvoiced or silence
+%this VDA is independent from the PDA process, and does not take advantage of the info derived from PDA
+%thus, it requires more computation load.
+if nargin<4
+    %minzcr=2500; %unit: hertz
+    minzcr=3000;
+end
+if nargin<3
+    noisefloor=0.01;
+end
+[nf, len]=size(x);
+voice=ones(nf,1);
+engergy=sum(x.^2,2);
+index=find(engergy<=noisefloor*3);
+voice(index)=0;
+
+%*****************************************************************************************
+%% --------------------------------- determine the energy threshold for silence-------------------------
+function thr=ethreshold(frames)
+%%%%% use Rabiner and Sambur (1975) method
+[nf,len]=size(frames);
+lastpoint=1;
+emax=0;
+emin=0;
+e=sum(frames.^2,2);
+emax=max(e);
+emin=min(e);
+I1=0.03*(emax-emin)+emin;
+I2=4*emin;
+thr=25*min(I1,I2);
+
+%*****************************************************************************************
+%% ------------------- split signal into frames ---------------
+function frames=toframes(input,curpos,segmentlen,wintype)
+len=length(input);
+numFrames=length(curpos);
+frames=zeros(numFrames,segmentlen);
+start=curpos-round(segmentlen/2);
+offset=(0:segmentlen-1);
+index_start=find(start<1); % find out those frames beyond the first point
+start(index_start)=1; % for those, just use the first frame
+endpos=start+segmentlen-1;
+index=find(endpos>len);
+endpos(index)=len; % duplicate the last several frames if window is over the limit.
+start(index)=len+1-segmentlen;
+frames(:)=input(start(:,ones(1,segmentlen))+offset(ones(numFrames,1),:));
+[nf, len]=size(frames);
+win=window(segmentlen,wintype);
+frames = frames .* win(ones(nf,1),:);
+%*****************************************************************************************
+%-------------- post voicing checking ---------------------------------------------
+function voicing=postvda(segment, curf0,Fs,r_threshold)
+%%% check voicing again using estimated F0, which follows Hermes, SHS algorithm, JASA, 1988
+if nargin<4
+    r_threshold=0.2;
+end
+estimated_period=1/curf0;
+mid_point=round(length(segment)/2);
+num_points=round(estimated_period*Fs); % number of points in each period
+start_point=mid_point-num_points;
+end_point=mid_point+num_points;
+if (start_point <1)
+    start_point=1;
+    mid_point=start_point+num_points;
+    if (mid_point>length(segment)) % this is unreasonable, set f0 to zero
+        voicing=0;
+        return;
+    end
+end
+segment1=segment(start_point:mid_point);
+if (end_point>length(segment))
+    end_point=length(segment);
+    mid_point=end_point-num_points;
+    if (mid_point<1) % this is unreasonable, set f0 to zero
+        voicing=0;
+        return;
+    end
+end
+segment2=segment(mid_point:end_point);
+len=min(length(segment1),length(segment2));
+r=corrcoef(segment1(1:len),segment2(1:len));
+r1=r(1,2);
+if (r1<r_threshold) % correlation threshold
+    voicing=0;
+else
+    voicing=1;
+end
+USE_ZCR=1;
+if(USE_ZCR & voicing)
+    zcr1=zcr(segment1,estimated_period);
+    zcr2=zcr(segment2,estimated_period);
+    %minzcr=2500;
+    minzcr=3500;
+    if (zcr1<minzcr | zcr2<minzcr)
+        voicing=1;
+    else
+        voicing=0;
+    end
+end
+%%*****************************************************************************************
+%--------------------- Compute zero-crossing rate -------------------------------------------  
+function zcr=zcr(x,dur)
+% function zcr=zcr(x,dur) : compute zero-crossing rate
+% x: input data
+% x: duration of the input data
+[nf,len]=size(x);
+zcr=sum(0.5*abs(sign(x(:,2:len))-sign(x(:,1:len-1))))/dur;
+%%*************************************************************************************
+%--------------------- Window function -------------------------------------------  
+function w = window(N,wt,beta)
+%
+%  w = window(N,wt)
+%
+%  generate a window function
+%
+%  N = length of desired window
+%  wt = window type desired
+%       'rect' = rectangular        'tria' = triangular (Bartlett)
+%       'hann' = Hanning            'hamm'  = Hamming
+%       'blac' = Blackman
+%		  'kais' = Kaiser	
+%
+%  w = row vector containing samples of the desired window
+% beta : used in Kaiser window
+
+nn = N-1;
+n=0:nn;
+pn = 2*pi*(0:nn)/nn;
+if wt(1,1:4) == 'rect',
+    w = ones(1,N);
+elseif wt(1,1:4) == 'tria',
+    m = nn/2;
+    w = (0:m)/m;
+    w = [w w(ceil(m):-1:1)];
+elseif wt(1,1:4) == 'hann',
+    w = 0.5*(1 - cos(pn));
+elseif wt(1,1:4) == 'hamm',
+    w = .54 - .46*cos(pn);
+elseif wt(1,1:4) == 'blac',
+    w = .42 -.5*cos(pn) + .08*cos(2*pn);
+elseif wt(1,1:4) == 'kais',
+    if nargin<3
+        error('you need provide beta!')
+    end
+    w =bessel1(beta*sqrt(1-((n-N/2)/(N/2)).^2))./bessel1(beta);   
+else
+    disp('Incorrect Window type requested')
+end