二次IIR滤波器

RY DSP开发 2007/09/26

在声音信号处理中,经常需要设计滤波器,二次IIR滤波器由于其参数设计简单,运算量少,运用的最为广泛。这里介绍几种常用的滤波器的参数计算。

Peak Equalizer

Peak Equalizer的设计参数为中心频率f,最大增益db,形状参数Q。下面这个Python小程序通过这三个参数计算滤波器的参数a和b,取样频率为44100Hz。

01def peak(f, db, Q):
02    A = 10.0**(db/40.0)
03    w0= 2*math.pi*f/44100.0
04    alpha = math.sin(w0)/2/Q
05    b0 = 1+alpha*A
06    b1 = -2*math.cos(w0)
07    b2 = 1-alpha*A
08    a0 = 1+alpha/A
09    a1 = -2*cos(w0)
10    a2 = 1-alpha/A
11    b0 /= a0
12    b1 /= a0
13    b2 /= a0
14    a1 /= a0
15    a2 /= a0
16    a0 /= a0
17    return [b0,b1,b2],[a0,a1,a2]
 


高、低通滤波器

二次高、低通滤波器的滤波倾斜度为-12dB/oct ,只有一个参数:滤波频率f。

01def lpf(f):
02    a = [1.0, math.sqrt(2), 1.0]
03    b = [0.0, 0.0, 1.0]
04    ad = [0.0] * 3
05    bd = [0.0] * 3
06    fs = 44100.0
07    f = 2*math.pi*f
08    b[2] *= f*f
09    a[1] *= f
10    a[2] *= f*f
11    ad[0] = a[2] - a[1]*2*fs + 4*fs*fs*a[0]
12    ad[1] = a[2]*2 - 8*fs*fs*a[0]
13    ad[2] = a[2] + 2*fs*a[1] + 4*fs*fs*a[0]
14    bd[0] = b[2] - b[1]*2*fs + 4*fs*fs*b[0]
15    bd[1] = b[2]*2 - 8*fs*fs*b[0]
16    bd[2] = b[2] + b[1]*2*fs + 4*fs*fs*b[0]
17    t = ad[2]
18    return [bd[2]/t, bd[1]/t, bd[0]/t], [ad[2]/t, ad[1]/t, ad[0]/t]
 

01def hpf(f):
02    a = [1.0, math.sqrt(2), 1.0]
03    b = [0.0, 0.0, 1.0]
04    ad = [0.0] * 3
05    bd = [0.0] * 3
06    fs = 44100.0
07    f = 2*math.pi*f
08    b[0],b[2] = b[2],b[0]
09    a[1] *= f
10    a[2] *= f*f
11    ad[0] = a[2] - a[1]*2*fs + 4*fs*fs*a[0]
12    ad[1] = a[2]*2 - 8*fs*fs*a[0]
13    ad[2] = a[2] + 2*fs*a[1] + 4*fs*fs*a[0]
14    bd[0] = b[2] - b[1]*2*fs + 4*fs*fs*b[0]
15    bd[1] = b[2]*2 - 8*fs*fs*b[0]
16    bd[2] = b[2] + b[1]*2*fs + 4*fs*fs*b[0]
17    t = ad[2]
18    return [bd[2]/t, bd[1]/t, bd[0]/t], [ad[2]/t, ad[1]/t, ad[0]/t]
 




shelving滤波器

Shelving滤波器提升或者抑制高频或者低频的增益,有三个参数: 中心频率f,增益db,类型type。type为0时,对低频进行增益处理,为1时对高频进行增益处理。

01def shelving(f, db, type):
02    if type==1:
03        f = 22050 - f
04    Wb = 2.0*f/44100.0*math.pi
05    K = math.tan(Wb/2)
06    g = 10.0 ** (db/20.0)
07    V = g ** (1.0/2.0) - 1
08    a0 = 1 + 2*K + K*K
09    a1 = 2*K*K - 2
10    a2 = 1 - 2*K + K*K
11    b0 = a0
12    b1 = a1
13    b2 = a2
14    b0 = b0 + 2*V*K*(K+1)
15    b1 = b1 + 2*V*K*(2*K)
16    b2 = b2 + 2*V*K*(K-1)
17    b0 = b0 + V*V*K*K*1
18    b1 = b1 + V*V*K*K*2
19    b2 = b2 + V*V*K*K*1
20    b0 = b0 / a0
21    b1 = b1 / a0
22    b2 = b2 / a0
23    a1 = a1 / a0
24    a2 = a2 / a0
25    a0 = a0 / a0
26    if type==1:
27        b1 = -b1
28        a1 = -a1
29    return [b0,b1,b2],[a0,a1,a2]
 


