runline.m 836 B

123456789101112131415161718192021222324252627282930313233
  1. function y_line=runline(y,n,dn)
  2. % Running line fit (local linear regression)
  3. %
  4. % Usage: y_line=runline(y,n,dn);
  5. %
  6. % Inputs:
  7. % y: input 1-d time series (real)
  8. % n: length of running window in samples
  9. % dn: stepsize of window in samples
  10. %
  11. % Outputs:
  12. % y_line: local line fit to data
  13. y=y(:);
  14. nt=length(y);
  15. y_line=zeros(nt,1);
  16. norm=y_line;
  17. nwin=ceil((nt-n)/dn);
  18. yfit=zeros(nwin,n);
  19. xwt=((1:n)-n/2)/(n/2);
  20. wt=(1-abs(xwt).^3).^3;
  21. for j=1:nwin,
  22. tseg=y(dn*(j-1)+1:dn*(j-1)+n);
  23. y1=mean(tseg);
  24. y2=mean((1:n)'.*tseg)*2/(n+1);
  25. a=(y2-y1)*6/(n-1); b=y1-a*(n+1)/2;
  26. yfit(j,:)=(1:n)*a+b;
  27. y_line((j-1)*dn+(1:n))=y_line((j-1)*dn+(1:n))+(yfit(j,:).*wt)';
  28. norm((j-1)*dn+(1:n))=norm((j-1)*dn+(1:n))+wt';
  29. end
  30. mask=find(norm>0); y_line(mask)=y_line(mask)./norm(mask);
  31. indx=(nwin-1)*dn+n-1;
  32. npts=length(y)-indx+1;
  33. y_line(indx:end)=(n+1:n+npts)'*a+b;