Finding formants in speech signal

I'm trying to find first 4 formants in signal with sliding window, but I'm getting some different results .

Code:

samples = importdata('samples.txt');
Fs=27100; 
signal_preemf = filter(1,[1 0.95],samples);  
section=round(0.03/(1/Fs));

for i=1:round(0.001/(1/Fs)):length(samples)-section
    signal_part=signal_preemf(i:section+i-1); 
    
    signal_hamm = signal_part.*hamming(length(signal_part));      
    
    signal_LPC = [lpc(signal_hamm,70) zeros(1,100)];                             
    
    DFT=fft(signal_LPC);
    N=length(DFT);                                
    df = (Fs/2)*linspace(0,1,N);  
    
    signal_PSD=(1/(N))*(abs(DFT)).^2; 
    
    plot(df(1:length(df)/2),signal_PSD(1:length(signal_PSD)/2));pause;
end

I should get frequencies around:

F1=654,44 Hz,

F2=911,56 Hz,

F3=2576 Hz,

F4=3921 Hz.

But I have mixed F1 with F2 and getting additional formants in higher frequencies.

or

Is it only guessing numbers of LPC coefficients with adjusting resolutoin of FFT ?

[Christopher McCausland]

Hi Piesiolaw,

With FFT's in MATLAB it is important to rescale the plot x-axis, the FFT appears to be plotted with respect to the number of samples (bins) used. While the shape will stay the same the rescaled axis should make much more sense,

Christopher

[Piesioslaw]

Thank you for answering,

I was experimenting with adding zeros to LPC but no progress there.

So how can i adjust/rescale the x-axis (frequencies) to final plot values?

Answer 1

[Christopher McCausland]

I Piesioslaw,

To recalculate your axis have a look at my example below;

% Compute a 1024 point fft for my signal
ecgcleanNoDCFFT = fft(ecgcleanNoDC, 1024);
% Convert axis to Hz 
Fs = 200;
f = (0:511)*((Fs/2)/512);
% Plot single side of spectrum
plot(f, abs(ecgcleanNoDCFFT(1:512)));
xlabel ('Frequency(Hz)'); % labels x-axis of spectral plot
ylabel ('Magnitude of Fourier Transform'); % labels y-axis of spectral plot

Here I have a sampling frequency of 200 Hz and a 1024 bin fft. I compute 'f' to calculate my bin frequencies which depends on the number of bins I want and the sampling frequency. As I am only plotting the single sided spectum I will only use half of the bins (as the other half give the same information but - x-axis). I can then use 'plot' to plot my newly calculated x-axis values 'f' against the absoloute value of my fft.

Let me know how you get on,

Christopher

[Piesioslaw]

Now it look like this

I still don't get it how it should be done. What is the best way to find these formants ?

Changed code :

    DFT=fft(signal_LPC,1024);
    f = (0:511)*((Fs/2)/512);
    N=length(DFT);                                
%     df = (Fs/2)*linspace(0,1,N);  
    
    signal_PSD=(abs(DFT)).^2; 
    
    plot(f,signal_PSD(1:length(signal_PSD)/2));pause;
end

Addinionaly i noticed that the pre-emphasis filter was set wrong, now i'm filtering higher freq to get lower. Still not what is should be.

signal_preemf = filter([1 -0.95],1,samples);