源程序高亮显示代码公布

经过一段时间的调试,高亮度显示Python源代码的程序似乎已经没有什么BUG了,于是就给网站添加了一个功能: 用这段显示本网站的代码,不过由于安全方面的考虑,目前只能浏览如下3个源文件,今后将用这种方式发布一些Python写的小程序。

colorit.py : 高亮度显示Python源代码

showpy.py : 调用colorit.py显示网站的源代码

sudokuSolver.py : 数独游戏求程序


Python源程序的高亮度化

最近博客中的Python源程序越贴越多,但是由于不支持源程序的高亮度显示,因此阅读大段的程序比较困难。

很早就想在博客中加入源程序的高亮度显示功能。也找到了不少Python的库可以完美地完成这些功能。可是有一个问题:因为服务器速度慢、内存少,为了一个高亮度显示而载入一大堆库实在是负担太大。

于是我就着手自己写一个最最简单的源程序高亮度化的Python小程序。这个问题初看上去似乎比较困难,可是仔细思考之后才发现原来很简单。不到1个小时这个程序就完成了,而且也只有九十多行,处理速度也很快, 因此我就不考虑缓存结果,而是每次显示的时候都重新计算一遍。

下面是一个高亮度化的例子,可以看看效果:
※程序还是有一些小BUG,所以暂时不公开 (源程序已经公布 2007/9/22)

class Dummy(object):
def __init__(self):
self.n = 1000

def go(self):
for i in range(1000):
print "number %s:" % i

用回溯法解迷宫谜题

按照下面的规则,找到从开始方格(红色)到目标方格(绿色)的路径。

2
6
4
1
2
2
 
      
      
      
      
      
      
3
3
2
1
3
5

规则

  1. 同一个方格不能经过两次
  2. 行走方向只能为横竖方向,斜对角方向不能行走
  3. 上方和右方的数字表示那一行或者列中,所经过的方格的个数

例如上讲这个题目的解答为:
2
6 4 1 2 2
 


   


   
    
     
   
 
3
3
2
1
3
5

下面是一个10×10的谜题,这个题目有些难度,手工解决花了不少时间,于是就想用Python写一个程序解解看。

4
1
9
4
6
5
3
3
4
8
 
          
          
          
          
          
          
          
          
          
          
4
6
7
4
4
3
3
7
1
8


解这种题目最简单的方法就是回溯法,回溯法最简单的实现就是用递归。这里先给出完整的程序,然后再简单介绍一下编程思路。

01from sets import Set
02 
03class Maze:
04    def __init__(self, N):
05        self.maze = Set()
06        self.sx = [0]*N
07        self.sy = [0]*N
08        self.N = N
09             
10    def add(self, pos):
11        self.maze.add(pos)
12        self.sx[pos[0]] += 1
13        self.sy[pos[1]] += 1
14         
15    def remove(self, pos):
16        self.maze.remove(pos)
17        self.sx[pos[0]] -= 1
18        self.sy[pos[1]] -= 1
19         
20    def check(self, pos, X, Y):
21        x,y = pos
22        if self.sx[x] > X[x]: return False
23        if self.sy[y] > Y[y]: return False
24        if self.sx[x] == X[x]:
25            for i in xrange(x):
26                if self.sx[i]!=X[i]: return False
27        if self.sy[y] == Y[y]:
28            for i in xrange(y):
29                if self.sy[i]!=Y[i]: return False
30        return True
31 
32    def printout(self):
33        for i in xrange(self.N):
34            for j in xrange(self.N):
35                if (j,i) in self.maze: print "*",
36                else: print ".",
37            print
38 
39def near_pos(pos, N):
40    x,y=pos
41    if x-1>=0: yield x-1,y
42    if x+1<=N-1:yield x+1,y
43    if y-1>=0: yield x,y-1
44    if y+1<=N-1:yield x,y+1
45        
46def maze_step(maze, pos, target, N):
47    maze.add(pos)
48    if maze.check(pos,X,Y):
49        if pos == target:
50            maze.printout()
51            return
52        for npos in near_pos(pos, N):
53            if npos not in maze.maze:
54                maze_step(maze, npos, target, N)
55    maze.remove(pos)
56 
57X = [4,1,9,4,6,5,3,3,4,8]
58Y = [4,6,7,4,4,3,3,7,1,8]
59maze_step(Maze(10),(0,0),(9,9),10)
 

