annotate Code/Descriptors/Matlab/MPEG7/FromWeb/VoiceSauce/func_GetHNR.m @ 4:92ca03a8fa99 tip

Update to ICASSP 2013 benchmark

function [HNR05, HNR15, HNR25, HNR35] = func_GetHNR(y, Fs, F0, variables)
% N = func_GetHNR(y, Fs, F0, variables, textgridfile
% Input:  y, Fs - from wavread
%         F0 - vector of fundamental frequencies
%         variables - global settings
%         textgridfile - this is optional
% Output: N vectors
% Notes:  Calculates the Harmonic to noise ratio based on the method
% described in de Krom, 1993 - A cepstrum-based technique for determining a
% harmonic-to-noise ratio in speech signals, JSHR.
%
% Author: Yen-Liang Shue, Speech Processing and Auditory Perception Laboratory, UCLA
% Copyright UCLA SPAPL 2009

N_periods = variables.Nperiods_EC;
sampleshift = (Fs / 1000 * variables.frameshift);

HNR05 = zeros(length(F0), 1) * NaN;
HNR15 = zeros(length(F0), 1) * NaN;
HNR25 = zeros(length(F0), 1) * NaN;
HNR35 = zeros(length(F0), 1) * NaN;


for k=1:length(F0)
    ks = round(k * sampleshift);

    if (ks <= 0 || ks > length(y))
        continue;
    end

    F0_curr = F0(k);
    N0_curr = 1 / F0_curr * Fs;

    if (isnan(F0_curr))
        continue;
    end

    ystart = round(ks - N_periods/2*N0_curr);
    yend = round(ks + N_periods/2*N0_curr) - 1;

    if (mod(yend - ystart + 1, 2) == 0)  % length is even, make it odd
        yend = yend - 1;
    end

    if (ystart <= 0)
        continue;
    end;

    if (yend > length(y))
        continue;
    end;

    yseg = y(ystart:yend);


    hnr = getHNR(yseg, F0_curr, Fs, [500 1500 2500 3500]);

    HNR05(k) = hnr(1);
    HNR15(k) = hnr(2);
    HNR25(k) = hnr(3);
    HNR35(k) = hnr(4);

end


% main processing function
function n = getHNR(y, F0, Fs, Nfreqs)

NBins = length(y);
N0 = round(Fs / F0);
N0_delta = round(N0 * 0.1);  % search 10% either side

y = y .* hamming(length(y));
Y = fft(y, NBins);
aY = log10(abs(Y));
ay = ifft(aY);

% find possible rahmonic peaks
peakinx = zeros(floor(length(y)/ 2 / N0), 1);
for k=1:length(peakinx)
    ayseg = ay(k*N0 - N0_delta : k*N0 + N0_delta);
    [val, inx] = max(abs(ayseg));
    peakinx(k) = inx + k*N0 - N0_delta - 1;

    % lifter out each rahmonic peak
    s_ayseg = sign(diff(ayseg));
    
    % find first derivative sign change
    l_inx = inx - find((s_ayseg(inx-1:-1:1) ~= sign(inx)), 1) + 1;
    r_inx = inx + find((s_ayseg(inx+1:end) == sign(inx)), 1);
    
    l_inx = l_inx + k*N0 - N0_delta - 1;
    r_inx = r_inx + k*N0 - N0_delta - 1;
    
    % lifter out the peak
    ay(l_inx:r_inx) = 0;
   
end

% put the signal back together
midL = round(length(y) / 2) + 1;
ay(midL : end) = ay(midL - 1: -1 : midL - 1 - (length(ay) - midL));

Nap = real(fft(ay));
N = Nap;
Ha = aY - Nap;  % approximated harmonic spectrum

% calculate baseline corrections
Hdelta = F0 / Fs * length(y);
for f=Hdelta+0.0001:Hdelta:round(length(y)/2)
    fstart = ceil(f - Hdelta);
    
    Bdf = abs(min(Ha(fstart:round(f))));
    N(fstart:round(f)) = N(fstart:round(f)) - Bdf;
end

% calculate the average HNR for Nfreqs
H = aY - N;  %note that N is valid only for 1:length(N)/2
n = zeros(length(Nfreqs), 1);
for k=1:length(Nfreqs)
    Ef = round(Nfreqs(k) / Fs * length(y));  % equivalent cut off frequency
    n(k) = 20*mean(H(1:Ef)) - 20*mean(N(1:Ef));
end