爱吃面条

今天中午海月自己吃了一带面条,面条里面放了一个西红柿和一个鸡蛋。其实每次面条都是这么做的,但她好象吃不厌的感觉,每次都吃的有滋带味的。她边吃还边说:这是小猫的尾巴,我把尾巴先吃了。这是蚂蚁,我最喜欢吃蚂蚁了。还时不时的给我倒出来一点,说:妈妈你也吃,嚼,咽下去。就这样我全听她的指挥。
最近海月在儿童馆里渐渐习惯了,以前跳舞的时候总是躲起来,现在虽然还不想一起跳,但已经站在离大家比较近的位置了。而且还在馆里跑来跑去,丝毫没有拘谨的感觉。

好看的双眼皮怎么少了一个。

今天我们到儿童馆的时候,其它小朋友都还没有来,所以我想玩什么就玩什么。


能动噪声控制 - 测量系统响应

RY 隐藏 2008/06/25

在做具体的算法实验之前,先来测量一下实验装置的脉冲响应。在获得系统的脉冲响应之后,也可以先在PC上做Offline算法测试。关于实验装置的介绍请参看: 能动噪声控制 - 实验装置

我们用下面的系统框图来描述此实验装置中的各个声音回路,其中

  • P为噪声麦克风到误差麦克风之间的传递函数
  • C为噪声抑制扬声器到误差麦克风的传递函数
  • B为噪声抑制扬声器到噪声麦克风的传递函数



这些传递函数包括扬声器和麦克风的特性在内。为了测量B,C的脉冲响应,可以让噪声抑制扬声器播放TSP信号,然后分别对误差麦克风和噪声麦克风所获取的声音信号进行分析,计算出B和C的脉冲响应。详细方法请参看:用TSP信号测量系统的脉冲响应

因为P为两个麦克风之间的传递函数,这里无法使用TSP信号测量其脉冲响应。因此我采用自适应滤波器算法(NLMS)来测量它。具体做法如下:

  1. 从噪声扬声器播放白色噪声
  2. 假设噪声麦克风和误差麦克风的输入为x, y
  3. 让x通过自适应滤波器得到输出d
  4. 计算误差e = d - y
  5. 用误差e更新自适应滤波器的系数,使得误差e逐渐减少
  6. 当e稳定地小于一定的数值时,我们可以认为自适应滤波器同定了传递函数P
  7. 这时的自适应滤波器的系数就是P的脉冲响应

通过以上两种方法分别对P, B, C进行测量,得到如下结果

这个结果很有些意思,让我们来分析一下:

B和C的响应时间较长,其频谱在220Hz附近出现峰值。这个峰值应该是由于管道的固有频率引起的,我们可以将其理解为:扬声器所输出的声音中,220Hz的成分很容易在管道中传递。
误差麦克风虽然离扬声器比较近,但是其获取的信号没有噪声麦克风强。这是因为噪声麦克风位于管内,更容易获取管道中的各种反射的声音。

P则是一个比较单纯的延时系统,没有看到共振频率。这是因为两个麦克风都会获取共振信号,相互抵消。并且噪声麦克风获取的共振信号较强,因此P的频率响应在220Hz附近出现一个谷底。

我们看到B和C在3.5kHz以上有较高能量,这是因为我采用的低延迟DAC不能够对输出信号进行插值运算所导致的。

本文介绍的所有步骤都是在6713DSK上实现的,我将结果保存在DSP的内存中,然后通过CCS的data输出功能将结果保存到文件中,文件的格式如下,每一行都是一个32bit的HEX值:

1651 1 1e38 0 200
0x3D23C8EC
0x3C183F0C
0xBC340353 
....

为了能快速处理这样的数据,我写了一个Python的小程序自动将其转换为浮点数,并绘制频谱图:

