ncl的butterworth滤波

小懒C

    最近被滤波纠缠，因为发现ncl自带的butterworth滤波最后结果会趋向于0，汇报的时候被老师问这个数据最后怎么会没波呢？？于是就搞了一下butterworth滤波，因为butterworth滤波官方解释是不丢失数据的，但趋向于0这个完全就是表示后面数据可能不可用。程序在最后
  我用到的资料：1.http://bbs.06climate.com/forum.php?mod=viewthread&tid=58216&highlight=butter的文献以及《气候变率诊断与预测方法》
                       2.李建平教授的Fortran文件（程序基于上面所提到的文献）                       其他reference：http://www.seismosoc.org/Publications/BSSA_html/bssa_96-2/05055-esupp/（ncl基于此）
  以下以《气候变率诊断与预测方法》p208的例子进行讨论：

    对于ncl滤波尾部数据趋于0，ncl-talk有人给出了解决方法：I talked with one of the professors in my department who knows quite a bit about numerical methods, and he suggested padding the end of the series with 1.5*(order of filter/fca) zeros.   https://mailman.ucar.edu/pipermail/ncl-talk/2015-June/003008.html 按照这个我试了一下，尽管尾部数据有改善，但产生了大相位偏差（图1最左侧，应该0的时候振幅最大），同时末端也存在相位误差（最末端按照书上应该看得见一个波谷）

   于是我试了试，在数据最左侧增加1.5*(order of filter/fca) 个最开始的数据，后方加上1.5*(order of filter/fca) 个最末端数据，使用2阶的BW程序，于此同时我也将李建平的程序翻译成了ncl语言，结果如图2，两个数据之间存在较小的振幅与相位差，需要注意的是李建平老师的程序是没有数据丢失的，而ncl画图时没要padding进去的数据，如果采用ncl的其他阶数的bw程序，和李建平老师的进行对比可以发现2阶差距最小。
   
    从理论来看，不好意思，我看不懂，找不到理论联系，但我还是在文献最后找到了一句话：This procedure results in zero phase shift for all frequency。也就是从理论上来讲李建平教授程序那一套方法是没有位相误差的。如果把书上p208的结果当做滤波后的真实情况的话，李老师的程序与ncl自带程序padding后取2阶BW是最接近该波形的。
   附：matlab的butter函数与这个ncl子函数对比（本来目的是对数据大小进行检验，就那样吧，懒得动它了）
  总结：1.ncl的滤波个人觉得是丢失数据的，在此建议假如你要夏季数据的话，扔进滤波器的数据你得取5-9月或全年滤波
            2.李老师的程序在理论上是没有相位误差的，个人觉得是非常适合资料较短的数据的
            3.ncl与李老师Fortran版之间的结果是存在0.01数量级的误差的，个人站Fortran精度更好，但从结果来看，滤波后数据数量级在1，影响应该较小
  
最后给程序，各位可以看看改进
；================== 更新 =========================
这个函数放在do循环里时要注意，必须在do循环里面对fca、fcb进行赋值
试了一下发现，ncl子函数可能和Fortran的子函数一致，函数虚参的改变会返回至相应变量（就是说子函数里fca数据变化，返回来fca也变了）



来自群组: 南气院校友总会

不会游泳的fish

同学你好，想请教一下，如果我是看逐日演变的波动，有必要做高通滤波吗？

小懒C

不会游泳的fish 发表于 2020-11-13 16:25
同学你好，想请教一下，如果我是看逐日演变的波动，有必要做高通滤波吗？

那个应该不用吧，逐日变化你用原始资料就可以看出来波动逐日咋变的（就像分析天气图一样），具体我也不知道，你可以问问老师或者其他做过这个的同学。

猎人

感谢，先下载学习一下！{:eb502:}

何源

请问，该带通滤波是不是不适合做单一断点的，即高通或低通的滤波，利用李老师的Fortran程序计算的时候，发现报错。或者是否可以将fcb参数调为1.0即为高通，原序列减去高通序列既为低通？

小懒C

何源 发表于 2020-12-14 16:47
请问，该带通滤波是不是不适合做单一断点的，即高通或低通的滤波，利用李老师的Fortran程序计算的时候，发 ...

这个我没试过.，因为滤波器并不是理想的（也就是响应函数呈现阶梯那种），响应函数你就可以理解某波段的保留多少的信号这样。减去高通这个想法从理想上来看好像可以，但我感觉误差就不知道怎么说了，，你可以拿matlab试一试，因为matlab自己的butterworth滤波函数是可以高通低通带通带阻都可以实现的。

何源

小懒C 发表于 2020-12-14 17:09
这个我没试过.，因为滤波器并不是理想的（也就是响应函数呈现阶梯那种），响应函数你就可以理解某波段的 ...

的确是这样，我是初次接触滤波，关于滤波增加误差与不确定性这个我也是纠结了很久，和老板讨论了几次未果，我去试试matlab，谢谢。此外还想咨询关于滤波还有一个问题，是否是利用滤波提取变量的年际分量，再计算相关关系，就是两个变量的年际关系了？如果对提取后的变量计算滑动相关就可以探讨其年际关系的年代际变化？

小懒C

何源 发表于 2020-12-14 17:20
的确是这样，我是初次接触滤波，关于滤波增加误差与不确定性这个我也是纠结了很久，和老板讨论了几次未果 ...

可以这么做，但最好要弄清楚你提取的周期是否有讨论意义

一大碗年糕

太感谢了，学习学习

行星光影

最近要用到滤波 学习一下