“须菩提！如恒河中所有沙数，如是沙等恒河，于意云何？是诸恒河沙宁为多不？”须菩提言：“甚多，世尊！但诸恒河尚多无数，何况其沙。”

现在是大数据和人工智能的时代，我们必须跟上时代的步伐，为处理大数据做些准备。学好大学普通物理学就是很好的开始，但是仅仅满足于传统课程里的各种解析解和微扰方法，恐怕是不够的。学一些计算方法，能够自己编写一些程序、实现一些算法，就是又前进了一步。同时你会发现，物理学课程里强调的建立物理图像、抓住主要矛盾，是解决任何实际问题的法宝，可以帮助你判断大数据的收集、建模和分析过程中的优势和不足。

很多领域里都有显而易见的相关项目。比如说，天文学里的“巡天计划”，物理学里的“高能粒子探测”，生物学里的基因组、蛋白质甚至大脑等各种计划，材料科学里也有“材料基因组计划”，更别说天气预报以及各种社会数据的实时监控了。但是，这些项目都不太适合初学者，原因各种各样，比如说，不一定有解（搜索外星文明的“SETI”计划），拿不到数据（实时交易信息，如股市、期货甚至淘宝和京东的销售数据），硬件跟不上（比特币，高频交易），没有合适的模型（人工智能的下棋程序，现在还仅仅是知其然而不知其所以然），等等。我觉得有一个比较适合初学者的项目，就是分析引力波的数据。

为什么说，引力波数据可能适合初学者学习计算方法呢？大致有这么几个原因：

这些数据是公开的，而且有科学意义，不像是传统课本中的习题，更多的只有教学意义。

这些数据比较大，但又是个人能够处理的。这些数据比任何中学课程中碰到的数据多得多，也比目前我们大多数大学基础课程中碰到的数据量大。

这些数据处理有标准答案（也就是过去几年里得到的6次引力波事件），但是一般人并不知道如何得到这个标准答案（我就不知道）。

涉及到的数据处理方法很多。因为引力波信号很弱，探测数据中包含了大量的噪音，如何消除这些噪音就是非常困难的问题。

这些数据处理有标准的教材（网络上的教学说明），也就是LIGO团队设立的引力波开放科学中心里的解说，采用的也是开放源代码的程序（比如Python）。每个人只要愿意，都可以检验其算法正确与否，进而可以自己实现其算法。比如说，可以试一试用SciLab实现它。在此过程中，自然就会熟悉各种常用的计算方法。

我对引力波探测是个外行，但是我很关心这件事，写过不少博文（引力波博文汇总)，还用GW170817引力波事件介绍过引力波数据分析的困难。最近，我写了个简单的SciLab程序，可以展示引力波数据，并做一些最基本的处理。

利用SciLab，可以从这些数据里了解很多东西。

比如说，用abs（取绝对值）和max（取最大值）很容易知道，原始数据的最大值是$1.18\times 10^{-18}$。类似地可以得到最小值是$1.24\times 10^{-25}$。数据横跨了7个数量级。差分数据给出的最小值也在这个数量级：一阶差分的最小值是$1.86\times 10^{-25}$，一阶差分的最小值是$1.46\times 10^{-25}$，似乎意味着能够测量的最小值大约是$2\times 10^{-26}$，对应着8个数量级的跨度。再考虑到最大值应该还小于探测器的上限（我记得有套数据的最大值是$10^{-17}$的量级），所以总跨度可能在9-10个数量级左右。探测器应该是高灵敏的CCD阵列，每个阵列有大约100万个像素，每个像素有12个有效灰阶，这还是有可能的。

再比如说，在0.1秒的时间里（引力波事件的典型持续时间），数据的最大值和最小值的差别仍然是$10^{-18}$的量级，而且相邻两个数据的差的最大值也达到了$10^{-19}$的量级（min(abs(diff(y)))），最小值则是$10^{-23}$的量级。要知道，在这个时间段里，只有1300多个数据，从中提取引力波的信号，自然就是非常艰巨的任务了。

而这还仅仅是开始。

采用一些简单的方法，可以去除低频的扰动（如下图所示，蓝色），但是离提取引力波信号还远得很呢。也有一些方法可以去除高频扰动（32点平均，红色）。

下面这个SciLab程序可以进行简单的读写文件，确定光滑的背景（低频扰动），确定高频峰的位置和幅度（确定局部极大值）。但是距离从原始数据得到引力波信号还远得很。引力波科学开放中心（ https://www.gw-openscience.org）有数据、有教程，感兴趣的读者可以去看看。当然，最好是自己可以编个程序来提取引力波信号。

更多的数据图可以参看我的另一篇博文计算方法之引力波数据的初步显示。

初步掌握了SciLab以后，就可以用它来实现各种常见的算法，对学习计算方法和大数据处理技术做些准备。

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-");