function [power, duration] = frequencySpectrum(signal, fs, pad, verbose, rangefr)
%%%%%%%%%%%%%%%%%%
%function power = frequencySpectrum(signal, fs, pad)
%
% Task: Display the power spectrum (lin and log scale) of a given signal
%
% Input:
%	- signal: the input signal to process
%	- fs: the sampling rate
%	-pad: boolean if true, signal is padded with 0 to the next power of 2 -> FFT instead of DFT
%
% Output:
%	- power: the power spectrum
%
%
% Guillaume Gibert, guillaume.gibert@ecam.fr
% 25/04/2022
%%%%%%%%%%%%%%%%%%

n = length(signal);        % number of samples

if (pad)
	n = 2^nextpow2(n);
end

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

power = abs(y).^2/n;    % power of the DFT

[val, ind] = max(power); % find the mx value of DFT and its index

rangemin = rangefr(1);
rangemax = rangefr(2);

if (verbose)
  % plots
  figure;

  subplot(1,3,1) % time plot
  t=0:1/fs:(n-1)/fs; % time range
  %pad signal with zeros
  if (pad)
    signal = [ signal; zeros( n-length(signal), 1)];
  end
  plot(t, signal)
  xticks(0:0.1*fs:n*fs);
  xticklabels(0:0.1:n/fs);
  xlabel('Time (s)');
  ylabel('Amplitude (a.u.)');

  subplot(1,3,2) % linear frequency plot
  f = (0:n-1)*(fs/n);     % frequency range
  plot(f,power, 'b*'); hold on;
  plot(f,power, 'r');
  xlabel('Frequency (Hz)')
  ylabel('Power (a.u.)')
  xlim([rangemin rangemax]);


  subplot(1,3,3) % log frequency plot
  power_db = 10*log10(power/power(ind));
  plot(f, power_db);
  xlabel('Frequency (Hz)')
  ylabel('Power (dB)')
  xlim([rangemin rangemax]);

  subplot(1,3,3) % log frequency plot
  power_db = 10*log10(power/power(ind));
  plot(f, power_db);
  xlabel('Frequency (Hz)')
  ylabel('Power (dB)')
  xlim([rangemin rangemax]);

  % Shift power_db to all positive values
  shift_amount = min(power_db) * -1;
  power_db_shifted = power_db + shift_amount;

  % Find peaks and their indices
  [pks, locs] = findpeaks(power_db_shifted);

  % Interpolate to get a smooth line
  xi = linspace(min(f(locs)), max(f(locs)), 1000);
  yi = interp1(f(locs), pks, xi, 'pchip');

  % Shift the interpolated values back down and up by 10
  yi = yi - shift_amount + 10;

  % Smooth the line using a moving average filter
  windowSize = 15; % Adjust this value to change the amount of smoothing
  b = (1/windowSize)*ones(1,windowSize);
  yi = filter(b, 1, yi);
  yi = filter(b, 1, flip(yi)); % Apply the filter in reverse direction
  yi = flip(yi); % Flip the data back to original direction

  % Shift yi to all positive values
  shift_amount_yi = min(yi) * -1;
  yi_shifted = yi + shift_amount_yi;

  % Find peaks of the smoothed line
  [pks_smooth, locs_smooth] = findpeaks(yi_shifted);

  % Plot the envelope
  hold on;
  plot(xi, yi, 'r-'); % plot the envelope

  % Plot points and labels at the peaks
  for i = 1:length(pks_smooth)
    plot(xi(locs_smooth(i)), pks_smooth(i) - shift_amount_yi, 'ko');
    text(xi(locs_smooth(i)), pks_smooth(i) - shift_amount_yi, ['F' num2str(i)], 'VerticalAlignment', 'bottom');
  end

  hold off;

end