clear all
close all
clc

[y, fs] = audioread("sound/modulator22.wav");
ranges = [17000, 21000; 30500, 36000; 41500, 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 = three;

n = length(word);   
f = (0:n-1)*(fs/n);
f1 = 0;%Hz
f2 = 2500;%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);%limit the signal to the frequency ROI 
[val, ind] = max(power);

%lowpass for the formant, moving average

for j = 1:length(idx)
  tot = 0;
  for k = j-18:j
    if k <=0
      tot = tot + power(1);
    else
      tot = tot + power(k);
    endif
  endfor
  power_avg(j) = 1/6*(tot);
endfor

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(power_avg/power(ind)), 'r');
xlabel('Frequency (Hz)')
ylabel('Power (dB)')