一个脚本测试：一阶Butterworth带通滤波

qxsw2016

根据教材“气候变率和诊断方法”（吴洪宝，吴蕾，2005，气象出版社）中first-order Butterworth带通滤波方法（p203-208），编写了一个脚本，对书上的例子进行测试，滤波迭代公式的y[0]和y[1]按照文献《一阶Butterworth递归式带通滤波器技术改进方案》（陈寅生等，1997），可取为0，或取y[0]=x[0]，y[1]=x[1]。关于书上所说逆序在计算一次和Murakami(1979)的原文里说还需要倒着算一遍("Then this output is reversed in time and processed again to obtain the final output")，并不理解其中的原因。

脚本：
import math
def Butterworth_filter(w0,w1):

    w2 = w0**2/w1
    deltat = 1
    delta_omega= 2*abs(sin(w1*deltat)/(1+cos(w1*deltat))-sin(w2*deltat)/(1+cos(w2*deltat)))
    omega02 = 4*sin(w1*deltat)*sin(w2*deltat)/((1+cos(w1*deltat))*(1+cos(w2*deltat)))

    deno = 4+2*delta_omega+omega02

    a=2*delta_omega/ deno
    b1=2*(omega02-4)/ deno
    b2=(4-2*delta_omega+omega02)/ deno
    return a,b1,b2

def cal_y(x,a,b1,b2):
    n = len(x)
    y = zeros(n)
    y[0] = x[0]
    y[1] = x[1]
    for i in arange(2,n):
        y = a*(x-x[i-2])-b1*y[i-1]-b2*y[i-2]

    return y

def x1(t):
    return 0.5*cos(2*math.pi*t/16)
def x2(t):
    return 1.0*cos(2*math.pi*t/45)
def x3(t):
    return 2.0*cos(2*math.pi*t/80)
def x(t):
    return x1(t)+x2(t)+x3(t)

n = 160
t = arange(n)
subplot(3,2,1)
plot(t,x1(t),'r--',linewidth=2)
title('$x1(t) = 0.5*\r{cos}(2 \pi t/16)$', fontsize=16)
xlabel('$time  \ t$')
ylabel('x1')

subplot(3,2,2)
plot(t,x2(t),'g--',linewidth=2)
title('$x2(t) = 1.0*\r{cos}(2 \pi t/45)$', fontsize=16)
xlabel('$time  \ t$')
ylabel('x2')

subplot(3,2,3)
plot(t,x2(t),'b--',linewidth=2)
title('$x3(t) = 2.0*\r{cos}(2 \pi t/80)$', fontsize=16)
xlabel('$time  \ t$')
ylabel('x3')
subplot(3,2,4)
plot(t,x(t),'k--',linewidth=2)
title('$x=x1(t)+x2(t)+x3(t)$')
xlabel('$time  \ t$')
ylabel('x')

w0=2*math.pi/45
w1=2*math.pi/35
c= Butterworth_filter(w0,w1)
a = c[0]
b1= c[1]
b2= c[2]
xx = zeros(n)
for i in arange(n):
    xx = x(i)
print a,b1,b2

y = cal_y(xx,a,b1,b2)
yy = cal_y(y[::-1],a,b1,b2)

subplot(3,2,5)
plot(t,yy[::-1],'k--',linewidth=2)
title('$filter \ x$')
#xlabel('x')
ylabel('x')

伽蓝鸟

正算一遍反算一遍是为了不要有位相移动，我之前也试过他们这个代码，可以看看我的帖子。

qxsw2016

伽蓝鸟 发表于 2020-12-1 21:31
正算一遍反算一遍是为了不要有位相移动，我之前也试过他们这个代码，可以看看我的帖子。

学习了，不知道二维怎么滤波

伽蓝鸟

qxsw2016 发表于 2020-12-3 09:36
学习了，不知道二维怎么滤波

啊？啥二维？

FDU.EDU

伽蓝鸟 发表于 2020-12-3 11:54
啊？啥二维？

可能意思是一个场变量，比如一年的风场，365*144*73这样的数据怎么滤波呢