function [S, H, C, V] = phasorlisa(SP, K)
% SUM OF PHASORS
%
%  [S, H, C, V] = phasor(SP, [K])
%
%  SP     vector containing starting points of each exponential.
%         From each starting point the length of a phasor and
%         the phase shift are computed. If no other arguments
%         are given, then indices of SP ..
%  K      If K is given, it can have two values:
%           -1      which uses indices of SP and inverts 
%           vector  each SP having a multiplier
%
%  S      sum signal
%  H      individual signals
%  C      cumulative sum of signals
%  V      visualization matrix
%
%
% Examples:
%   phasor([1 1 1]/3,-1);   % three-point MA
%                           % 0.33 [ e^(-jw) + e^(-j2w) + e^(-j3w) ]
%   phasor([2 4 -3], -[1 3 5]); 
%                           % 2 e^(-jw) + 4 e^(-j3w) - 3 e^(-j5w)
%
% A nice Java Applet:
% http://www.jhu.edu/~signals/phasorlecture2/indexphasorlect2.htm
%
% You can interrupt and change values:
% See also pauseButton.m in MathWorks File Exchange.
%
% TODO:
% See also arrow.m in MathWorks File Exchange.
%

% Jukka Parviainen, parvi@tkk.fi, 10-July-2007

if (nargin == 1)
  K = [1 : length(SP)];
elseif ((length(SP) > 1) & (K == -1))
  %disp('vastakkainen suunta');
  K = -[1 : length(SP)];
elseif ((length(SP) > 1) & (length(SP) ~= length(K)))
  error('Vector K must be either -1 or the same length as SP');
else
  %warning('** is everything ok?**');
end;

% INIT
res = 64;
per = 1;
pausmom = 0.2;

A = abs(SP);
T = angle(SP);

% H(w) = A(1) . exp(j K(1) w + T(1)) + 
%        A(2) . ...

% INIT: resolution and how many periods
%res = 64;
%per = 2;
w   = [0 : (2*pi)/res : per*(2*pi)];
L   = length(w);

%
A = vector2col(A); % see in the end of this file
T = vector2col(T);
K = vector2col(K); 

H = repmat(A, 1, L) .* exp(j * (K * w + repmat(T, 1, L)));
S = sum(H);
C = cumsum(H);

V = zeros(2*size(C,1),size(C,2));
V(2:2:end,   :) = C;
V(3:2:end-1, :) = C(1:end-1, :);
V = reshape(V, 2, size(V,1)/2, L);


SRmin = min(real(S));
SRmax = max(real(S));
SRDiff = SRmax-SRmin;
SRlimA = SRmin - 0.1*SRDiff;
SRlimB = SRmax + 0.1*SRDiff;

SImin = min(imag(S));
SImax = max(imag(S));
SIDiff = SImax-SImin;
SIlimA = SImin - 0.1*SIDiff;
SIlimB = SImax + 0.1*SIDiff;

SAmin = min(abs(S));
SAmax = max(abs(S));
SADiff = SAmax-SAmin;
SAlimA = SAmin - 0.1*SADiff;
SAlimB = SAmax + 0.1*SADiff;

STmin = min(angle(S));
STmax = max(angle(S));
STDiff = STmax-STmin;
STlimA = STmin - 0.1*STDiff;
STlimB = STmax + 0.1*STDiff;



if (nargout == 0)
  figure;
  subplot(2,4,[1 5]);
  p18 = plot(real([zeros(1,size(H,1)); H(:,1).']), imag([zeros(1,size(H,1)); H(:,1).']));
  axis([SRlimA SRlimB SIlimA SIlimB]);
  axis equal
  axis manual
  grid on;
  hold on;
  title({'Each phasor','in complex plane'});


  subplot(2,4,[2 6]);
  p1 = plot(real(V(:,:,1)), imag(V(:,:,1)));  
  axis([SRlimA SRlimB SIlimA SIlimB]);
  axis equal
  axis manual
  grid on;
  hold on;
  title({'Phasor sum','in complex plane'});
  ww = 1;  % täplää varten alustus
  
  subplot(2,4,3);
  p2 = plot(w(1), real(S(1)));
  axis([w(1) w(end) SRlimA SRlimB]);
  %axis equal
  axis manual
  title('real part -- x -- cos(\omega)');
  grid on;
  hold on;
  addlines_phasor([0.5 : 0.5 : 4], SRlimA, SRlimB);  
  
  subplot(2,4,7);
  p3 = plot(w(1), imag(S(1)));  
  axis([w(1) w(end) SIlimA SIlimB]);
  %axis equal
  axis manual
  title('complex part -- y -- sin(\omega)');
  grid on;
  hold on;
  xlabel('\omega');
  addlines_phasor([0.5 : 0.5 : 4], SIlimA, SIlimB);  
  
  % abs
  subplot(2,4,4);
  p4 = plot(w(1), abs(S(1)));
  axis([w(1) w(end) SAlimA SAlimB]);
  %axis equal
  axis manual
  title('absolute part -- abs(H(\omega))');
  grid on;
  hold on;
  addlines_phasor([0.5 : 0.5 : 4], SAlimA, SAlimB);  
    
  % angle
  subplot(2,4,8);
  p5 = plot(w(1), angle(S(1)));  
  axis([w(1) w(end) STlimA STlimB]);
  %axis equal
  axis manual
  title('angle part -- angle(H(\omega))');
  grid on;
  hold on;
  xlabel('\omega');
  addlines_phasor([0.5 : 0.5 : 4], STlimA, STlimB);  
  
  for n = 2 : L

    subplot(2,4,[1 5]);
    delete(p18);
    p18 = plot(real([zeros(1,size(H,1)); H(:,n).']), imag([zeros(1,size(H,1)); H(:,n).']));
    p17 = plot(real([H(:,n-1).'; H(:,n).']), imag([H(:,n-1).'; H(:,n).']),':');
    xlabel(['\omega = ' num2str(w(n))]);

    subplot(2,4,[2 6]);
    delete(p1);
    p1 = plot(real(V(:,:,n)), imag(V(:,:,n)));
