雪传说

      一路走来，无论结果如何，都能无悔这一次。
      不是xy，真的很佩服决心考研并且坚持下来的人们，祝你们顺利完成最后一步！
      祝大D姐姐，大D弟弟，马甲达人，dfl师兄，瘦子和我们班其他筒子们，亲爱的阿芳、巍哥、FoxFly、胖子、风华，以及所有其他奋斗至今的筒子们，马到成功！
      ＹＯＵ　ＡＲＥ　ＯＵＲ　ＨＥＲＯＥＳ．

      Knuth爷爷果然是超级大牛人- - //废话还用你说
      下午去理图，偶然瞟到一本TAOCP第四卷第四册，无比敬仰地翻了一下，一本就只有7.2.1一小节，我是当然看不懂了。其中7.2.1.7小小节讲的是历史和参考文献，从易经开始，很多篇幅讲八卦（是八卦，不是八卦...）和音乐。这一小小节后面只有三十几道习题（TAOCP的习题答案篇幅比正文多- -），我现在只记得一道：莱布尼茨于1666年发表On the Art of Combinations（也喜欢用Art= =），问为什么这是一个非常吉祥的年份？答案是，因为在罗马数字里，1666=MDCLXVI，其中M>D>C>L>X>V>I。
      刚瞟了几眼高爷爷的网站，写作风格果然很xxxx....在解释他为什么不用email感谢大家理解之后，还来了句"In return, I promise not to send unwelcome email requests to you."

我最早看到的提及高爷爷的文章，不知道出处
我心目中的编程高手
光华上看到的专门八卦高爷爷的文章
About Knuth
某人对TAOCP和Knuth的评论
建筑的永恒之道
从这里知道了The Stanford GraphBase
高德纳的二十年计划
The Stanford GraphBase的pdf下载地址
The Stanford GraphBase
最后，高爷爷的个人网站
http://www-cs-faculty.stanford.edu/~knuth/

      刚把网络课pro的新版文档发过去，瞬间就收到了老师的回复，原来这么晚还在看我们的报告...对比一下某些老板，sigh
      紧赶慢赶的，终于把所有东西完成了。周末写了四份实验报告，大概就是本科最后的实验课的实验报告了。医电的两份就算以我的标准来看也算是精品了//shy，DSP实验报告则无疑是垃圾中的垃圾...自我bs- -
      不知道这是不是第一门对拿A没有任何兴趣的课。卫星大概是第二门？不过从今天的考试情况来看还是有可能的sigh@@
      大四了，本来还想任性一下，放弃复习一门自己不喜欢的课的，结果终究还是无法违背自己的性格。
      不尽力，会很难受。（这什么破性格= =）
      就像我会花很长时间对每一篇报告进行排版，会再三对比封面标题放在那一行最美观，会手工为有几十节的文档添加书签，一字一句地敲。
      有些人则不同，他们秉承“能做就做，能省则省”的原则，需要做的事情，马上去做，可以不做的就忽略。两点看似矛盾，其实不然。我不知道完美主义的对立面是什么，我也不想谈什么主义，但是我需要学习，需要改善。如果一个人不再追求完美，如果一个人觉得差不多就可以了，那么他做事必然不会碍手碍脚，他可以立即开始做下一件事，而不必过分考虑这个开端是否是最好的；他不再有压力，他对自己的成果总是满意，或者不置可否的；他有了许多时间做自己的事情；他的生活轻松得多。
      反过来，是否完美主义就一定和低效、踌躇联系在一起呢？完美主义者在需要填写一个词语、设置一处标记、确定一个位置的时候，看到别人不假思索地填上一个显然不完美的答案，心里不由嗤笑。完美主义者最终花费心思给出了完美的答案，而这个答案除了自己却往往没人在乎。
      完美主义者有时不在乎是否有人在乎，完美主义者有时会告诫自己没有人在乎却也有时难以克服完美主义的心思。
      完美主义者总是想等到万事俱备的时候再顺理成章地开始做一件事情，哪想这世上从来没有过万事俱备，有的只是火烧眉毛。有篇文章说，喜欢规划To do list的人往往完不成自己的list，嗯，进程调度总是比中断处理复杂得多。   
      完美主义者每天都很忙碌，好像总是有做不完的事。完美主义者把别人在看电影、玩游戏的时间用在“正事”上，却好像也没有做出更多的“正事”。完美主义者的心理压力总是很大，于是很多人把“完美主义”和“抑郁”联系在一起。
      如果你认为我上面所提到的“完美主义者”侵犯了你的权益，那么请与我联系。
      以下是不相干的引申。
      做技术的通常都会有完美主义的倾向，所以
      1.海归往往做不过山寨，因为对完美的追求大大增加了研发成本，更重要的是延误了市场时机。
      2.技术出身的往往做不好管理，因为完美主义会导致对他人工作的不信任。
      以下不是不相干的引申，是不相干的冷笑话。
      k借用了我的洗衣卡，完后对我说，咱们楼的洗衣机果然坏了，我用之前显示剩余5次，用完变成6次了，bg我吧。
      我说，那个显示的5，其实是9。
      这个故事告诉我们，七段数码管的最小码距太小。

    我写了个prod(size(W))，它居然给了个warning说"NUMEL(x) is usually faster than PROD(SIZE(x))"，囧

    [转载请保留出处]
    http://snowflurry.ycool.com/post.3025548.html
    作者: snowflurry (http://snowflurry.cn/)

    最近做DSP作业，需要用到不少复杂的表达式化简，所以尝试了一下Matlab的符号运算。用起来还是很方便的，有兴趣的可以尝试一下: )

