% Minimum Roughness with CET, % see http://climateaudit.org/2011/01/07/uc-on-mannian-smoothing/ % UC load c2010.dat % from http://hadobs.metoffice.com/hadcet/ , clear text & rename to c2010.dat c=c2010(:,end); X=NaN(352,252); Xr=X; for i=101:352 % series ending at i, 1759...2010 x=lowpassmin(c(1:i),1/50); % smooth from year c2010(1,1) to c2010(i,1), lowpassmin xr=lowpass(c(1:i),1/50,2,2); % smooth from year c2010(1,1) to c2010(i,1), minimum roughness X(1:i,i-100)=x; % store lowpassmin series Xr(1:i,i-100)=xr; % store mr series end close all figure plot(c2010(:,1),X) % plot all lpm smooths hold on kk=plot(c2010(:,1),X(:,end),'k'); % plot 1659 to 2010 smooth in bold set(kk,'LineWidth',3) title('lowpassmin') % result in http://www.climateaudit.info/data/uc/lpcet.png figure plot(c2010(:,1),Xr) % plot all mr smooths hold on kk=plot(c2010(:,1),Xr(:,end),'k'); % plot 1659 to 2010 smooth in bold set(kk,'LineWidth',3) title('minimum roughness') % result in http://www.climateaudit.info/data/uc/mr_cet.png figure cm=c2010(:,2:end-1); cmt=cm'; t=(1:36)-12; k=plot(1:12,cmt(:,82),'bo'); % coldest year set(k,'LineWidth',2) hold on k2=plot(1:12,cmt(:,348),'ro'); % warmest year set(k2,'LineWidth',2) k3=plot(1:12,cmt(:,end),'go'); % 2010 set(k3,'LineWidth',2) legend('1740','2006','2010') plot(t,[cmt;cmt;cmt],'k') % repeat, just to make figure look better hold on k=plot(1:12,cmt(:,82),'bo'); set(k,'LineWidth',2) k2=plot(1:12,cmt(:,348),'ro'); set(k2,'LineWidth',2) k3=plot(1:12,cmt(:,end),'go'); set(k3,'LineWidth',2) % result in http://www.climateaudit.info/data/uc/monthly_cet.png