01import sys
02from struct import *
03from binascii import *
04import pylab
05from numpy import *
06from Tkinter import *
07import glob
08import os.path
09files = glob.glob("*.dat")
10 
11def draw(filename):
12    def _plot():
13        r = []
14        for line in file(filename):
15            if line.startswith("0x"):
16                line = line.strip()
17                s = unhexlify(line[2:])
18                x = unpack(">f", s)[0]
19                r.append(x)
20        file(filename + ".txt","w").write(",\n".join([str(x) for x in r]))
21        pylab.subplot(211)
22        pylab.plot(r)
23        R = 10*log10(abs(fft.fft(r)))[:len(r)/2+1]
24        F = linspace(0, 4000, len(r)/2+1)
25        title.append(os.path.splitext(filename)[0])
26        pylab.title(title)
27        pylab.xlim(0, 512)
28        pylab.xlabel("Sample (8kHz sampling rate)")
29        pylab.grid(True)
30        pylab.subplot(212)
31        pylab.plot(F,R)
32        pylab.grid(True)
33        pylab.xlabel("Frequency (Hz)")
34        pylab.show()
35    return _plot
36 
37 
38root = Tk()
39for f in files:
40    b = Button(root, text=f, command=draw(f))
41    b.pack()
42mainloop()


长新本领了

我刚收拾好的房间,海月一会的功夫,就会又搞得到处是玩具,实在让我头疼。最近我想了个办法,就说:谁谁家的小孩很会收玩具,海月应该也可以,可看起来好象不行。妈妈先到别的房间,你来收玩具,看海月能不能收得干干净净。这么一说还真有效,海月就主动把我关在卧室里,自己收起客厅。收好后迫不及待地要我和她爸爸去看。收拾得还真干净,就是不仅仅是玩具,所有她认为碍眼的东西,一律扔进了玩具箱里。即使是这样还是受到了我们的欢迎,一起鼓掌以示鼓励。

爸爸妈妈带我故地重游,对于他们来说是故地。我可是一点印象也没有。


能动噪声控制 - 实验装置

RY 隐藏 2008/06/20

最近一直在忙这个能动噪声控制实验,实验在下面这个装置上进行。

播放器播放噪声从【噪声源喇叭】播放出,DSP通过【噪声检出麦克风】和【误差检出麦克风】的输入信号,计算出驱动【噪声抑制喇叭】的信号,使得其输出声音正好在【误差麦克风】处与噪声相互抵消。在【误差检出麦克风】处能够获得最大的噪声抑制,为了评价其它位置的噪声抑制情况,另外添加一个【性能评价麦克风】对其它位置的噪声进行录音。

这个实验的目的大致如下:

  • 先采用传统的filtered-x法进行实验,获取一些用于性能比较的参数
  • 然后用联立方程法(我们要研究的新算法)做实验,并和filtered-x的结果进行比较
  • 最后评价是否能在实际的工程环境(回转机器的降噪处理)稳定工作 

感觉此项目实在是很困难,即使在实验设备上获得较好的结果,在实际环境中还有许多需要考虑的因素。不过现在刚刚起步,目前已经调通了filtered-x和无反馈的理想情况下的联立方程法。今后的将陆续对此实验进行详细介绍。


善于雄辩的海月

上星期除了星期一没有出去,其它几天都在外面参加活动,渐渐地海月对一些活动场合也熟悉起来,能够自如地玩了。
现在海月非常会为自己辩解,她最常使用的一句话就是:"小孩子就是这样的。"每次她在哭闹过后,在总结经验的时候,总是最后来这么一句话。而且这句话是非常有出处的。其实这句话是大概在两个月前,海月不乖的时候,她爸爸用来安慰我的。结果被她灵活运用,时不时的用来为自己辩解,而且还总是强调:是爸爸这样说的。
海月在商店里会要买一些糖果,遇到这种情况不得不用些计策,比如:说那是药;或者说:家里有什么的。昨天我们一起去商店买东西,她爸爸拿着一袋糖,问她这是什么,她第一反应想说那是糖,然后马上改口说:是药。想想海月也挺可怜的,但是没办法,因为她还不能控制自己的馋嘴。

看天上月亮都出来了


用TSP信号测量系统的脉冲响应

RY DSP开发 2008/06/10

