%Exemplo para uso atraves de script simples no Matlab
%Depois de produzir os vetores dos sinais analogicos, filtrados,
%condicionados e digitalizados, execute o codigo a seguir para emular o
%funcionamento de um IED

%Define as saidas do bloco no início de cada loop
tIed = 0;
PMed = 0;
PhasorV = complex(0,0);
PhasorI = complex(0,0);
PhasorS = complex(0,0);
DigOut = 0;

fprintf(1, 'IED inicializando\r\n');
N = 16;     %amostras por ciclo
Fa = N*60;  %frequencia de amostragem multipla da fundamental 60 [Hz]

%Coeficientes da TDF para a componente fundamental
FCos1 = zeros(1,N); %inicializa vetores de coeficientes com zero
FSen1 = zeros(1,N);
for i=1:N
    FCos1(i) = (sqrt(2)/N)*cos(2*pi*i/N); %Coefs. para valores eficazes
    FSen1(i) = (sqrt(2)/N)*sin(2*pi*i/N);
end;

%Coeficientes da TDF para a componente da terceira harmônica
FCos3 = zeros(1,N); %inicializa vetores de coeficientes com zero
FSen3 = zeros(1,N);
for i=1:N
    FCos3(i) = (sqrt(2)/N)*cos(2*pi*3*i/N);
    FSen3(i) = (sqrt(2)/N)*sin(2*pi*3*i/N);
end;

%Inicializa buffers de amostras de tensao e corrente para filtros da TDF
LastI = zeros(N,1);
LastV = zeros(N,1);
LastPInst = zeros(N,1);

%Zera contador de iteracoes e tempo absoluto
cont = 0;
tsim = 0;

%Sinaliza que o setup foi feito uma vez
fprintf(1, 'IED configurado\r\n');

%Aqui vamos iterar por todo o conjunto de dados de entrada
%Atenção:
%Todos os dados de entrada devem ser indezados com (iter).
%Todos os resultados que devem ser armazenados em cada instante de tempo
%devem ser feitos colocando (iter) em cada variavel de saida.
%Vejam a ocorrencia de (iter) no codigo a seguir
for iter=1:size(Ch1_VAN,1)
    %Converte entrada 1 de corrente para numero de engenharia
    %Os ganhos abaixo consideram:
    % - ADC de 12 bits, unipolar, 3,3V de fundo de escala
    % - Canal de corrente CH1 com fundo de escala até 8A de pico
    % - Canal de tensao CH2 com fundo de escala até 350V de pico
    IInput = double(Ch4_IA(iter)) * (3.3/4095);
    IInput = IInput - 1.65;
    IInput = (IInput * 8)/1.65;

    %Converte entrada 2 de tensao para numero de engenharia 
    VInput = double(Ch1_VAN(iter)) * (3.3/4095);
    VInput = VInput - 1.65;
    VInput = (VInput * 350)/1.65;

    %Computa o FIR da TDF para a fundamental e 3a harmonica - da corrente
    LastI(1) = IInput;  % Coloca a amostra recebida mais nova no buffer
    %Convolucao com coeficientes
    IReal1=0; IImag1=0;
    for i=1:N
        IReal1 = IReal1 + FCos1(i)*LastI(i);
        IImag1 = IImag1 + FSen1(i)*LastI(i);
    end
    %Convolucao com 3a. harmonica
    IReal3=0; IImag3=0;
    for i=1:N
        IReal3 = IReal3 + FCos3(i)*LastI(i);
        IImag3 = IImag3 + FSen3(i)*LastI(i);
    end
    %Scroll da fila das entradas
    for i=N:-1:2 %loop for em sentido reverso
       LastI(i)=LastI(i-1); 
    end

    %Computa o FIR da TDF para a fundamental e 3a harmonica - da tensao
    LastV(1) = VInput;  % Coloca a amostra recebida mais nova no buffer
    %Convolucao com coeficientes
    VReal1=0; VImag1=0;
    for i=1:N
        VReal1 = VReal1 + FCos1(i)*LastV(i);
        VImag1 = VImag1 + FSen1(i)*LastV(i);
    end
    %Convolucao com 3a. harmonica
    VReal3=0; VImag3=0;
    for i=1:N
        VReal3 = VReal3 + FCos3(i)*LastV(i);
        VImag3 = VImag3 + FSen3(i)*LastV(i);
    end
    %Scroll da fila das entradas
    for i=N:-1:2 %loop for em sentido reverso
       LastV(i)=LastV(i-1); 
    end

    %calculo de potencia media pela definicao
    PInst = VInput*IInput;
    %Scroll da fila das potencias instantaneas
    for i=N:-1:2 %loop for em sentido reverso
       LastPInst(i)=LastPInst(i-1); 
    end
    LastPInst(1) = PInst;
    Pmedio=0;
    for i=1:N
        Pmedio = Pmedio + LastPInst(i);
    end;
    Pmedio = Pmedio / N; %potencia media em 1 ciclo

    %Um exemplo de algoritmo automatico para ligar alguma coisa em uma saida
    %digital, se a potencia for superior a um dado valor
    if Pmedio > 250 
        DigOut = 1;
    end;

    %Designa as saidas do bloco nesse intervalo de execucao
    tIed(iter) = tsim;
    PhasorI(iter) = complex(IReal1, IImag1);
    PhasorV(iter) = complex(VReal1, VImag1);
    PhasorS(iter) = PhasorV(iter)*conj(PhasorI(iter)); %potencia complexa considerando apenas fasor da 1a harmonica
    PMed(iter) = Pmedio;

    if(mod(cont,120)==0) %a cada 120 amostras de 1/Fa --> imprime algo
        fprintf(1,'t(%1.3f) V = %3.1f /_ %3.1f°\t', tsim, abs(PhasorV(iter)), angle(PhasorV(iter))*180/pi - angle(PhasorV(iter))*180/pi);
        fprintf(1,'I = %3.1f /_ %3.1f°\t', abs(PhasorI(iter)), angle(PhasorI(iter))*180/pi - angle(PhasorV(iter))*180/pi);
        fprintf(1,'S = %3.1f /_ %3.1f°\n', abs(PhasorS(iter)), angle(PhasorS(iter))*180/pi);
        cont=0;
    end

    %Incrementa tempo de simulacao para proxima iteracao
    tsim = tsim + 1 / Fa;
    cont = cont + 1;
end

tsim = tIed;