[符号运算能做什么]

    举个例子，用BZT方法设计IIR滤波器时，在选定原型滤波器之后，需要从H(s)得到H'(s)再得到H(z)，这个过程手算的话费时费力，用Matlab的符号运算功能只需要简单的几条命令：

 >> syms s z
 >> hs=1/(s+1);         % 以一阶Butterworth为例
 >> hsp=subs(hs,s/wp);  % lowpass to lowpass
 >> hz=subs(hsp,(z-1)/(z+1))

    wp是之前已经求出的，假设wp=1.23，那么上面最后一句会输出

 hz =
 1/(100/123*(z-1)/(z+1)+1)

    这就是最后得出的H(z)！当然，还可以通过各种方法来把它化简成易读的格式，后面再讲。

[符号运算的基本命令]

    下面我们就来解释一下刚才的程序。

 >> syms s z

    这一句定义了两个符号对象s和z。符号对象是Matlab数据类型的一种，至于到底是怎么回事大可不必理会，会用就行了: )

 >> hs=1/(s+1);

    输入原型滤波器的传递函数。由于s是符号对象，所以出来的hs也是符号对象，在Matlab的Workspace里用一个小正方体表示～

 >> hsp=subs(hs,s/wp);
 >> hz=subs(hsp,(z-1)/(z+1))

    这两句很关键，subs函数的作用是替换表达式里的变量，R=subs(S,new)会把S里的默认变量全部替换为new，也就是通常说的代入。subs函数的完整使用方法可以参阅帮助文档（doc subs）。

    但是，subs只是执行了替换，如何把替换结果进行化简呢？可以用simple函数：

 >> simple(hz)
 ans =
 123*(z+1)/(223*z+23)

    可以看到结果变成了两个分式相除的一般形式。Matlab在进行符号运算的时候会把数字也作为符号处理，像我们刚才的wp=1.23就会变成123/100，这样最后给出的是解析解，不会有任何精度损失。如果要把最终结果转换为小数表示形式的话，可以使用vpa函数，它的第二个参数还可以指定有效数字位数。