在测量系统的脉冲响应的时候,如果直接使用脉冲信号的话,信噪比会很低。因此一般都采用TSP(Time-Stretched Pulse)信号测量。TSP信号将脉冲的能量平均分散到时间轴上,将TSP信号输入到系统中获取其响应之后,再将分散的能量还原,重新构成系统的脉冲响应。能量分散的规则为:脉冲的相位和频率的平方成比例。重构脉冲响应的时候,只需要将系统的TSP响应与逆TSP信号进行卷积即可。

TSP信号和逆TSP信号的设计在频域进行。TSP的频域特性为:

H(k) = exp(j*4*pi*M*k^2/N^2)  0<=k<=N/2
H(k) = H(N-k)                 N/2

逆TSP的频域特性为

G(k) = 1/H(k)

其中N为TSP信号的长度,M为TSP实际数据长度的2倍(整数)。由于M为整数,我们很容易看到当k=N/2的时候,H(k)是一个实数,并且频域上有共轭对称的关系,所以对H(k)进行IFFT运算之后将得到一组实数,这就是TSP信号。

下面是Python求TSP信号的程序和TSP(512, 0.5)信号的波形图。

01from numpy import *
02def TSP(N, M):
03    M = int(N*M/2)
04    freq = array(range(N/2+1), "h")
05    H = exp(1j * 4*M*pi*(freq)**2.0/(N**2.0))
06    H = hstack( (H, conj( H[-2:0:-1] ) ))
07    x = fft.ifft(H)
08    x = real( x )
09    x = hstack( (x[N/2-M:], x[:N/2-M]) )
10    return x
11     
12if __name__ == "__main__":
13    import pylab
14    tsp = TSP(512, 0.5)
15    for item in tsp:
16        print item, ","
17    pylab.plot(tsp)
18    pylab.show()

参数M决定了TSP信号的有效参数,有效长度越大,信噪比越好。然而由于TSP信号在时间轴上的收敛需要一定的时间,有效长度太大的话,TSP信号的收敛特性将受到影响。

采用上述的TSP信号我实际测量了一下排气管的脉冲响应,测量用的DSP程序大致流程如下:

  1. 从音箱连续输出REPEAT次长度为N的TSP信号
  2. 被测量系统的脉冲响应长度为I
  3. 两次TSP信号之间的间隔(上一个TSP信号结束到下一个TSP信号开始的时间)为T, T>I
  4. 从麦克风输入系统对TSP信号的响应,长度为N+T,对REPEAT次TSP响应求平均
  5. 把平均之后的TSP响应和逆TSP信号就卷积,要计算逆TSP信号,只需要反转TSP信号即可
  6. 卷积结果中的N到N+I部分就是系统的脉冲响应


下面是这部分的示例程序

01#define I 512
02#define N 512
03#define T 1024
04#define L (N+T)
05#define RI (N+L-1)
06#define REPEAT 10
07 
08extern const float TSP_DATA[];
09 
10float TSP[N];
11float TSP_INV[N];
12float TSP_RESPONSE[L];
13float IMPULSE_RESPONSE[I];
14int tsp_index = 0;
15int tsp_count = 0;
16 
17void InitTSPData(){
18    int i;
19    for(i=0;i<N;i++){
20        TSP[i] = TSP_DATA[i];
21        TSP_INV[i+1] = TSP_DATA[N-1-i];
22    }
23    TSP_INV[0] = 0;
24    for(i=0;i<L;i++){
25        TSP_RESPONSE[i] = 0;
26    }
27    tsp_index = 0;
28    tsp_count = 0;
29}
30 
31void conv(float x[], float y[], float z[], int lenx, int leny, int startz, int lenz)
32{
33    int k, j;
34    int l = lenx+leny-1;
35    int min, max;
36    float s;
37    for(k=startz;k<startz+lenz;k++){
38        s=0;
39        min = (k-leny+1<0) ? 0 : (k-leny+1);
40        max = (k<lenx-1) ? k : (lenx-1);
41        for(j=min;j<=max;j++){
42            s += x[j]*y[k-j];
43        }
44        z[k-startz] = s;
45    }
46}
47 
48void main()
49{
50    InitTSPData();
51}
52 
53void DoSample()
54{
55    float out;
56    TSP_RESPONSE[tsp_index] += InputSample();
57     
58    if(tsp_index < N){
59        out = TSP[tsp_index];
60    }
61    else{
62        out = 0;
63    }
64    tsp_index ++;
65    if(tsp_index == L){
66        tsp_index = 0;
67        tsp_count ++;
68        if(tsp_count == REPEAT){
69            for(i=0;i<L;i++){
70                TSP_RESPONSE[i] /= tsp_count;
71            }
72            conv(TSP_RESPONSE, TSP_INV, IMPULSE_RESPONSE, L, N, N, I);
73            asm(" nop");
74        }
75    }
76    OutputSample(out);
77}

 

