function [ z ] = potenza( )
%POTENZA Calcolo potenza in regime non sinusoidale
Fs=10000;
t=[0:1/Fs:1-1/Fs];
N=length(t);
F=50;
v=5*sin(2*pi*F*t)+2*rand*sin(2*pi*2*F*t+pi/(10*rand))+2*rand*sin(2*pi*3*F*t+pi/(10*rand))+
2*rand*sin(2*pi*4*F*t-pi/(10*rand))+2*rand*sin(2*pi*5*F*t+pi/(10*rand));
i=10*sin(2*pi*F*t+pi/(10*rand))+2*rand*sin(2*pi*2*F*t+pi/(10*rand))+
2*rand*sin(2*pi*3*F*t+pi/(10*rand))+2*rand*sin(2*pi*4*F*t-pi/(10*rand))+
2*rand*sin(2*pi*5*F*t+pi/(10*rand));
plot(t(1:1000),v(1:1000));
xlabel('t [s]');
ylabel('v [V]');
title('Tensione');
pause;
plot(t(1:1000),i(1:1000));
xlabel('t [s]');
ylabel('i [mA]');
title('Corrente');
pause;
Vn=fft(v);
In=fft(i);
f=(0:N-1)*Fs/N;
stem(f(1:500),abs(Vn(1:500))/N,'-b.');
xlabel('f [Hz]');
title('Spettro della tensione');
pause;
stem(f(1:500),abs(In(1:500))/N,'-b.');
xlabel('f [Hz]');
title('Spettro della corrente');
pause;
close;
%Potenza attiva
Pint=sum(v.*i)/N;
%Potenza attiva
P=sum(abs(Vn).*abs(In).*cos(abs(angle(Vn)-angle(In))))/N2;
%Potenza reattiva secondo Budeanu
Q=sum(abs(Vn).*abs(In).*sin(abs(angle(Vn)-angle(In))))/N2;
%Potenza apparente
S=sqrt(sum(abs(Vn).2)*sum(abs(In).2))/N2;
St=sqrt(P2+Q2);
%Potenza di distorsione
D=sqrt(S2-St2);
%Scomposizione
j=find(f==F);
V1=Vn([j N-j+2]);
I1=In([j N-j+2]);
Vh=Vn([1:j-1 j+1:N-j+1 N-j+3:N]);
Ih=In([1:j-1 j+1:N-j+1 N-j+3:N]);
%Potenza attiva fondamentale
P1=sum(abs(V1).*abs(I1).*cos(abs(angle(V1)-angle(I1))))/N2;
%Potenza reattiva fondamentale
Q1=sum(abs(V1).*abs(I1).*sin(abs(angle(V1)-angle(I1))))/N2;
%Potenza apparente fondamentale
S1=sqrt(sum(abs(V1).2)*sum(abs(I1).2))/N2;
S1t=sqrt(P12+Q12);
%Potenza armonica totale
Ph=sum(abs(Vh).*abs(Ih).*cos(abs(angle(Vh)-angle(Ih))))/N2;
%Potenza armonica apparente
Sh=sqrt(sum(abs(Vh).2)*sum(abs(Ih).2))/N2;
%Potenza armonica totale non attiva
Nh=sqrt(Sh2-Ph2);
%Prodotti in croce
V1Ih=sqrt(sum(abs(V1).2)*sum(abs(Ih).2))/N2;
VhI1=sqrt(sum(abs(I1).2)*sum(abs(Vh).2))/N2;
%Potenza apparente non fondamentale
Sn=sqrt(V1Ih2+VhI12+Sh2);
%Potenza apparente
S2=sqrt(S12+Sn2);
stringa=['Potenza attiva con il metodo dell''integrazione discreta: Pa = ' num2str(Pint) ' mW'];
disp(stringa);
stringa=['Potenza attiva con il metodo della FFT: Pa = ' num2str(P) ' mW'];
disp(stringa);
stringa=['Potenza reattiva secondo Budeanu: Qb = ' num2str(Q) ' mW'];
disp(stringa);
stringa=['Potenza apparente: S = ' num2str(S) ' mW'];
disp(stringa);
stringa=['Potenza apparente ottenuta dal triangolo delle potenze: St = ' num2str(St) ' mW'];
disp(stringa);
stringa=['Potenza di distorsione secondo Budeanu: D = ' num2str(D) ' mW'];
disp(stringa);
disp(' ');
stringa=['Potenza attiva fondamentale: P1 = ' num2str(P1) ' mW'];
disp(stringa);
stringa=['Potenza reattiva fondamentale: Q1 = ' num2str(Q1) ' mW'];
disp(stringa);
stringa=['Potenza apparente fondamentale: S1 = ' num2str(S1) ' mW'];
disp(stringa);
stringa=['Potenza apparente fondamentale ottenuta con il triangolo
delle potenze: S1 = ' num2str(S1t)
' mW'];
disp(stringa);
stringa=['Potenza armonica totale: Ph = ' num2str(Ph) ' mW'];
disp(stringa);
stringa=['Potenza armonica apparente: Sh = ' num2str(Sh) ' mW'];
disp(stringa);
stringa=['Potenza armonica totale non attiva: Nh = ' num2str(Nh) ' mW'];
disp(stringa);
stringa=['Potenza apparente non fondamentale: Sn = ' num2str(Sn) ' mW'];
disp(stringa);
stringa=['Potenza apparente: S = ' num2str(S2) ' mW'];
disp(stringa);