clear all;
close all;
a=input("Inserisci un vettore nel formato [f A0 A1 ... An]")
if length(a) <3 
disp("Errore, devi inserire almeno 3 valori f, A0 e A1")
return
end
fa=a(1); %frequenza fondamnetale
Ta=1./fa; %Periodo fondamentale 
t_min=-50*Ta; %Consideriamo un numero sufficientemente elevato di periodi
t_max=50*Ta;
A0=a(2); %componente continua (o componente costante)
n=length(a); 
f_max=(n-2)*fa; %frequenza dell'armonica più elevata
Fs=1000*f_max; 
%frequenza di campionamento (simulazione di un segnale continuo)
Ts=1./Fs; %Periodo di campionamento
t=t_min:1/Fs:t_max; %vettore tempo - s
L=length(t);  %nuero di campioni 
f=Fs/L*(-L/2:L/2-1); %vettore frequenza - Hz
x=A0;
for indice=3:length(a)
    k=indice-2; % N.B. il primo elemento di a è la frequenza, il secondo elemento è la componente continua
  x=x+a(indice)*cos(2*pi*k*fa*t); 
end
figure(1)
subplot(2,1,1)
plot(t,x)   % segnale nel dominio del tempo
title("Segnale x(t)")
xlabel("t (s)")
ylabel("X(t)")
axis([min(t) max(t) min(x) max(x)]);
grid on

Y=fft(x,L)./L;  % trasformata FFT 

figure(1)
subplot(2,1,2)
plot(f,abs(fftshift(Y)),"LineWidth",2)
title("Spettro del segnale x(t)")
xlabel("f (Hz)")
ylabel("|fft(X)|")
Y_max=max(abs(Y));
axis([-2*f_max 2*f_max 0 Y_max]);
grid on;
P=trapz(t,x.^2)./(t_max-t_min);
sprintf("La potenza del segnale è: %d",P)
p_f=Y*Y';
sprintf("Y*Y' fornisce lo stesso valore di potenza: %d",p_f)
p=abs(Y).^2;
p=sum(p);
sprintf("sum(abs(Y).^2) fornisce lo stesso valore di potenza: %d",p)
fid=fopen("potenza.txt","w");
for indice=1:n
    fprintf(fid,"%d ",a(indice));
end
fprintf(fid,"%d ",P);
fclose(fid);