用此程序实际测量的排气管脉冲响应图如下:


点灯游戏及其求解

点灯游戏的规则:有个N行N列的灯板,当你开关其中一盏灯,它和上下左右的灯的状态全部反转,目标是将全部的灯点亮。

这个程序采用递归回溯法对任意状态的灯板进行求解。游戏中左点灯将按照游戏规则动作,右点灯只对一盏灯进行操作,以便您测试任何状态的自动求解。

我们很容易知道,对于每盏灯至多只需要点击一次即可以达到目标状态,因为点击偶数次相当于没有点,点奇数次相当于点一次。clickBoard 保存的就是这个点击信息。-1表示尚未处理,0表示不点,1表示点。当某个灯及其周围所有与其连通的灯都被处理之后,这个灯的状态就不能再被改变了,因此可以通过判断其状态来确定当前的 clickBoard 中的点击信息是否正确,函数CheckError(self, pos)就是进行这个判断。

LightBoard类定义了灯板的所有操作,以及上述所有求解所需的辅助工具。真正的求解部分在函数 Solve(board) 中进行。它依次对board中的每个待检测的灯进行如下动作:

  • 如果点击此灯没有造成错误,则递归调用Solve
  • 如果Pass此灯没有造成了错误,则递归调用Solve
  • 如果上述都失败,则表示当前状态无解,回溯

下面是源程序,关于源程序的讨论请去 {CodeLibrary}点灯游戏及其求解

