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

2017年05月10日 18:25:47

在查看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

gmx_smxvg.bsh

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

脚本支持两个选项:

-n指定平滑时所用的点数, 必须是奇数. 默认为-n 5, 表示使用数据点及其前后两点的平均值作为平滑值. 点数越大, 平滑后曲线越光滑.
-c指定平滑数据所在的列编号. 默认为-c 2, 表示平滑xvg文件的第2列数据.

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

下图是一个简单的例子