Maze类的maze属性是一个集合(Set),用来记录当前走过的所有的方格的坐标。sx和sy属性分别用来记录maze中每列以及每行的方格数。add和remove方法用来添加和删除某个方块,我们看到添加和删除方块的时候,同时更新sx和sy的相应坐标的值。

check方法用来检测当前的走法是否符合规则。pos参数为最后一个添加进去的方格坐标,也就是需要检查的对象。我们先看看这个方格所在的行和列的和是否满足条件:

        if self.sx[x] > X[x]: return False
if self.sy[y] > Y[y]: return False

由于入口在左上,而出口在右下,因此一旦某一行或者列的方格数等于了规则所给出的方格数,那么它左边的列或者上面的行就必须也符合规定的数目,因为如果要走回去改变方格数的话,那么就必须再经过一次已经等于规定数目的行或者列,这样就不能满足条件了。这个条件很关键,如果去除掉这个条件检测的话,那么等几个小时可能都出不来结果。

        if self.sx[x] == X[x]:
for i in xrange(x):
if self.sx[i]!=X[i]: return False
if self.sy[y] == Y[y]:
for i in xrange(y):
if self.sy[i]!=Y[i]: return False

maze_step完成递归回溯功能,注意在回溯(也就是函数返回)之前,需要通过执行remove方法删除掉对maze的修改。

最后这个程序的输出为:

* . . * * * . . . .
* . * * . * . . * *
* . * . . * * * * *
* * * . . . . . . *
. . * * * . . . . *
. . * . * . . . . *
. . * . * . . . . *
. . * . * * * * * *
. . * . . . . . . .
. . * * * * * * * *

方便起见,程序没有记录行走路径,只输出路径所经历过的方块,如果需要记录行走路径的话,那就需要在Maze类中再添加一个path属性,它是一个列表,在add和remove中分别对其进行添加和删除的操作即可。


爸爸回家看海月

今年8月份是我最高兴的一个月, 因为爸爸妈妈都回家来看我了。特别是爸爸,我已经有一年多没有和他见面了。在这一年多的时间里,我们只能靠网络聊天相互熟悉对方的声音,靠照片相互熟悉对方的相貌。因此当我看见爸爸风尘仆仆地走进来的时候,我就躲在妈妈身后,打量这个既陌生又熟悉的爸爸。

妈妈问我:“海月,这个人是谁?”

我扭扭捏捏,挤出一句“爸爸”,就逃到别的房间去了。

我能一眼就认出爸爸来,令他很欣慰,不过接下来我的表现,就让他有些着急了。我想,这么久你都没有来看我,绝对不能马上和你靠近乎,所以我对他的态度若即若离。先是一见到爸爸就假装哭,令他手足无措,只好拿好吃的来哄我。我看爸爸的态度很诚恳,就答应让他亲我了。

不过亲归亲,爸爸想抱我,我可是不让的。以至于好几次我们全家出去玩都是妈妈抱着我,把她累了个够呛。我越不让爸爸抱我,他就越想抱我,想尽各种办法逗我开心。这样自然而然地,我们父女之间的感情就一下子好了起来。爸爸也就在不知不觉中将我抱了起来。

爸爸抱我和别人不一样,他让我跳蹦极(我轻轻一跳,爸爸就把我高高地举起来),玩飞机(把我举在胸前,满屋子乱跑),这样的活动令我很兴奋和高兴,爸爸也是好久没有这么开心地运动了。