function formant_features = calc_formants(in,fs)
%function formant_features = calc_formants(in,fs)
%
% CALC_FORMANTS  --  Estimates the frequencies (in Hz) of the two most prominent
%                    formants in an vowel signal by calculating the
%                    smoothed spectrum and finding the most prominent
%                    peak both above and below a cutoff frequency of 1500 Hz
%
% Input Parameters:
%   in           --  The input vowel signal.  If 'in' is a matrix, the function 
%                    assumes one signal per row and estimates formant frequencies
%                    for each signal.
%   fs           --  Sampling frequency of the input signal(s)
%
% Output Parameters:
%   formant_features -- A matrix of formant feature vectors with one row per
%                    row of the input matrix, 'in'.  The first column contains the
%                    frequency of the most prominent peak below 1500 Hz; the
%                    second column is the frequency of the most prominent peak above 
%                    1500 Hz

if nargin < 2
    fs = 8192;
end

fft_size = 2048;
frqs = (0:fft_size-1)*fs/fft_size;

% Define the impulse response for spectral smoothing
h1 = hanning(100)'.^2;

% Define the cutoff frequency
cutoff = 1500;                      % In Hertz
cutoff = round(cutoff/fs*fft_size);

% Make sure the input is oriented properly
if any(size(in)==1)
    in = in(:)';
end

for i=1:size(in,1);
    % Calculate the spectrum
    spec = abs(fft(in(i,:),fft_size));
    sm_spec = 20*log10(filter2(h1,spec));

    % Find the two formant peaks
    i1 = find_peaks(sm_spec(1:cutoff));
    i2 = find_peaks(sm_spec(cutoff:end/2));
    formant_features(i,:) = [frqs(i1) frqs(i2+cutoff)];
end



function out = find_peaks(in);
% Returns the index of the largest peak (not including edges) in 'out'
% (A peak is defined as any sample that is higher than both of its neighbors)

% Find all peaks
ind = find(in(2:end-1) > in(1:end-2) & in(2:end-1) > in(3:end)) + 1;

% Find peak with maximum amplitude
[a,i] = max(in(ind));

% Return the index of the largest peak
out = ind(i);