[一个实际例子]

    现在我们已经知道了符号运算到底是怎么回事，下面就以一个具体例子来熟悉它的使用方法。

    例题：DSP课本(E.C.Ifeachor & B.W.Jervis, Digital Signal Processing: A Practical Approach, 2nd Edition)第八章，Problem 8.11:

  The requirements for a certain lowpass filter are given below:
    passband              0-30Hz
    stopband edge         50Hz
    stopband attenuation  >40dB, at f>50Hz
    sampling frequency    256Hz
  Assuming that the filter has a Butterworth response, determine via the bilinear transformation the transfer function, H(z), of the filter. Obtain the realization diagram for the filter in cascade form using second- and/or first-order sections.

    首先，根据分析可以知道此题中需要以8阶Butterworth滤波器为原型。分析过程都是简单的计算，在此不再赘述，下面我们直接进入滤波器设计环节。

    Stage 1. 原型滤波器设计

    这一步里需要生成8阶Butterworth滤波器的传递函数H(s)。首先我们根据Butterworth滤波器的极点定义式求出所有极点{p(k)}，然后再把它们带入到式子H(s)=1/∏(s-p(k))中得到传递函数。程序如下：
   
  p=exp(j*pi*(2*(1:8)+7)/16); 
  % 求出所有极点。也可以直接用[z,p]=buttap(8)
  syms s z
  % 定义后面需要用到的符号s和z
  hs=1/sum(real(poly(p)).*(s.^(8:-1:0)));
  % 得到H(s)

    最后一句中，poly(p)的作用是计算以数组p中元素为根的多项式的系数，例如poly([1 2])返回[1 -3 2]，也就是多项式x^2-3x+2的系数。把poly(p)返回的结果与s^k(k=8,...,0)逐项相乘(.*)，再求和(sum)，就得到了H(s)的分母表达式。根据Butterworth滤波器规则，H(s)分子为1。

    real函数的作用是取poly(p)返回结果的实部。因为p中元素是以共轭对形式存在的，所以多项式系数肯定是实数，但是实际计算中精度有限引起的误差会造成poly(p)的返回值中带有很小的虚数部分，在这里将它去掉可以使之后的计算结果看起来更清爽。

    如果不用poly函数也可以很方便地得到H(s)，就是直接用前面的连乘式子H(s)=1/∏(s-p(k))：
    hs=1/prod(s-p);
    但是这样的话我就不知道怎样去掉系数中的虚部了……

    Stage 2. 双线性z变换

    这一步就是要用BZT变换得到最终的传递函数H(z)。程序很简单：

  wp=tan(pi*30/256);
  hsp=subs(hs,s/wp);
  hz=subs(hsp,(z-1)/(z+1));

    经过两步替代，已经得到了H(z)的表达式，下面来看怎么化简。

    首先用numden函数分别取出H(z)的分子和分母，这个函数会自动把H(z)化简成两个关于默认变量z的多项式相除的形式，然后分别返回分子和分母多项式：

  [b,a]=numden(hz);

    现在分子和分母多项式分别保存在变量b和a中，但是由于前面我们已经提到过的Matlab在进行符号运算的过程中会把小数全部做为分数进行处理，最后得到的b和a都包含非常大的整数系数。而根据我们的习惯，最后H(z)的表达式里分母中z最高次项的系数都要化为1，所以我们需要对b和a同除以a中最高次项的系数进行归一化：

  Ca=sym2poly(a);
  bfinal=b/Ca(1);
  afinal=a/Ca(1);

    sym2poly是一个很有用的函数，它可以输入符号多项式，返回一个包含其系数（按幂次从高到低排列）的数组。相应地也有poly2sym函数，用于从系数数组构建符号多项式。所以前面我们用到的
    hs=1/sum(real(poly(p)).*(s.^(8:-1:0)))
    也可以写成
    hs=1/poly2sym(real(poly(p)),s)

    现在已经得到了H(z)分子和分母各自的表达式，但是它们实在是太长了，因为所有的系数仍然是分数的形式。要转换为小数形式进行显示，就需要用我们之前提到过的vpa函数了：

  B=vpa(bfinal,4)
  A=vpa(afinal,4)

    其中第二个参数指定所有输出数值均包含4位有效数字。好了，现在就可以把输出结果抄下来了: )

    其实，如果只是要系数的话，用sym2poly(bfinal)比用vpa更方便，因为它直接把系数按顺序排列出来，不用再从长长的表达式里找幂次！

    Stage 3. 滤波器实现

    现在我们已经得到了IIR滤波器的传递函数，但是要通过2阶级联形式实现它的话，还需要计算每一个section的系数，也就是要把H(z)分解为四个二阶分式相乘的形式。H(z)的分子afinal的形式是a0*(z+1)^8，很容易分解，而bfinal给出的是8次多项式，要分解的话必须先求出它的根，也就是H(z)的极点。我们用roots函数：

  hzpoles=roots(sym2poly(afinal))

    roots函数可以直接从多项式的系数求出它的所有根，所以我们先用sym2poly来得到afinal的系数，作为roots的输入参数，得到的就是afinal的八个共轭复根：

 hzpoles =
    0.6551 + 0.5824i
    0.6551 - 0.5824i
    0.5396 + 0.4067i
    0.5396 - 0.4067i
    0.4755 + 0.2394i
    0.4755 - 0.2394i
    0.4467 + 0.0790i
    0.4467 - 0.0790i

    最后让每一对共轭复根对应的因式相乘，就得到了bfinal的四个二阶分式：

  Bsections=vpa(expand((z-hzpoles(1:2:7)).*(z-hzpoles(2:2:8))),4)
 
 Bsections =
  z^2-1.310*z+.7683
  z^2-1.079*z+.4566
  z^2-.9509*z+.2834
  z^2-.8934*z+.2058

    其中expand函数的作用是把多项式展开，用simple也可以。

