sdtrend remove annually cycle in matlab

function [y,psave,ssave]=sdtrend(x,T,n,plfg,descriptor)

%

%   [y,psave,ssave]=sdtrend(x,T,n,plfg,descriptor)

%

% x = input series

% T = period

% n = order of polynomials to fit

% plfg > 0 to plot series, residuals and seasonal trends

% descriptor = string for plots

if plfg

    �=figure;

    subplot(211);

    plot(x);

    title([descriptor,' with seasonal polys of degree : ',int2str(n)]);

    hold on;

end

colors=['y';'m';'c';'r';'g';'b'];

ncols=6;

colper=ncols;

if T < ncols

    colper=T;

end

nx=length(x);

resid0=zeros(size(x));

ssave=[];

psave=[];

%legstr=descriptor;

legstr=[];

cnum=1;

for j=1:T

    indexes=[j:T:nx];       % get the indexes starting at j and skipping by T

    xs=x(indexes);          % xs is x sampled at these indexes

    [pout,sout]=polyfit(indexes',xs,n);

    ssave=[ssave;sout];

    psave=[psave;pout];

    y=polyval(pout,indexes)'; % make it a column

    if plfg

        %figure(f1);

        hold on;

        plot(indexes',y,colors(cnum));

%         legstr=[legstr;sprintf('s = %0.2d',j)];

        cnum=rem(cnum,colper)+1;

    end

    resid0(indexes)=x(indexes)-y;

    

end

y=resid0;

if plfg

    % put on the legend for last plot

%     legend(legstr,2);   % 2 puts legend in upper left corner

    hold off;

    

    %figure;

    subplot(212);

    plot(resid0);

    title([descriptor,'resids from sdtrend with polys of degree : ',int2str(n)]);

end