Accompanying MATLAB code for https://github.com/project-spectra/test-app/wiki/Formants-and-their-extraction

formants_and_their_extraction_code.m

[snd, originalFs] = audioread('pooling_transfusion_oo_0_95_1_0.wav');
Fs = 10000;

% downsample to 10 KHz
snd = resample(snd, Fs, originalFs);

plot(snd);
figure;
segmentlen = 100;
noverlap = 90;
NFFT = 4096;

% spectrogram(snd,segmentlen,0,NFFT,Fs,'yaxis');

% https://dsp.stackexchange.com/questions/36346/how-can-i-plot-the-spectrum-of-a-signal-in-matlab
spectrum = 10*log(abs(fftshift(fft(snd))) / length(snd));   %compute the FFT
precision = Fs/length(snd);
f = linspace(-Fs/2+precision/2, Fs/2-precision/2, length(snd)); % Create the frequency axis and put the measure in the middle of the bin.
plot(f,spectrum);

excerpt_audio_player = audioplayer(snd, Fs);

% warning: preEmphasisFrequency must be below the nyquist frequency!
% http://www.fon.hum.uva.nl/praat/manual/Sound__Filter__pre-emphasis____.html
preEmphasisFrequency = 50; % the default in praat, hz
preEmphasized = zeros(size(snd));
alpha = exp(-2.0 * pi * preEmphasisFrequency * (1/Fs) ); 

for i = length(preEmphasized):-1:2
    preEmphasized(i) = snd(i) - alpha * snd(i-1);
end

snd = preEmphasized;

snd = snd.*gausswin(length(snd));

[A, err] = lpc(snd,10);

[H, F] = freqz(err^0.5,A,NFFT,Fs);
figure;
plot(F,log(abs(H)));

rts = roots(A);

rts = rts(imag(rts)>=0); % retain only the roots with positive imaginary component.
angz = atan2(imag(rts),real(rts));

[frqs,indices] = sort(angz.*(Fs/(2*pi)));
bw = -1/2*(Fs/(2*pi))*log(abs(rts(indices)));

nn = 1;
for kk = 1:length(frqs)
    if (frqs(kk) > 90 && bw(kk) <400)
        formants(nn) = frqs(kk);
        nn = nn+1;
    end
end