百闻不如一见

下面的数据来自于引力波开放科学中心（http://losc.ligo.org，新网址是https://www.gw-openscience.org）。这是GW170608引力波事件在Hanford观测站得到的结果（减去噪音之前，before noise substraction）。每个文件里包含了32秒钟的引力波测量数据，每秒钟16384个数据点，总计524288个数据（网站还有全套的数据，太多了。）。应变数据是无量纲量$\delta L /L$，数值非常小，乘以$10^{-20}$是为了便于处理。

数据处理采用了自己编写的SciLab程序（见本文的最后部分）。
32秒钟的数据如下，可以看到，应变数据在$\pm 1.2\times 10^{-18}$之间振荡。随后的几张图是同样的数据，但是观看的窗口变小了（注意横坐标）。

原始数据和用来进行光滑处理的初始值。以下四张图是同样的数据，差别只在于观察窗口的大小。

原始数据和经过插值处理后的光滑背景。以下四张图是同样的数据，差别只在于观察窗口的大小。

减去光滑背景后的结果（蓝色），以及经过32点平均后的结果（红色）。以下四张图是同样的数据，差别只在于观察窗口的大小。注意，即使经过了这些处理，数据的典型大小仍然是$10^{-20}$的量级，比引力波信号还要大1-2个数量级。

下面这几张图是为了说明如何求得引力波数据的包络线。关键在于求得局部的极大值。根据观察可以发现，极大值的周期大约是15个点，所以，我选取11个相连的点（红色点），用最小二乘法求出其抛物线近似（蓝色点）。可以看出来，局部极大值有明显的起伏，这与原始数据中的拍频现象是一致的。以下四张图是同样的数据，差别只在于观察窗口的大小。注意，在这种处理以后，包络线的起伏变小了，不考虑拍频的话，典型值小于$10^{-20}$了，但仍然比引力波信号大1个数量级左右。

下方包络线可以用类似方法求得。拍频也可以用类似方法去除（把上方包络线再求一次包络）。

这张图单独给出了极大值的位置。

GWtest.sce

//双斜杠表示注释行，其后的文字不做编译。
//注释
//每个文件里包含了32秒钟的引力波测量数据，每秒钟16384个数据点，总计524288个数据。
//引力波数据来自于 http://losc.ligo.org
//新网址是 https://www.gw-openscience.org/about/ 

clear ;//清除此前的所有变量
DataFileName="d:\SciLab6\GWtest\H-H1_LOSC_C01_16_V1-1180922478-32.txt"; //文件名
//这是GW170608引力波事件在Hanford观测站得到的结果（减去噪音之前，before noise substraction）。
Data2d=fscanfMat(DataFileName);//读文件
y1=Data2d'*10**20;//应变数据（strain data）保存在y1里，单位是10^-20
//应变数据是无量纲量， delta L /L，数值非常小，乘以10^20是为了便于处理
x1=[1:524288]/16384;//横坐标是时间，每秒钟16384个数据

figure;//定义一幅新图
plot(x1,y1);//画出引力波的数据
xlabel("时间（秒）", "fontsize", 4);//标出x轴
ylabel("应变（10^-20）", "fontsize", 4, "color", "blue");//标出y轴
title("引力波数据","fontsize", 6);//给出标题

//得到缓慢变化的背景
x2=([1:4096]*128-64)/16384;//每128个点为一组
[y2, d2] = lsq_splin(x1, y1, x2);//最小二乘法的样条函数

figure;//定义一幅新图
plot(x1,y1);//画出引力波的数据
plot(x2,y2,"rd");//光滑背景的初选点
//这两组数据出现在同一张图上
xlabel("时间（秒）", "fontsize", 4);//标出x轴
ylabel("应变（10^-20）", "fontsize", 4, "color", "blue");//标出y轴
title("引力波数据","fontsize", 6);//给出标题
legend(["引力波原始数据";"光滑背景初值"]);//标出每条曲线

yy1=interp1(x2,y2,x1,'spline');//对光滑背景进行插值，得到与应变数据一样多的结果
//光滑背景只有应变数据的很小一部分，所以要插值
figure;//定义一幅新图
plot(x1,y1,"bo-");//画出引力波的数据
plot(x1,yy1,"rd-");//插值后的光滑背景
//这两组数据出现在同一张图上
xlabel("时间（秒）", "fontsize", 4);//标出x轴
ylabel("应变（10^-20）", "fontsize", 4, "color", "blue");//标出y轴
title("引力波数据","fontsize", 6);//给出标题
legend(["原始数据";"光滑背景"]);//标出每条曲线

//对得到的结果再进行平均，每32个点求平均。
s1=y1-yy1;

for i=1:16384
    ii=32*(i-1)+1;
    iii=ii+31;
    rx(i)=(32*(i-1)+16)/16384;
    ry(i)=mean(s1(ii:iii));
end

figure;
plot(x1,s1);//应变数据减去插值后的光滑背景
plot(rx,ry,"ro-");//32点平均后的数据
xlabel("时间（秒）", "fontsize", 4);//标出x轴
ylabel("应变（10^-20）", "fontsize", 4, "color", "blue");//标出y轴
title("减去背景的引力波数据","fontsize", 6);//给出标题
legend(["减去背景后的数据";"32点平均值"]);//标出每条曲线

//保存数据
DataFileName2="d:\SciLab6\GWtest\savetest.txt";
rr=[rx ry];
fprintfMat(DataFileName2, rr, "%5.9f");

figure;
//寻找局部极大值

AA=[11,0,110;0,110,0;110,0,1958];
//假定11个初值的x坐标为[-5:5]，就得到二次型拟合的A矩阵如上
//参见王新民 董小刚《计算方法简明教程》2.5最佳平方近似，第72-75页
//或者任何一本计算方法的教材

ya=s1(40001:45000);//随便选了5000个数据
nmax=length(ya);
xa=[1:nmax]/16384;
plot(xa,ya);//画出这些数据
nn=1;// 峰的位置：初始猜测值

ii=0;
while (nn+30)<nmax
    ii=ii+1;
    nn=nn+15;
    za=ya((nn-10):(nn+10));//以初始猜测值为中心，选取21个数据
    [yt xt]=max(za); // 在上述21个数据里，yt是最大值，xt是最大值的位置
    xt1=xa((nn+xt-10-5):(nn+xt-10+5));//以xt为中心，选取11个位置
    yt1=ya((nn+xt-10-5):(nn+xt-10+5));//以及相应的数据
//    yt1=yt1';
    plot(xt1,yt1,"rd-");//画出待拟合的11个数据
    uu0=sum(yt1);
    uu1=sum([-5:5].*yt1);
    uu2=sum([-5:5].*[-5:5].*yt1);
    BB=[uu0;uu1;uu2];
    stemp=[-5:5];
    xc=linsolve(AA,-BB);//以上用最小二乘法进行二次函数的拟合
    yt2=xc(3)*stemp.*stemp+xc(2)*stemp+xc(1);//给出拟合值
    plot((stemp+nn+xt-10)/16384,yt2,"bd-");//画出拟合数据
    [ymax xmax]=max(yt2);
    nn=nn+xt-10+xmax-6;
    rrx(ii)=nn/16384;
    rry(ii)=ymax;
end
plot(rrx,rry,"bo-");

figure;
plot(rrx,rry,"bo-");