% Piecewise linear approximation (PLA)
% See: HJLM. Vullings, Automated ECG segmentation with Dynamic Time Warping
function [y1 PLA] = PLA(input,s,th)
% input: source data
% s: step
% th: threshold
if nargin<3
th=10;
end
if nargin<2
s=10;
end
n=length(input);
% Low pass filter
% [b,a] = cheby1(3,0.5,40/62.5);
% y1=filter(b,a,input);
% b = fir1(31,40/62.5);
% y=conv(b,input);
% y1=y(int16((length(y)-length(input))/2):(int16((length(y)-length(input))/2+length(input)-1)));
s1=s;
pp=1;
PLA(pp)=1;
pp=pp+1;
i=1;
while i<n
if (i+s1>=n)
i_plus_s=n;
else
i_plus_s=i+s1;
end
interrupt=0;
while (interrupt==0)
j=i+1;
while j<=i_plus_s
distance=input(i_plus_s)-input(i);
dcur=input(j)-input(i)-(distance*(j-i)/(i_plus_s-i));
% distance is larger than threshold, truncate old line
if abs(dcur)>th
s1=j-i;
i_plus_s=i+s1;
j=i+1;
interrupt=1;
continue;
end
j=j+1;
end
if interrupt==1
PLA(pp)=j-1;
pp=pp+1;
i=j-2;
s1=s;
else
% distance is smaller than threshold, expand old line
if (i_plus_s>=n)
i_plus_s=n;
break;
else
if ((i_plus_s+s1)>=n)
i_plus_s=n;
else
i_plus_s=i_plus_s+s1;
end
end
end
end
i=i+1;
end
PLA(pp)=n;
y1=input;
% PLA=y1(int16((length(y1)-length(input))/2):(int16((length(y1)-length(input))/2+length(input)-1)));