% The MATLAB program for p-q simulation algorithm
%
% -----------------------------------------------------------------------------------
%
% (C) 2017, Krzysztof Sozanski, Digital Signal Processing in Power Electronics Control 
% Circuits, second edition, Springer-Verlag, London, 2017. 
% 
% -----------------------------------------------------------------------------------
%
% 
  clear all; close all;
  NM=2^8; % number of samples per line period 
  fs=50*NM; % sampling frequency
  fM=50; % line frequency
  N_period=4; % number of periods
  N=N_period*NM; % number of samples
  fi=0; %phase shift
  t=(0:N-1)/fs; % time
  s1=sin(2*pi*fM*t+fi); c1=cos(2*pi*fM*t+fi);
  % ------- Orthogonal signals -------
  fi2=-5*pi/6;
  inA=sin(2*pi*fM*t+0+fi2); 
  inB=sin(2*pi*fM*t-2*pi/3+fi2); 
  inC=sin(2*pi*fM*t-4*pi/3+fi2);
  in1=sign(inA)/2; in2=sign(inB)/2; in3=sign(inC)/2; 
  iL3=in1-in2; iL2=in2-in3; iL1=in3-in1; 
  % ------- Clark transformation -------
  ialfa=(iL1 - 0.5*iL2 - 0.5*iL3)*(2/3)^0.5; 
  ibeta =( (3^0.5)/2*(iL2-iL3))*(2/3)^0.5;
  % ------- Park transformation -------
  p=c1.*ialfa-ibeta.*s1;
  q=s1.*ialfa+c1.*ibeta;
  % ----- Butterworth first-order high-pass filter
  fcr=10; % crossover frequency 10 Hz
  Fcr=fcr/fs*2; % relative crossover frequency 
  [b a]=butter(1,Fcr,'high');
  pf=filter(b,a,p); qf=filter(b,a,q);
  % ------- Park inverse transformation -------
  iCalfa = pf.*c1 + qf.*s1;
  iCbeta = -pf.*s1 + qf.*c1;
  % ------- Clark inverse transformation -------
  iC1 = (iCalfa)*(2/3)^0.5;
  iC2 = (-1/2*iCalfa + (3^0.5)/2*iCbeta)*(2/3)^0.5;
  iC3 = (-1/2*iCalfa - (3^0.5)/2*iCbeta)*(2/3)^0.5;
  % ------- Line current calculation -------
  iM1=iL1-iC1; iM2=iL2-iC2; iM3=iL3-iC3;
  