clear all
close all
clc

[y, fs] = audioread("sound/modulator22.wav");
ranges = [17000, 20000; 29000, 37000; 41000, 46000];
one = y(ranges(1,1):ranges(1,2));
two = y(ranges(2,1):ranges(2,2));
three = y(ranges(3,1):ranges(3,2));

word = one;

n = length(word);   
f = (0:n-1)*(fs/n);
f1 = 0;%Hz
f2 = 4000;%Hz
idx = find(f >= f1 & f <= f2); %define the index of the freq range
f = f(idx);
y = fft(word, n);% compute DFT of input signal
power = abs(y).^2/n;
power = power(idx);
[val, ind] = max(power);

%lowpass for the formant

Fc = 2000; % define the cutoff frequency of the low-pass filter
[b, a] = butter(6, Fc/(fs/2), 'low'); % design a 4th-order Butterworth low-pass filter
Pxx_filt = filter(b, a, power); % apply the filter to the power spectrum
length(Pxx_filt)
length(f)


figure;

subplot(1,2,1) % time plot 
plot(0:1/fs:(length(word)-1)/fs,word);
xlabel('Time (s)');
ylabel('Amplitude (a.u.)');

subplot(1,2,2) % freq range plot 
plot(f,10*log10(power/power(ind))); hold on;
plot(f, 10*log10(Pxx_filt), 'r');
xlabel('Frequency (Hz)')
ylabel('Power (dB)')