%% guide_to_wia.m - scripts for web page
%  KHP  (28/12/12)

%% load data
P=load('..\data\Ao_Pressure.txt');
U=load('..\data\Ao_Velocity.txt');
E=load('..\data\Ao_ECG.txt');
SR=1000;                           % sampling rate

%% look at data
pp=(P-min(P))/(max(P)-min(P));
uu=U/max(U);
ee=E/max(E);
N=length(P);
t=(0:N-1)/SR;
% figure, grid on, hold on;
% plot(t,pp+1,t,uu,'r',t,ee-1,'g');
% axis tight;
% ylabel '        E                                 U                      P';
% xlabel 't (s)'
% title 'Normalised P, U and E'
% zoom xon;

%% power spectrum
fp=fft(P-mean(P));
fu=fft(U-mean(U));
f=(0:N-1)*SR/N;
% figure, grid on, hold on;
% plot(f(1:round(N/2)),abs(fu(1:round(N/2))),'r',f(1:round(N/2)),abs(fp(1:round(N/2))),'b');
% zoom xon;
% xlabel 'f (Hz)';
% ylabel 'power';
% title 'Power spectra of P and U';

%% savitsky-golay differences
Npoly=2;                 % Order of polynomial fit
F=51;                % Window length
[b,g]=sgolay(Npoly,F);   % Calculate S-G coefficients
HalfWin=((F+1)/2) -1;
p=P-min(P); u=U;       % retain P and U as the original data
dp=zeros(N,1); du=dp; ps=dp; us=dp;
for n=(F+1)/2:N-(F+1)/2
  % Zero-th derivative (smoothing only)
  ps(n)=dot(g(:,1),p(n-HalfWin:n+HalfWin));
  us(n)=dot(g(:,1),u(n-HalfWin:n+HalfWin));
  % 1st differential
  dp(n)=dot(g(:,2),p(n-HalfWin:n+HalfWin));
  du(n)=dot(g(:,2),u(n-HalfWin:n+HalfWin));
end
di=dp.*du;
% figure(30); clf; zoom xon;
% subplot(3,1,1);
% plot(t,ps),grid on, axis tight;
% ylabel 'P-Pd';
% title 'net wave intensity analysis'
% subplot(3,1,2);
% plot(t,us),grid on, axis tight;
% ylabel 'U';
% subplot(3,1,3);
% plot(t,di),grid on, axis tight;
% ylabel 'dI';
% xlabel 't (s)';

%% ensemble average
nb=Rwave(E);               % returns the indices of the R-wave peaks in the ECG signal
Nb=length(nb);             % the number of beats detected
nperiod=mean(diff(nb));    % find the average cardiac period
nint=round(1.2*nperiod);   % chose a beat interval slightly longer than nperiod
tens=(0:nint-1)/SR;
pens=zeros(Nb-1,nint);
for n=1:Nb-1
    pens(n,:)=ps(nb(n):nb(n)+nint-1);
    uens(n,:)=us(nb(n):nb(n)+nint-1);
end
pensavg=mean(pens);
penssd=std(pens);
uensavg=mean(uens);
uenssd=std(uens);
% figure, grid on, hold on;
% plot(tens,pens');
% plot(tens,pensavg,'k',tens,pensavg+penssd,'k:',tens,pensavg-penssd,'k:','LineWidth',2);
% axis tight;
% ylabel 'P-Pd';
% xlabel 't (s)';
% title 'ensemble average pressure'
% figure, grid on, hold on;
% plot(tens,uens');
% plot(tens,uensavg,'k',tens,uensavg+uenssd,'k:',tens,uensavg-uenssd,'k:','LineWidth',2);
% axis tight;
% ylabel 'U';
% xlabel 't (s)';
% title 'ensemble average velocity'

%% wave intensity of the ensemble averages
dpe=zeros(nint,1); due=dpe;
for n=(F+1)/2:nint-(F+1)/2
  % 1st differential
  dpe(n)=dot(g(:,2),pensavg(n-HalfWin:n+HalfWin));
  due(n)=dot(g(:,2),uensavg(n-HalfWin:n+HalfWin));
end
die=dpe.*due;
% figure(40); clf; zoom xon;
% subplot(3,1,1);
% plot(tens,pensavg),grid on, axis tight;
% ylabel 'Pe-Pd';
% title 'net wave intensity analysis'
% subplot(3,1,2);
% plot(tens,uensavg),grid on, axis tight;
% ylabel 'Ue';
% subplot(3,1,3);
% plot(tens,die),grid on, axis tight;
% ylabel 'dIe';
% xlabel 't (s)';

%% PU loop
figure(50); clf;
plot(uensavg,pensavg), grid on;






