似乎已经完成了phase lag的校正, Oh Yeah,

Reversivle Heat capacity  = best linear fitting slope 

跟另外一个近似方法实质完全一样;

应该用Reversible Heat capacity = best linear fitting slope * cos( phase lag)

这样才行.

我尝试了一下, 果然如此. 结果得到了大大改善.

其中找 最佳线性拟合的算法, (对 dQm/dt  和 -dTm/dt 两个向量, 也就是热流的调制部分, 以及(降温)速率的调制部分)

每次40个点, 逐点向后移动, 也就是40点的窗口宽度

每一次的计算:

40个x=dQm/dt, 40g个对应的y=-dTm/dt

1. 对两个40维向量, x(1:40), y(1:40), y对x做线性拟合, 得到: 斜率, 截距, 残差,相关系数

2. 变换两个向量为39维, x(1:39), y(2:40), y对x做线性拟合, 得到: 斜率, 截距, 残差,相关系数

3. 进一步变为38维....

...

30. 差不多, 变成11维的时候, 计算完毕就可以停下来了

找residue的平均值的模最小的点, 对应的slope是可逆热容的复数形式的模

phase lag = k*2*pi/N, k是移动的次数, 或slope向量的下标, N是调制的每个周期中数据点数.

然后slope(k) * cos(phaselag(k)) 就得到经过phase lag correction的可逆热容.

%function RevHeatCapacity;

close all

clear

clc

%

tic

N=120;

load MDSCdata

t=data(:,1).*60;

Tm=data(:,12)-smooth(data(:,12),N);

dTmdt=mmspder(t,Tm,t); % Modulated Heating rate

dQmdt=data(:,13)-smooth(data(:,13),N); % Modulation of Modulated Heat Flow

%  x and y for linear regression

% y= dQmdt;

% x=-dTmdt;

%

L=length(t);

n=N/2; % maximum dimension of y and x

for j=N:L-N

  for k=1:n-20

    x=-dTmdt(j:n-k+j);

    y=dQmdt(k+j-1:n+j-1);

    [p,s]=polyfit(x,y,1);

    Slope(k)=p(1);

    Intercept(k)=p(2);

    Residue(k)=s.normr;

    if 0

    if k>1 & Slope(k)>=0 & Residue(k)<0.3

        if Residue(k)<Residue(k-1)

        HeatCapacity(j)=Slope(k);

        break

        end

    end % find desired Heat capability

    end

  end % of k, to determine phase lag

index=find(Slope<0);

Residue(index)=inf;

if 1

index=find(Residue==min(Residue));

index=index(1);

HeatCapacity(j)=Slope(index);

end

phaselag(j)=index/N*2*pi;

end % to Loop advanced point by point along the vectors

%for kk=1:2

HC=smooth(HeatCapacity, 2*N);

phaselag=smooth(phaselag,N);

HC=smooth(HC.*cos(phaselag),N);

%end % smooth obtained HeatCapacity

figure, plot(HC), hold on, plot(data(:,4),'r')

legend('Reversible Heat Capacity by Matlab', 'Reversible Heat Capacity by Q100')

legend boxoff

figure, plot(Residue),hold on, plot(1:40, Residue,'ro')

['how long does it run: ',num2str(toc),' seconds']

%end % of function