%    p1 = arrow(real(V(:,:,n)), imag(V(:,:,n)));
    xlabel(['\omega = ' num2str(w(n))]);
    p6 = plot([real(S(n-1)) real(S(n))], [imag(S(n-1)) imag(S(n))],  'Color', [0 0 1]); 

    % lisätään täplä joka pi/2:n välein
    if (w(n) >= (ww * pi/2 - 0.001))
      if (w(n) < 2.2*pi)
        p32 = plot(real(S(n)), imag(S(n)), 'o', 'MarkerFaceColor', [0 0 0]);
        t32 = text(real(S(n)), imag(S(n)), ['\omega = ' num2str(w(n))]);
        ww  = ww+1;
      end;
    end;
    
    subplot(2,4,3);
%
%    p2 = plot([w(n-1) w(n)], [real(S(n-1)) real(S(n))]);    
%
    p2 = plot([w(n-1) w(n)], [real(S(n-1)) real(S(n))],'LineWidth',2);
    p8 = plot(repmat([w(n-1) w(n)]',1,size(H,1)), real([H(:,n-1).'; H(:,n).']));
    
    subplot(2,4,7);
%
%    p2 = plot([w(n-1) w(n)], [imag(S(n-1)) imag(S(n))]);    
%
    p2 = plot([w(n-1) w(n)], [imag(S(n-1)) imag(S(n))],'LineWidth',2);
    p8 = plot(repmat([w(n-1) w(n)]',1,size(H,1)), imag([H(:,n-1).'; H(:,n).']));    

    subplot(2,4,4);
%    p2 = plot([w(n-1) w(n)], [abs(S(n-1)) abs(S(n))]);    
    p2 = plot([w(n-1) w(n)], [abs(S(n-1)) abs(S(n))],'LineWidth',2);
    p8 = plot(repmat([w(n-1) w(n)]',1,size(H,1)), abs([H(:,n-1).'; H(:,n).']));        
    subplot(2,4,8);
%    p2 = plot([w(n-1) w(n)], [angle(S(n-1)) angle(S(n))]);    
    p2 = plot([w(n-1) w(n)], [angle(S(n-1)) angle(S(n))],'LineWidth',2);
    p8 = plot(repmat([w(n-1) w(n)]',1,size(H,1)), angle([H(:,n-1).'; H(:,n).']));            

    % INIT
    pause(pausmom);
  
    pauseButton;
  end;
  
end;

% axis auto  % return back



%% Functions vector2col/row.m
%
%function y = vector2col(x)
%% VECTOR2COL forces vector x to be a column vector y
%%
%%  y = vector2col(x)
%
%% Jukka Parviainen, parvi@tkk.fi, 10.7.2007
%
%if (prod(size(x)) ~= length(x))
%  error('x must be a vector, not a matrix');
%end;
%
%if (size(x,1) < size(x,2))
%  y = x.';
%else
%  y = x;
%end;
%
%
%function y = vector2row(x)
%% VECTOR2ROW forces vector x to be a row vector y
%%
%%  y = vector2row(x)
%
%% Jukka Parviainen, parvi@tkk.fi, 10.7.2007
%
%if (prod(size(x)) ~= length(x))
%  error('x must be a vector, not a matrix');
%end;
%
%if (size(x,1) > size(x,2))
%  y = x.';
%else
%  y = x;
%end;
%


function addlines_phasor(interV, limA, limB)
% a subfunction for addlines_phasor

for m = interV
  line(m*pi*[1 1], [limA limB],'Color',[0 0 0], 'LineStyle', '-.');
end;
line(pi*[1 1], [limA limB],'Color',[0 0 0]);
