一阶Butterworth带通滤波

肖凯

{:eb302:}{:eb302:}{:eb302:}

游泳的鱼

谢谢楼主！楼主可以讲一讲butter+filtfilt的用法吗？里头的参数设置多是信号处理的专业术语，不知在气象上该如何使用呢？

伽蓝鸟

游泳的鱼 发表于 2019-1-22 11:41
谢谢楼主！楼主可以讲一讲butter+filtfilt的用法吗？里头的参数设置多是信号处理的专业术语，不知在气象上 ...

比如说我有个逐日OLR数据，想滤出准双周振荡（10-25天），可以这么干
[b,a]=butter(4,[1/25,1/10]/0.5,'bandpass');
OLR12=filtfilt(b,a,OLR);
4是四阶的意思，然后是所需要的频段，带通滤波可以换成高通低通等其他滤波方式。
不要再下我那个脚本了，能算，但是只是一阶的，也就是Murakami的那篇文章里的

游泳的鱼

伽蓝鸟 发表于 2019-1-22 12:44
比如说我有个逐日OLR数据，想滤出准双周振荡（10-25天），可以这么干
=butter(4,[1/25,1/10]/0.5,'bandp ...

谢谢答复！我现在有逐年降水数据p，一共50年。想提取年际变化的信号，即把10年以下的信号。这么写对吗？
[b,a]=butter(3,10/50,'high');
p_high=filt(b,a,p);
我疑惑的是‘频段’的表达，你是用倒数表示的，这与时间序列长度无关吗？又除以了0.5表示什么呢？
谢谢！

游泳的鱼

伽蓝鸟 发表于 2019-1-22 12:44
比如说我有个逐日OLR数据，想滤出准双周振荡（10-25天），可以这么干
=butter(4,[1/25,1/10]/0.5,'bandp ...

谢谢答复！我现在有逐年降水数据p，一共50年。想提取年际变化的信号，即把10年以下的信号。这么写对吗？
[b,a]=butter(3,10/50,'high');
p_high=filt(b,a,p);
我疑惑的是‘频段’的表达，你是用倒数表示的，我理解的是想要的周期除以时间序列总长度。你又除以了0.5表示什么呢？
谢谢！

游泳的鱼

伽蓝鸟 发表于 2019-1-22 12:44
比如说我有个逐日OLR数据，想滤出准双周振荡（10-25天），可以这么干
=butter(4,[1/25,1/10]/0.5,'bandp ...

谢谢答复！我现在有逐年降水数据p，一共50年。想提取年际变化的信号，即把10年以下的信号。这么写对吗？
[b,a]=butter(3,10/50,'high');
p_high=filt(b,a,p);
我疑惑的是‘频段’的表达，你是用倒数表示的，我理解的是想要的周期除以时间序列总长度。你又除以了0.5表示什么呢？
谢谢！

伽蓝鸟

游泳的鱼 发表于 2019-1-24 15:52
谢谢答复！我现在有逐年降水数据p，一共50年。想提取年际变化的信号，即把10年以下的信号。这么写对吗？
...

[b,a]=butter(3,1/10/0.5,'high');
这样才对吧？你可以对滤波结果算个periodogram，看看是不是的确是10年以下的信号。
参数为啥这么设置，你还是老老实实看它的帮助文件吧。。
这几天感冒，迟复为歉。。

游泳的鱼

伽蓝鸟 发表于 2019-1-25 11:05
=butter(3,1/10/0.5,'high');
这样才对吧？你可以对滤波结果算个periodogram，看看是不是的确是10年以下 ...

对的，应该是[b,a]=butter(3,1/10/0.5,'high');
谢谢回复，保重~

伽蓝鸟

游泳的鱼 发表于 2019-1-25 16:41
对的，应该是=butter(3,1/10/0.5,'high');
谢谢回复，保重~

[Pxx,f]=periodogram(p_high,[],[],1);
plot(1./f,Pxx)
来确认一下的确是滤出了10年以下的高频信号？
此外呢，一种简单的滤出年际变率的方法是，对原始数据做9年滑动平均，再拿原始数据减去这个9年滑动的，就能得到较为高频的年际变率了。但是这个的响应函数可能不如Butterworth那么干脆。

游泳的鱼

伽蓝鸟 发表于 2019-1-25 19:46
=periodogram(p_high,[],[],1);
plot(1./f,Pxx)
来确认一下的确是滤出了10年以下的高频信号？

画出图来了，看不懂这个图呢