一阶Butterworth带通滤波

伽蓝鸟

不知道怎么取消贡献，别下载附件了啊喂!!后来发现butter和filtfilt这俩内置函数就足够搞定了。。好久不登气象家园，我都忘了还有这个帖子。。抱歉抱歉。。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         最近的研究兴趣转向了季节内振荡，需要做带通滤波。气象家园网站上好像没多少MATLAB带通滤波的资料（至少二爷的MATLAB精华帖里没有）。看文献的时候似乎蛮多前辈都用过first-order Butterworth带通滤波。很多气象统计的教材都介绍了这种滤波方法（e.g. 黄嘉佑和李庆祥，2014，气象出版社；吴洪宝，吴蕾，2005，气象出版社）。原始文献应该是Murakami (1979, Mon Wea Rev)。这里也就按照Murakami(1979)用MATLAB编写了滤波代码。

      NCL里有类似的函数，但我不会用。MATLAB内置函数也有butter，明天就去学。      Fortran版的butterworth滤波可以在这里找到http://www.seismosoc.org/Publica ... 55-esupp/bfilter1.f
      我蛮久不用Fortran了，本想把它直译成MATLAB但没成功。很早前二爷告诉过我有软件可以把Fortran翻译成MATLAB，这回再问他说忘了。。啧。。
      Murakami(1979)的原文里说还需要倒着算一遍("Then this output is reversed in time and processed again to obtain the final output")，我不确定是否理解正确了，应该是的吧。如果有啥错误，请不吝赐教，先谢过啦~稍微做了点测试，应该是可以的。
      测试1：吴洪宝和吴蕾，2005，p208      这本书气象家园里有，我就不再发了，也介绍了此滤波方法，其实就是翻译了下Murakami(1979)

      测试2：Qiao et al. (2015, Adv. Atmos. Sci.)

      乔云亭老师教的我们气候统计，她当时在课上介绍的也是这种方法。基本上是一样的，但在数据开头的地方还是有点差异。我不知道原因，比如，是因为递归开头用零补齐的问题？还是计算精度的问题？又或者我码错了代码？一定要倒着再算一遍，不然的话差异很大。Murakami(1979)关于倒着算的描述有点含糊，我不确定是不是理解对了。另外，在110E, 18N并不是OLR数据的整点啊。。最接近的格点是110E, 17.5N。

      测试3：Sun and Ding (2008, Adv. Atmos. Sci.)
      除了数据开头那里大小稍微有点差，曲线形状还是蛮像的。感觉在滤波的时候应该往前往后多取一些数据，尽量避免端点的影响。

      测试4：Mao and Chen (2005, J. Clim.)
      这篇文章用的是小波来滤波，所以结果不完全一样也是可以预期的。但效果还蛮好的，走势很接近。
      
      测试5：和EEMD（集合经验模分解）做了下对比，数据是2010年南海地区的850hPa纬向风，是我之前做过的一个分析。EEMD是自适应的滤波器，不用预先给定阈值，这里的Butterworth分别设了12-24day和30-60day，相当接近了。
     EEMD的资料之前上传过，这方法计算太耗时了。。

      如果嫌麻烦，也可以用两个滑动相减来设计带通滤波器，气候统计的书上也有讲，又比如Cao and Wu (2017, AOSL)。图就不放上来了，不是很满意，但也可以用来测试结论的鲁棒性。
      不是很经常逛气象家园的网站，倒是会看MATLAB交流群。如果哪位朋友发现程序有错，可以在那里找到我，或者用扣扣邮箱793209554@qq.com，谢谢~
       啧……午睡之后看了下MATLAB内置函数butter，人家早就算好了。。果然遇到问题先看内置函数啊。。butter+filtfilt就能搞定了。。

伽蓝鸟

排版怎么都弄不好。。放弃了。。请凑合着看吧。。

liyf

不知道怎么给你送金钱。暂时给你点个赞！

看云淡风清

不错的帖子，学习下！

lleoiu

写的不错，给你点个赞！

weinihou

有理有据！

（不过果然应该先看说明，都知道有函数了还去自己编写的话，只能安慰自己“加深巩固理解和应用”了……）

zhangsanbao

学到了，谢谢楼主

zdcgl

不错，谢谢楼主分享！

学生党少华

                                 看下

学生党少华

看看