[结语]

    我们只介绍了Matlab符号运算的一小部分内容，Matlab符号工具箱还有许多其它的强大的功能，比如可以直接执行Laplace变换和z变换，还能求微分方程的解析解，有兴趣的可以阅读Matlab帮助文档的Symbolic Math Toolbox部分。我刚刚在看的时候就发现了一个有趣的函数pretty，大家可以试试，比如pretty(hz)，呵呵。


附: 例题中用到的程序完整代码
（在MATLAB 7.6.0 (R2008a)中测试通过）
--------------------------------------------------------------------
% p811.m  snowflurry  2008.11.25
wp = tan(pi*30/256);
% Stage 1. Design of prototype filter
p = exp(j*pi*(2*(1:8)+7)/16);
syms s z
hs = 1 / sum(real(poly(p)) .* (s.^(8:-1:0)));
% Stage 2. The BZT
hsp = subs(hs, s/wp);
hz = subs(hsp, (z-1)/(z+1));
[b,a] = numden(hz);
Ca = sym2poly(a);
bfinal = b/Ca(1);
afinal = a/Ca(1);
clear s Ca a b hs hsp p wp;
B = vpa(bfinal, 4)
A = vpa(afinal, 4)
% Stage 3. Realization
hzpoles = roots(sym2poly(afinal));
Bsections = vpa(expand((z-hzpoles(1:2:7)) .* (z-hzpoles(2:2:8))), 4)
--------------------------------------------------------------------

      备份 & 给小喵

      firefox下载
      mozilla官网首页(www.mozilla.org)上就有链接。

      设置
      windows版的设置在工具下面，linux版的设置在编辑下面。
      一般需要设置如下内容：
      1.主要-启动-首页；
      2.主要-下载-总是询问保存文件的位置；（大多数人的习惯）
      3.内容-默认字体-中文和英文的最小字号。
      4.高级-网络-连接-设置，如果要设置代理的话（装了FoxyProxy就可以取消这里的设置了）。

      插件
      需要手动安装flashplayer，去Adobe的官网找。可能在这里。

      扩展
      安装扩展的方法：
      工具-附加组件-获取附加组件-浏览全部附加组件。一般不需要从其它地方下载插件。
      以下是我安装的扩展。
      1. Easy DragToGo: 拖拽链接打开网页，我现在已经离不了它了。
      2. FireGestures: 鼠标手势，同样，自从抛弃IE之后就已经习惯了用鼠标手势来关闭页面和撤销关闭。
      3. FoxProxy: 强大的代理管理工具，不用代理的人可以无视。
      4. Tab Mix Plus: 想通过双击关闭标签页？想让没浏览过的标签页用不同颜色显示？想自定义标签的排放顺序？TMP可以满足几乎所有的标签栏定制要求。
      其它常用扩展：
      IE Tab: 用IE内核打开网页，用于对付firefox显示不正常的网页，当然，它只支持windows版本；
      Flashgot: 下载管理器。
      其它在此不再赘述，可以在附加组件下载页分类查看最热门的扩展。Google一下，有很多介绍常用扩展的文章。

      这两天被输入法搞得不行了，总结一下。

      fcitx:
      用了两年了，以前和scim相比还是蛮好用的，现在大家都基于网络词库了，早停止开发的fcitx显然跟不上时代了。而且那个界面真的很难看，and 没有大字符集连囧都打不出(囧)  兼容性还是不错的。

      scim:
      重新装回scim，对Qt支持不错，GTK下有些问题，Gnome的Brower里也不能定位了。

      scim-bridge:
      装了scim-bridge，Qt和GTK下都没问题了，可是java程序不行……可怜的永中Office -.-
      scim拼音的操作还是一如既往地差，已经打出的拼音不显示，前面出错的话不好改，而且一行9个备选词无法调整，不知道土人我一眼看不过来那么多 + 超过5的数字不会盲打么……总之对scim很不满意。

      fitx:
      据说很不错，一直都想尝试下，每次都以失败告终。1家依赖的包太多，要求的版本太新，花了数个小时把它依赖的包和它依赖的包依赖的包和它依赖的包依赖的包依赖的包们都一个个找到，发现还是通不过。有若干个原有的包次版本号不满足要求，出于安全性考虑最终放弃它了。

      ibus:
      偶然看到的，正好又在装永中office的时候加了cn99的源，直接apt-get装上了，发现真的很好用，而且拼音的设置很丰富很强大。现在才0.1.1版，个人感觉很有前途XD

      ibus现在的版本支持GTK和Qt4，不过在Qt3下面倒也可以用，要设置环境变量QT_IM_MODULE=xim，所以启动qterm的话可以在命令行里用
      $ QT_IM_MODULE=xim qterm
      每次这样很不爽，Gnome面板上的launcher又不能识别这样的命令，尝试了N久，最后突发奇想：
      $ cd /usr/bin
      $ sudo mv qterm qterm.bin
      $ sudo vi qterm
      "qterm":
          QT_IM_MODULE=xim qterm.bin
      :x
      $ sudo chmod a+x qterm
      哈。
      某些GTK+程序下貌似也有问题，不能shift切换中英文。不过没关系，大不了以后不用PCManX了，反正它也不好用XD

