## Copyright (C) 2023 Loic
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <https://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {} {@var{retval} =} freqRangeSpectrum (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn

## Author: Loic <Loic@LAPTOP-1F7SF0F7>
## Created: 2023-04-13

function [val, ind, f] = freqRangeSpectrum (signal, fs, f1, f2, formant, pad)
n = length(signal);  
if (pad)
	n = 2^nextpow2(n);
end 

f = (0:n-1)*(fs/n);

idx = find(f >= f1 & f <= f2); %define the index of the freq range
f = f(idx);

tic
y = fft(signal, n);% compute DFT of input signal
duration = toc

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/18*(tot);
endfor

figure;

subplot(1,2,1) % time plot 
if (pad)
	signal = [ signal, zeros( n-length(signal), 1)'];
end

plot(0:1/fs:(n-1)/fs,signal);
xlabel('Time (s)');
ylabel('Amplitude (a.u.)');

subplot(1,2,2) % freq range plot
plot(f,10*log10(power/power(ind)));
if formant == 1
  plot(f, 10*log10(power_avg/power(ind)), 'r');
endif
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
endfunction