001def InitTable(size, value):
002    return [[value for i in range(size)] for j in range(size)]
003 
004class LightBoard:
005    def __init__(self, size, status=""):
006        self.clickBoard = InitTable(size, -1)
007        self.lightBoard = InitTable(size, 0)
008        self.toCheck = []
009        for i in range(size):
010            for j in range(size):
011                self.PushToCheck((i,j))
012        for pos, value in enumerate(status):
013            if value == "1":
014                self.lightBoard[pos/size][pos%size] = 1
015        self.size = size
016 
017    def PopToCheck(self):
018        if len(self.toCheck) > 0:
019            pos = self.toCheck[-1]
020            del self.toCheck[-1]
021        else:
022            pos = None
023        return pos
024 
025    def PushToCheck(self, pos):
026        self.toCheck.append(pos)
027 
028    def ConnectedLight(self, pos):
029        x, y = pos
030        yield pos
031        if x - 1 >= 0: yield x-1,y
032        if x + 1 <self.size: yield x+1,y
033        if y -1 >= 0: yield x, y-1
034        if y+1 < self.size: yield x, y+1
035 
036    def NotCheckLight(self, pos):
037        return [(x,y) for x, y in self.ConnectedLight(pos) if self.clickBoard[x][y] == -1]
038 
039    def SetCheck(self, pos, value):
040        x, y = pos
041        self.clickBoard[x][y] = value
042        for x, y in self.ConnectedLight(pos):
043            self.lightBoard[x][y] = 1 - self.lightBoard[x][y]
044 
045    def Pass(self, pos):
046        x, y = pos
047        self.clickBoard[x][y] = 0
048 
049    def CancelPass(self, pos):
050        x, y = pos
051        self.clickBoard[x][y] = -1
052 
053    def CheckError(self, pos):
054        for p in self.ConnectedLight(pos):
055            x, y = p
056            if self.lightBoard[x][y] == 0 and len(self.NotCheckLight(p)) == 0:
057                return True
058        return False
059 
060def Solve(board):
061    clickpos = board.PopToCheck()
062    if clickpos == None:
063        return board.clickBoard
064    pos = board.NotCheckLight(clickpos)[0]
065    board.SetCheck(pos, 1)
066    if not board.CheckError( pos ):
067        r = Solve(board)
068        if r: return r
069    board.SetCheck(pos, -1)
070    board.Pass(pos)
071    if not board.CheckError( pos ):
072        r = Solve(board)
073        if r: return r
074    board.CancelPass(pos)
075    board.PushToCheck(clickpos)
076    return False
077 
078from Tkinter import *
079 
080class LightButton(Button):
081    def __init__(self, parent, pos, board):
082        Button.__init__(self, parent, width=2)
083        self.board = board
084        self.status = 0
085        self.x, self.y = pos
086        self.SetColor()
087        self.bind("<Button-1>",self.onClick)
088        self.bind("<Button-3>",self.onRightClick)
089 
090    def SetColor(self):
091        if self.status == 0:
092            self.config(bg="#cccccc")
093        if self.status == 1:
094            self.config(bg="#ffcccc")
095 
096    def onClick(self,event):
097        toclick = [(0, 0), (-1,0), (1,0), (0, 1), (0, -1)]
098        self.config(text="")
099        for x,y in toclick:
100            try:
101                b = self.board[(self.x-x, self.y-y)]
102                b.status = 1 - b.status
103                b.SetColor()
104            except:
105                pass
106 
107    def onRightClick(self, event):
108        self.status = 1 - self.status
109        self.SetColor()
110 
111if __name__ == "__main__":
112    from tkSimpleDialog import askinteger
113    def SolveIt():
114        boardStatus = ["0"] * Size * Size
115        for pos, light in buttons.items():
116            if light.status == 1:
117                boardStatus[pos[1]*Size+pos[0]] = "1"
118        clickBoard = Solve(LightBoard(Size, "".join(boardStatus)))
119        for x, r in enumerate(clickBoard):
120            for y, c in enumerate(r):
121                if c == 1:
122                    buttons[(y,x)].config(text="*")
123    def Clear():
124        for button in buttons.values():
125            button.status = 0
126            button.SetColor()
127    root = Tk()
128    root.title("Light Solver")
129    f = Frame(root)
130    f.pack(side = TOP)
131    Button(f, text="Solve", command=SolveIt).pack(side=LEFT)
132    Button(f, text="Clear", command=Clear).pack(side=LEFT)
133    buttons = {}
134    Size = askinteger("Light Solver", "Please Input Board Size", initialvalue=5)
135    for y in range(Size):
136        f = Frame(root)
137        f.pack(side = LEFT)
138        for x in range(Size):
139            pos = x, y
140            b = LightButton(f, pos, buttons)
141            b.pack()
142            buttons[pos] = b
143    mainloop()


海月上了一个新的台阶

妈妈给我换新发型了。

因为海月经常在公园的草丛中玩,为了不让蚊子亲她,我们特意给她买了装在小盒子里可以带在身上的除蚊药。她爸爸给她挂在身上的时候开玩笑地说:"海月带上这个,蚊子就不咬海月,咬爸爸妈妈了。"海月听了说:"我挡着爸爸妈妈,把蚊子赶走。"我们听了都觉得很诧异,没想到她能说出这样的话来。
最近几天明显感觉得出来海月又懂事了不少,渐渐地可以聊天,发表自己的意见,反驳我的话,还经常给我讲故事听。我也感觉比刚回日本的时候轻松了很多。
但并不意味着不忙了,某种意义上来说是更忙些。因为日本的幼儿园,托儿所每周都会有固定的开放日,就是学龄前的小孩子可以自由进去玩。我家附近的儿童馆每周三也有活动。另外区役所还专门给外国人开设了每月一次的交流活动。这样一来就算是遇上下雨的天气,每周也至少可以出去参加一次活动,当然都是免费参加的。只有以这种方式让海月多和其它的小朋友接触一下,扩展一下视野,为顺利进入日本的幼儿园做好准备。