Dawn@4: function [HNR05, HNR15, HNR25, HNR35] = func_GetHNR(y, Fs, F0, variables)
Dawn@4: % N = func_GetHNR(y, Fs, F0, variables, textgridfile
Dawn@4: % Input:  y, Fs - from wavread
Dawn@4: %         F0 - vector of fundamental frequencies
Dawn@4: %         variables - global settings
Dawn@4: %         textgridfile - this is optional
Dawn@4: % Output: N vectors
Dawn@4: % Notes:  Calculates the Harmonic to noise ratio based on the method
Dawn@4: % described in de Krom, 1993 - A cepstrum-based technique for determining a
Dawn@4: % harmonic-to-noise ratio in speech signals, JSHR.
Dawn@4: %
Dawn@4: % Author: Yen-Liang Shue, Speech Processing and Auditory Perception Laboratory, UCLA
Dawn@4: % Copyright UCLA SPAPL 2009
Dawn@4: 
Dawn@4: N_periods = variables.Nperiods_EC;
Dawn@4: sampleshift = (Fs / 1000 * variables.frameshift);
Dawn@4: 
Dawn@4: HNR05 = zeros(length(F0), 1) * NaN;
Dawn@4: HNR15 = zeros(length(F0), 1) * NaN;
Dawn@4: HNR25 = zeros(length(F0), 1) * NaN;
Dawn@4: HNR35 = zeros(length(F0), 1) * NaN;
Dawn@4: 
Dawn@4: 
Dawn@4: for k=1:length(F0)
Dawn@4:     ks = round(k * sampleshift);
Dawn@4: 
Dawn@4:     if (ks <= 0 || ks > length(y))
Dawn@4:         continue;
Dawn@4:     end
Dawn@4: 
Dawn@4:     F0_curr = F0(k);
Dawn@4:     N0_curr = 1 / F0_curr * Fs;
Dawn@4: 
Dawn@4:     if (isnan(F0_curr))
Dawn@4:         continue;
Dawn@4:     end
Dawn@4: 
Dawn@4:     ystart = round(ks - N_periods/2*N0_curr);
Dawn@4:     yend = round(ks + N_periods/2*N0_curr) - 1;
Dawn@4: 
Dawn@4:     if (mod(yend - ystart + 1, 2) == 0)  % length is even, make it odd
Dawn@4:         yend = yend - 1;
Dawn@4:     end
Dawn@4: 
Dawn@4:     if (ystart <= 0)
Dawn@4:         continue;
Dawn@4:     end;
Dawn@4: 
Dawn@4:     if (yend > length(y))
Dawn@4:         continue;
Dawn@4:     end;
Dawn@4: 
Dawn@4:     yseg = y(ystart:yend);
Dawn@4: 
Dawn@4: 
Dawn@4:     hnr = getHNR(yseg, F0_curr, Fs, [500 1500 2500 3500]);
Dawn@4: 
Dawn@4:     HNR05(k) = hnr(1);
Dawn@4:     HNR15(k) = hnr(2);
Dawn@4:     HNR25(k) = hnr(3);
Dawn@4:     HNR35(k) = hnr(4);
Dawn@4: 
Dawn@4: end
Dawn@4: 
Dawn@4: 
Dawn@4: % main processing function
Dawn@4: function n = getHNR(y, F0, Fs, Nfreqs)
Dawn@4: 
Dawn@4: NBins = length(y);
Dawn@4: N0 = round(Fs / F0);
Dawn@4: N0_delta = round(N0 * 0.1);  % search 10% either side
Dawn@4: 
Dawn@4: y = y .* hamming(length(y));
Dawn@4: Y = fft(y, NBins);
Dawn@4: aY = log10(abs(Y));
Dawn@4: ay = ifft(aY);
Dawn@4: 
Dawn@4: % find possible rahmonic peaks
Dawn@4: peakinx = zeros(floor(length(y)/ 2 / N0), 1);
Dawn@4: for k=1:length(peakinx)
Dawn@4:     ayseg = ay(k*N0 - N0_delta : k*N0 + N0_delta);
Dawn@4:     [val, inx] = max(abs(ayseg));
Dawn@4:     peakinx(k) = inx + k*N0 - N0_delta - 1;
Dawn@4: 
Dawn@4:     % lifter out each rahmonic peak
Dawn@4:     s_ayseg = sign(diff(ayseg));
Dawn@4:     
Dawn@4:     % find first derivative sign change
Dawn@4:     l_inx = inx - find((s_ayseg(inx-1:-1:1) ~= sign(inx)), 1) + 1;
Dawn@4:     r_inx = inx + find((s_ayseg(inx+1:end) == sign(inx)), 1);
Dawn@4:     
Dawn@4:     l_inx = l_inx + k*N0 - N0_delta - 1;
Dawn@4:     r_inx = r_inx + k*N0 - N0_delta - 1;
Dawn@4:     
Dawn@4:     % lifter out the peak
Dawn@4:     ay(l_inx:r_inx) = 0;
Dawn@4:    
Dawn@4: end
Dawn@4: 
Dawn@4: % put the signal back together
Dawn@4: midL = round(length(y) / 2) + 1;
Dawn@4: ay(midL : end) = ay(midL - 1: -1 : midL - 1 - (length(ay) - midL));
Dawn@4: 
Dawn@4: Nap = real(fft(ay));
Dawn@4: N = Nap;
Dawn@4: Ha = aY - Nap;  % approximated harmonic spectrum
Dawn@4: 
Dawn@4: % calculate baseline corrections
Dawn@4: Hdelta = F0 / Fs * length(y);
Dawn@4: for f=Hdelta+0.0001:Hdelta:round(length(y)/2)
Dawn@4:     fstart = ceil(f - Hdelta);
Dawn@4:     
Dawn@4:     Bdf = abs(min(Ha(fstart:round(f))));
Dawn@4:     N(fstart:round(f)) = N(fstart:round(f)) - Bdf;
Dawn@4: end
Dawn@4: 
Dawn@4: % calculate the average HNR for Nfreqs
Dawn@4: H = aY - N;  %note that N is valid only for 1:length(N)/2
Dawn@4: n = zeros(length(Nfreqs), 1);
Dawn@4: for k=1:length(Nfreqs)
Dawn@4:     Ef = round(Nfreqs(k) / Fs * length(y));  % equivalent cut off frequency
Dawn@4:     n(k) = 20*mean(H(1:Ef)) - 20*mean(N(1:Ef));
Dawn@4: end
Dawn@4: 
Dawn@4: 
Dawn@4: 
Dawn@4: