Dawn@4: function [A1, A2, A3] = func_GetA1A2A3(y, Fs, F0, F1, F2, F3, variables, textgridfile)
Dawn@4: % [A1, A2, A3] = func_GetA1A2A3(y, Fs, F0, F1, F2, F3, variables, textgridfile)
Dawn@4: % Input:  y - input signal
Dawn@4: %         Fs - sampling frequency
Dawn@4: %         F0 - fundamental frequency (vector)
Dawn@4: %         F1-F3 - formant frequencies (vectors)
Dawn@4: %         variables - settings from VS
Dawn@4: %         textgridfile - textgrid filename if available (optional)
Dawn@4: % Output: A1, A2, A3 - uncorrected amplitudes of the spectra at the formant
Dawn@4: % frequencies.
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;
Dawn@4: sampleshift = (Fs / 1000 * variables.frameshift);
Dawn@4: 
Dawn@4: A1 = zeros(length(F0), 1) * NaN;
Dawn@4: A2 = zeros(length(F0), 1) * NaN;
Dawn@4: A3 = zeros(length(F0), 1) * NaN;
Dawn@4: 
Dawn@4: if (nargin == 7)
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:         
Dawn@4:         if (isnan(F0_curr) || F0_curr == 0) 
Dawn@4:             continue; 
Dawn@4:         end
Dawn@4: 
Dawn@4:         N0_curr = 1 / F0_curr * Fs;
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 (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:         % proceed if no NaNs detected
Dawn@4:         if (~isnan(F1(k)) && ~isnan(F2(k)) && ~isnan(F3(k)) && ~isnan(N0_curr))
Dawn@4:             [A1_, fmax] = ana_GetMagnitudeMax(yseg, F1(k), Fs, 8192);
Dawn@4:             [A2_, fmax] = ana_GetMagnitudeMax(yseg, F2(k), Fs, 8192);
Dawn@4:             [A3_, fmax] = ana_GetMagnitudeMax(yseg, F3(k), Fs, 8192);
Dawn@4: 
Dawn@4:             A1(k) = A1_;
Dawn@4:             A2(k) = A2_;
Dawn@4:             A3(k) = A3_;
Dawn@4:         end    
Dawn@4:     end
Dawn@4: else
Dawn@4:     % use textgrids - get ignore string list
Dawn@4:     tbuffer = variables.tbuffer;
Dawn@4:     ignorelabels = textscan(variables.TextgridIgnoreList, '%s', 'delimiter', ',');
Dawn@4:     ignorelabels = ignorelabels{1};   
Dawn@4:     
Dawn@4:     [labels, start, stop] = func_readTextgrid(textgridfile);
Dawn@4:     labels_tmp = [];
Dawn@4:     start_tmp = [];
Dawn@4:     stop_tmp = [];
Dawn@4:     
Dawn@4:     for k=1:length(variables.TextgridTierNumber)
Dawn@4:         inx = variables.TextgridTierNumber(k);
Dawn@4:         if (inx <= length(labels))
Dawn@4:             labels_tmp = [labels_tmp; labels{inx}];
Dawn@4:             start_tmp = [start_tmp; start{inx}];
Dawn@4:             stop_tmp = [stop_tmp; stop{inx}];
Dawn@4:         end
Dawn@4:     end
Dawn@4:     
Dawn@4:     labels = labels_tmp;
Dawn@4:     start = start_tmp;
Dawn@4:     stop = stop_tmp; 
Dawn@4:     
Dawn@4:     
Dawn@4:     for m=1:length(start)
Dawn@4:         switch(labels{m})
Dawn@4:             case ignorelabels
Dawn@4:                 continue;  % skip anything that is within the ignore list
Dawn@4:         end
Dawn@4:         
Dawn@4:         kstart = round((start(m) * 1000 - tbuffer) / variables.frameshift);
Dawn@4:         kstop = round((stop(m) * 1000 + tbuffer) / variables.frameshift);
Dawn@4:         kstart(kstart <= 0) = 1;
Dawn@4:         kstop(kstop > length(F0)) = length(F0);
Dawn@4:         
Dawn@4:         ystart = round(kstart * sampleshift);
Dawn@4:         ystop = round(kstop * sampleshift);
Dawn@4:         ystart(ystart <= 0) = 1;
Dawn@4:         ystop(ystop > length(y)) = length(y);
Dawn@4:         
Dawn@4:         ysegment = y(ystart:ystop);
Dawn@4:         F0seg = F0(kstart:kstop);
Dawn@4:         F1seg = F1(kstart:kstop);
Dawn@4:         F2seg = F2(kstart:kstop);
Dawn@4:         F3seg = F3(kstart:kstop);
Dawn@4:                
Dawn@4:         for k=1:length(F0seg)
Dawn@4:             ks = round(k*sampleshift);
Dawn@4:             
Dawn@4:             if (ks <= 0 || ks > length(ysegment))
Dawn@4:                 continue;
Dawn@4:             end
Dawn@4:             
Dawn@4:             F0_curr = F0seg(k);
Dawn@4:             if (isnan(F0_curr) || F0_curr == 0)
Dawn@4:                 continue;
Dawn@4:             end
Dawn@4:             
Dawn@4:             N0_curr = Fs/F0_curr;
Dawn@4:             ysstart = round(ks - N_periods/2 * N0_curr);
Dawn@4:             ysend = round(ks + N_periods/2 * N0_curr) - 1;
Dawn@4:             
Dawn@4:             if (isnan(F1seg(k)) || isnan(F2seg(k)) || isnan(F3seg(k)))
Dawn@4:                 continue;
Dawn@4:             end
Dawn@4:             
Dawn@4:             if (ysstart <= 0) 
Dawn@4:                 continue; 
Dawn@4:             end;
Dawn@4:             if (ysend > length(ysegment)) 
Dawn@4:                 continue; 
Dawn@4:             end;
Dawn@4:             
Dawn@4:             yseg = ysegment(ysstart:ysend);
Dawn@4:             
Dawn@4:             [A1_, fmax] = ana_GetMagnitudeMax(yseg, F1seg(k), Fs, 8192);
Dawn@4:             [A2_, fmax] = ana_GetMagnitudeMax(yseg, F2seg(k), Fs, 8192);
Dawn@4:             [A3_, fmax] = ana_GetMagnitudeMax(yseg, F3seg(k), Fs, 8192);
Dawn@4:             
Dawn@4:             inx = kstart + k - 1;
Dawn@4:             if (inx <= length(A1))
Dawn@4:                 A1(inx) = A1_;
Dawn@4:                 A2(inx) = A2_;
Dawn@4:                 A3(inx) = A3_;
Dawn@4:             end                
Dawn@4:         
Dawn@4:         end
Dawn@4:     
Dawn@4:     end
Dawn@4: end
Dawn@4: 
Dawn@4: 
Dawn@4: %------------------------------------------------------------------------
Dawn@4: function [M,fmax] = ana_GetMagnitudeMax(x, Fx, Fs, fftlen)
Dawn@4: % Get maximal spectral magnitude of signal x around position Fx Hz in dB
Dawn@4: % Fx can be a vector of frequency points
Dawn@4: % Note that the bigger fftlen the better the resolution
Dawn@4: 
Dawn@4: if (isnan(Fx))
Dawn@4:   M = NaN;
Dawn@4:   fmax = NaN;
Dawn@4: else
Dawn@4: 
Dawn@4:   len = length(x);
Dawn@4:   hamlen = min(fftlen, len);
Dawn@4:   %X = fft(hamming(hamlen).*x(1:hamlen), fftlen);
Dawn@4:   factor = 1;  %/length(x);  % not exactly Kosher but compatible to dfs_magn()
Dawn@4:   X = fft(x,fftlen);
Dawn@4:   X(X==0) = 0.000000001;  % guard against log(0) warnings
Dawn@4:   X = 20*log10(factor*abs(X(1:fftlen/2, :)));
Dawn@4:   fstep = Fs / fftlen;
Dawn@4:   lowf = Fx-0.1*Fx;
Dawn@4:   if (lowf<0) lowf=0; end
Dawn@4:   highf = Fx+0.1*Fx;
Dawn@4:   if (highf>Fs/2-fstep) highf=Fs/2-fstep; end
Dawn@4: 
Dawn@4:   for cnt = 1:length(Fx)
Dawn@4:       [M(cnt),pos] = max(X(1+round(lowf(cnt)/fstep):1+round(highf(cnt)/fstep), :));
Dawn@4:       fmax(cnt) = (pos-1+round(lowf(cnt)/fstep))*fstep;
Dawn@4:   end
Dawn@4:   
Dawn@4: end
Dawn@4: 
Dawn@4: % f=0:fstep:(Fs/2-fstep);
Dawn@4: % plot(f,X); hold;
Dawn@4: % plot(pos, M, 'rx'); hold off;
Dawn@4: % fprintf('End get_magnitude_max()\n');