.bashrc:
PS1='${debian_chroot:+($debian_chroot)}\[3[01;37m\]\u@\h\[3[00m\]:\[3[01;34m\]\w\[3[01;37m\]$\[3[00m\] '
这样前面的命令找起来就比较方便了，不会被一堆程序输出行淹没。
ps. 原来这里的颜色编码和bbs下的ansi一样哈哈

      中午在共青，在吃了好多糊不拉及的肉之后，我们的摊位终于没肉了，也没柴了，也没人了。趁着人多混乱，去另外一个摊位那里偷了一大盒羊肉，回来用报纸和满地的竹签生了火，偷偷地把一盒肉全烤了吃了=.=
      反正留下也是浪费嘛……
      ============== 我是黄昏分割线 ===============
      晚上去光华楼看神经网络，完全不懂……本来想大致浏览下关键的内容，忽略掉大部分数学推导的，结果发现这样的话就没几个字可以看了。这样的书真不好，概念全都是用公式来给的，真是的这么严谨作甚= = 于是去了理图，找到两本国外的教材，其中一个好像还有点好玩emmm 不过，总归是要仔细看了。
      下了两集MIT的Intro to Linear Algebra教学视频，那老师实在是太丑了-,-

      本来只是想去尝尝曲奇和匹萨的……遇见好多同学（叔叔的功劳），教育展也还蛮好玩的。开始填调查单，要写最想去的学校还是什么来着，我一个都不知道，后来Jc不知哪里搞了份宣传册，Aalborg University，看封面和名字都挺个性的就把它写上了。之后到了里面，拿到了不少宣传册，还在Aalborg那里留了email，不知道会不会收到什么东西...
      拿了件StudyInDenmark的Tshirt，现场的倒数第三件。Jc拿了倒数第二件，calla拿了the last one...自从穿上这衣服拍照的丹麦mm就总在对着我拍，囧。
      今天的丹麦人讲英文都很慢很清晰，赞。看了宣传册们搞得我又好想去……总是对欧洲有那么一种憧憬，田园式的浪漫与清淡。如果有机会可以考虑一下交流项目……好像不太可能@@
      想起大一时候的经典台词：北欧中心，我的最爱。
      好久没见yb了。
      提问：为什么我和Jc同时同地拿了同样的白色饮料，她的杯子里是可乐而我的是纯净水呢？
      回答：医生说我不能喝可乐。XD