%% Data processing for v(t),i(t) measurements of mains power 
%% consumption by computing equipment.

% read the raw data file, in the current directory:
eval('raw_data_vi');

% list of measurements to be processed
list = { ...	
	'backup_standby'
	'backup_3spundown'
	'backup_normal'
	'backup_fullcpu'
	'backup_fullcpu_muchdisk'
	'solar_off'
	'solar_onguiidle1disk'
	'solar_onguiidle4disk'
	'solar_onguiidle1diskthrash'
	'solar_onguifullcpu4disk'
	'p3fn_off'
	'p3fn_idle'
	'p3fn_fullcpu'
	'p3fn_full3raid'   
        'penguin_idle3pow'
        'penguin_idle1pow'
        'penguin_diskload'
        'penguin_fullcpu'
        'penguin_fullcpudiskload'
	'big_crtmon_off'
	'big_crtmon_onblank'
	'big_crtmon_ongui'
	'tftmon_off'
	'tftmon_ongui'
};


% groups to plot on same figure
groups = { [1:5], [6:10], [11:14], [15:19], [20:22], [23:24] };

N = numel(list); % number of measurements
Vdiv = 0.018036; % voltage divider ratio
Idiv = 10; % amps per out-volt of current sensor [NOT TRUE for item `backup': Idiv=1]
Nc = 2; % cycles
Nt = 1e4;  % timepoints
Nh = 32;  % max harmonic order
V    = NaN*ones(Nt,N);
I    = NaN*ones(Nt,N);
Vft  = NaN*ones(Nh+1,N);
Ift  = NaN*ones(Nh+1,N);
Vpk  = NaN*ones(1,N);
Vrms = NaN*ones(1,N);
Ipk  = NaN*ones(1,N);
Irms = NaN*ones(1,N);
P    = NaN*ones(1,N);

disp(' ')
disp('Name                             Vrms   Vpk     Irms   Ipk      Srms   Preal   PF')
for n=1:N,
	% read data
	eval(['v = ', list{n}, '(:,1);']);
	eval(['i = ', list{n}, '(:,2);']);
	v = v/Vdiv;
	% ugly, because of `backup' being measured with series 1ohm instead
	% of with a 100 mV / A current probe
	if ~strcmp( list{n}(1:7), 'backup_' ),
		i = 10*i;
	end
	
	% perhaps filter it:
	dofilt = 1;
	if dofilt,	
		fn = 8;
		V(:,n) = filtfilt( ones(1,fn)/fn, 1, v );
		I(:,n) = filtfilt( ones(1,fn)/fn, 1, i );
	else
		V(:,n) = v;
		I(:,n) = i;
	end
	% make harmonics 0:Nh of frequency-domain view
	Vft_tmp    = fft(V(:,n))*2/length(V(:,n));
	Vft_tmp(1) = Vft_tmp(1)/2;
	Vft(:,n)   = Vft_tmp(1:2:(2*Nh+1));
	Ift_tmp    = fft(I(:,n))*2/length(I(:,n));
	Ift_tmp(1) = Ift_tmp(1)/2;
	Ift(:,n)   = Ift_tmp(1:2:(2*Nh+1));
	% remove '_' from legend text
	legends{n} = strrep(list{n},'_',' ');
	% calculate 1-number indices (peak, RMS, power)
	Vrms(n)    = sqrt(sum(V(:,n).^2)/size(V,1));
	Vpk(n)     = max(abs(V(:,n)));
	Irms(n)    = sqrt(sum(I(:,n).^2)/size(I,1));
	Ipk(n)     = max(abs(I(:,n)));
	S(n)       = Vrms(n)*Irms(n);
	P(n)       = sum(V(:,n).*I(:,n))/size(V,1);
	pf(n)      = P(n)/S(n);
	
	% display text data
	fprintf('%-30s   %5.1f  %5.1f   %5.3f  %5.3f   %5.1f  %5.1f   %5.2f\n', ...
		['''', list{n}, ''''], Vrms(n), Vpk(n), Irms(n), Ipk(n), ...
		 S(n), P(n), pf(n) );
end

% timepoints
t = linspace(0,40e-3,Nt);

% plot some example voltage curves, along 
% with a pure sinusoid of the same peak
vplots = [1, 10, 11];
figure;
subplot(2,1,1);
plot(t,V(:,vplots));
hold on;
plot(t, max(max(abs(V(:,vplots))))*sin(2*pi*50*t), 'k--');
hold off;
title('voltage time-curves');
xlabel('t [s]');
ylabel('V(t) [V]');
subplot(2,1,2);
bar(abs(Vft([1,3:9],:))./(ones(8,1)*abs(Vft(2,:))));
set(gca, 'XTickLabel', [0,2:Nh]);
xlabel('harmonic order, starting with 0 (d.c.)');
ylabel('harmonic ratio (value/fundamental)');

% plot
for g=1:numel(groups),
	figure;
	subplot(2,1,1);
	plot(t, I(:,groups{g}));
	hold on;
	plot(t, max(max(abs(I(:,groups{g}))))*sin(2*pi*50*t), ':k');
	xlabel('time [s]');
	ylabel('current [A]');
	legend(legends(groups{g}));
	grid
	subplot(2,1,2);
	bar(abs(Ift(1:8,groups{g})));
	xlabel('harmonic order');
	set(gca, 'XTickLabel', 0:Nh);
	ylabel('harmonic magnitude [A]');
	legend(legends(groups{g}));
	set(gca,'YMinorGrid','on'); 
end