Finding formants in speech signal

17 views (last 30 days)
Piesioslaw
Piesioslaw on 3 Feb 2022
Commented: Piesioslaw on 4 Feb 2022
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 ?
  2 Comments
Piesioslaw
Piesioslaw on 3 Feb 2022
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?

Sign in to comment.

Answers (1)

Christopher McCausland
Christopher McCausland on 4 Feb 2022
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
  1 Comment
Piesioslaw
Piesioslaw on 4 Feb 2022
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);

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by