xvg曲线中心移动平均平滑脚本

类别:    标签: gmx bash   阅读次数:   版权: (CC) BY-NC-SA

在查看GROMACS分析程序生成的xvg文件时, 有时候曲线不够平滑, 上下波动很厉害. 为美观或其他目的有时候需要对曲线进行平滑. 常见的绘图软件一般都支持曲线平滑, 可选方法也很多. 其中最简单的是移动平均算法(moving average). 但简单的移动平均算法平滑出来的曲线存在滞后性, 对于有峰和谷的曲线不适合. 一个解决方法是采用基于中心点的移动平均算法. 简单地说, 每个点的平滑值是其自身及前后若干点数值的平均值. 下面这段脚本专门用于xvg文件的平滑, 供需要的人参考.

gmx_smxvg.bsh
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
usage="\
>>>>>>>>>>>>>>>>    gmx_smxvg      <<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>    Jicun Li       <<<<<<<<<<<<<<<<
>>>>>>>>>>>>    2017-05-10 14:44:21    <<<<<<<<<<<<
>> Usage: gmx_smxvg {-n 5} {-col 2} <File.xvg...>
>>    -n: Number of Points to Run Smoothing Average, Shoule be ODD
>>    -c: Column Number of the Data to Smooth"
[[ $# -lt 1 ]] && { echo "$usage"; exit; }

pnt=5;    col=2
opt=($@); N=${#opt[@]}
for((i=0; i<N; i++)); do
  arg=${opt[$i]}; j=$((i+1))
  [[ $arg =~ -n ]] && { pnt=${opt[$j]}; opt[$i]=""; opt[$j]=""; }
  [[ $arg =~ -c ]] && { col=${opt[$j]}; opt[$i]=""; opt[$j]=""; }
done

[[ $((pnt%2)) -eq 0 ]] && { echo 'ERROR!!! Number of Points should be ODD!'; exit; }

for file in ${opt[@]}; do
  awk -v pnt=$pnt -v col=$col ' \
  /^[#@]/  {print; next}
  /^[^#@]/ {
    n++; j=n%pnt
    if(n<pnt) {
      sma += $col
      if(n<(pnt+1)/2) {
        for(i=1; i<=col; i++) printf "%12.6f", $i
        printf "%12.6f", sma/n
        for(i=col+1; i<=NF; i++) printf "%12.6f", $i
        print ""
      }
    } else {
      sma += $col-f[j]
      now=(n-(pnt-1)/2)%pnt
      split(txt[now], dat)
      printf "%12.6f", dat[1]
      for(i=2; i<=col; i++) printf "%12.6f", dat[i]
      printf "%12.6f", sma/pnt
      for(i=col+1; i<=NF; i++) printf "%12.6f", dat[i]
      print ""
    }
    f[j]=$col; txt[j]=$0
  }
  ' $file >${file%.xvg}~sm.xvg
done

脚本支持两个选项:

原始数据及其平滑后的数据会另写入原文件名~sm.xvg文件中. 其中第n列数据的平滑值会添加到n+1列中.

下图是一个简单的例子

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 随机抽样一致性算法matlab示例代码
后一篇: 四种计算自由能方法的原理示例教程

访问人次(2017年1月27日起): | 最后更新: 2017-10-25 23:46:46 UTC | 版权所有 © 2008 - 2017 Jerkwin