clear all
close all
clc

%%%%%%%%%%%%%%%%%%
%function power = frequencySpectrum(signal, fs)
%
% Task: Display the power spectrum (lin and log scale) of a given signal
%
% Input:
%	- signal: the input signal to process
%	- fs: the sampling rate
%
% Output:
%	- power: the power spectrum
%
%
% Thomas Périn, thomas.perin@ecam.fr
% 20/04/2023
%%%%%%%%%%%%%%%%%%

signal = csvread('unknownsignal.csv');
fs = 300;                   %Sampling frequency
n = length(signal);
t = 0:1/fs:(n-1)/fs;
windowDuration = 1;

figure;
plot(t, signal);
title('Original Signal');
xlabel('time (s)');
ylabel('amplitude (a.u.)');

%%%% Windowing %%%%
rectangularWin = zeros(1, n);
for l_sample=1:windowDuration*fs
	rectangularWin(l_sample) = 1;
end

% plot rectangular window
figure;
plot(t, rectangularWin);
title('Rectangular Window');
xlabel('time (s)');
ylabel('amplitude (a.u.)');

% apply the rectangular window
for l_sample=1:n
	signal_rect(l_sample) = signal(l_sample) * rectangularWin(l_sample);
end

% plot rectangular signal
figure;
plot(t, signal_rect);
title('Signal with Rectangular Windowing');
xlabel('time (s)');
ylabel('amplitude (a.u.)');


%%%% Spectral Analisis %%%%
% compute DFT of input signal
y = fft(signal_rect, n);
% power of the DFT
power = abs(y).^2/n;

figure;
f = (0:n-1)*(fs/n);         % frequency range
plot(f,power);
title('Frequency Plot');
xlabel('frequency (Hz)');
ylabel('amplitude (a.u.)');

idx = find(power(61:81) == max(power(61:81)));
f(idx+60)