圖靈機 電腦是如何計算的？它的原理是什麼？

發表： 2011 年 07 月 28 日 | 作者： vinsonhsieh | Filed under: 計算機相關 | 發表留言

電腦是如何計算的？它的原理是什麼？

它為什麼和希爾伯特第十個問題、康托爾集合論、羅素悖論、哥德爾不完備定理及布林邏輯等息息相關。

自從20世紀30年代以來，圖靈機、計算這些重要的概念在科學的天空中就一直閃爍著無限的光彩。尤其是近年來量子電腦、生物電腦、DNA計算等領域的創新工作引起了世人的廣泛關注。我們不禁問這樣的問題，國外究竟為什麼能發明出這些各式各樣的電腦呢？這些意味著什麼呢？其實這一切的源頭都來源於計算理論。國內在介紹計算理論方面的教材雖然有不少，但一般都比較深奧難懂。所以我覺得很有必要對這些內容進行科普。於是嘗試寫下這麼一篇文章，希望我的文章能夠讓你更加清楚、透徹的理解圖靈機、計算等等一些基本而重要的概念，並洞悉到這些概念的本質和深遠涵義。
本篇文章大體上可以分成四部分：首先給大家講一講關於圖靈、哥德爾等科學家的故事；然後正式引入圖靈機的概念，為了對這個概念有比較直觀的理解，我採用了一個人工生命：“小蟲”的比喻來敘述。接下來，文章介紹了跟圖靈機有關的概念：什麼是類比，什麼是“萬能電腦”等等；最後是關於圖靈停機問題的探討，我個人認為很有可能未來對科學的重大突破都來源於對圖靈停機問題的深入理解。在行文過程中，我除了用自己的方式介紹一些現有的基本概念之外（為了儘量表達得清楚明白，我不得不放棄理論論證的嚴格性），還探討了很多我認為別人沒有探討的問題，這些問題多是我自己的思考結果，而它們沒有經過科學的驗證。在這部分內容上我都標上了*號，希望你能有選擇的看待這些問題和觀點。

一、故事

任何科學思想、科學概念的誕生都有它的背景，在背景中往往有很多迷人的故事。關於計算理論可以追溯到1900年，當時著名的大數學家希爾伯特在世紀之交的數學家大會上給國際數學界提出了著名的23個數學問題。其中第十問題是這樣的：存在不存在一種有限的、機械的步驟能夠判斷“丟番圖方程”是否存在解？這裡就提出來了有限的、機械的證明步驟的問題，用今天的話說就是演算法。但在當時，人們還不知道“演算法”是什麼。實際上，當時數學領域中已經有很多問題都是跟“演算法”密切相關的，因而，科學的“演算法” 定義呼之欲出。之後到了30年代的時候，終於有兩個人分別提出了精確定義演算法的方法，一個人是圖靈，一個人是丘奇。而其中圖靈提出來的圖靈機模型直觀形象，於是很快得到了大家的普遍接受。
不知道你是否聽說過圖靈這個名字。可能有些人知道牛頓，知道愛因斯坦，甚至知道馮諾依曼，但不知道圖靈。然而圖靈的貢獻絕對不亞於這些科學大師。圖靈最大的貢獻就是把演算法這樣一個基本的、深刻的概念用他的圖靈機模型講清楚了。正是因為圖靈奠定的理論基礎，人們才有可能發明20世紀以來甚至是人類有史以來最偉大的發明：電腦。因此人們稱圖靈為：電腦理論之父。
圖靈生活的年代經歷了第二次世界大戰。在二戰期間他曾經為英國政府效力成功破譯了德國的密碼，因而為英國做出了突出貢獻。其實也正是因為二戰，英國政府才肯掏錢讓圖靈製造最原始的電腦，當然這種電腦是專門用來破譯密碼用的，而不是我們現在用的通用電腦。（有一部片子叫《密碼迷情》英文名是《enigma》就是根據圖靈當時破譯德國密碼的故事改編的，大家有興趣可以去找一找。）
圖靈這個人很古怪，只喜歡自己一個人悶頭研究，不喜歡與別人交流。並且據說他還是一個同性戀者。要知道在當時的英國，同性戀行為可是大逆不道的。最後，在他事業剛剛達到頂風的時候，他自殺了。為了紀念這個偉大的學者，電腦界設立了最高榮譽獎：ACM圖靈獎。
圖靈機的產生一方面奠定了現代數位電腦的基礎（要知道後來馮諾依曼就是根據圖靈的設想才設計出第一台電腦的）。另一方面，根據圖靈機這一基本簡潔的概念，我們還可以看到可計算的極限是什麼。也就是說實際上電腦的本領從原則上講是有限制的。請注意，這裡說到電腦的極限並不是說它不能吃飯、掃地等硬體方面的極限，而是僅僅就從資訊處理這個角度，電腦也仍然存在著極限。這就是圖靈機的停機問題。這個問題在圖靈看來更加重要，在他當年的論文中，其實他是為了論證圖靈停機問題才“捎帶手”提出了圖靈機模型的。
提到了圖靈停機問題，我不禁又要提一提哥德爾定理、羅素悖論、康托爾(Georg Cantor)的集合論等等一系列大事兒。早在19世紀末的時候，康托爾為集合論做了奠基性的研究。要知道，數學雖然五花八門，但是人們發現，運用集合這個概念可以概括所有的數學，也就是說集合是一切數學的基礎。因而如果為集合論奠定了公理化的基礎，也就等於為數學奠定了基礎。康托爾就是做了這方面的貢獻。另外，他為了證明實數的個數比自然數多這個結論，發明瞭一種被稱為“對角線刪除”的證明方法。沒想到的是，這個方法影響非常深廣，直到後來的圖靈停機問題、哥德爾定理其實都是該方法的不同延伸。
19世紀末的人們忙於為基於集合論的數學建立公理體系大廈。然而正當這座大廈即將完工的時候，一件可怕的事情發生了，羅素提出來的羅素悖論粉碎了數學家的夢想。
羅素在一九Ｏ一年發現了關於自我指涉的集合論悖論，羅素指出，集合可分為兩類：一類以自身為元素，另一類則不以自身為元素。假設所有不以自身為元素的集合所構成的一個集合，稱做「羅素集」。

由此便引出這樣一個問題：該集合是否以自身為元素？

於是推出：如果羅素集以自身為元素，則有羅素集不以自身為元素；如果羅素集不以自身為元素，則有羅素集以自身為元素。

這又是熟悉的Ｐ←→非Ｐ的形式。大家要昏頭了嗎？羅素悖論的通俗的語義版本是理髮師悖論：在薩維爾村有一個理髮師，他掛出了一塊招牌規定著：「我給而且只給村民中不給自己刮鬍子的人刮鬍子。」於是有人就問他：「你給不給自己刮鬍子呢？」無論這個理髮師怎麼回答都會產生矛盾，如果他不給自己刮鬍子，那麼按照招牌，他應該給自己刮鬍子；如果他給自己刮鬍子，按照招牌所言，他只給村裡不給自己刮鬍子。

如果你嘗試回答這個問題就會發現奇怪的事情：這個問題本身似乎是不可能的 ！正是因為這種奇怪的邏輯，哲學家羅素才顛覆了整個數學大廈的基礎 ！

因為集合論中存在著矛盾，所以整個數學體系存在著根本性的矛盾，然而具有諷刺意味的是，數學一直以嚴格著稱 ！這種感覺就好像正當你得意洋洋的時候突然有人狠狠打了你一個耳光 ！

人們在一陣慌亂之後開始逐漸穩下腳步尋找避免羅素悖論產生的方法，並總希望通過細心的選擇，數學公理能夠將類似羅素悖論這樣的怪物一勞永逸的排除精確數學體系之外。這就是後來希爾伯特提出來一套數學綱領的原因，他希望找到一套公理體系能夠排除悖論，並挽救純粹而美麗無暇的數學。雖然希爾伯特沒能完成他的夢想，但他堅信夢想是對的。
然而，沒過多少年，一個名不見經傳的年輕邏輯學家哥德爾提出來的定理卻徹底粉碎了希爾伯特的夢。這就是後來著名的哥德爾定理 ！該定理大致上說：任何一個數學的公理化體系都不是“完美的”。換個角度說，任何數學公理化系統都是死的，它總需要人為地從外界注入新的公理進去才能讓它日趨完善，而它自己並不能完全自動避免矛盾產生。哥德爾定理可以說整個扭轉了人們的世界觀。因為被認為最“完美”最“純粹”的數學都是不完全的，那麼純粹完美的世界也應該不存在。哥德爾定理還指出了“理性”和“分析”方法的極限。這才是後來人們步入到了“綜合”時代的部分原因。更有趣的是，哥德爾用來證明他定理的方法正和康托爾證明實數比自然數多、圖靈停機問題以及羅素悖論的方法是一脈相承的
所有這些都是20世紀上半個世紀發生的大事，後來發生了什麼？電腦出現了、資訊時代來臨了，似乎科學技術是萬能的，它們總會改善我們的生活，滿足我們的欲望。漸漸的，人們似乎淡忘了圖靈、歌德爾、康托爾等等大師們的思想了。然而近年來隨著複雜性科學的研究，人們卻又有了重拾這些相對古老而根本問題的跡象了 ！

二、圖靈機

下面言歸正傳，我們開始講圖靈機的概念。我先把圖靈機的模型給你，雖然有些無趣，不過請堅持看下去，我會在下面運用大家比較好理解的形式重新解釋的。在這裡你僅僅需要認識它的輪廓。一個圖靈機是形如下麵的一個裝置：

這個裝置由下面幾個部分組成：一個無限長的紙帶，一個讀寫頭。（中間那個大盒子），內部狀態（盒子上的方塊，比如A,B,E,H），另外，還有一個程式對這個盒子進行控制。這個裝置就是根據程式的命令以及它的內部狀態進行磁帶的讀寫、移動。
它工作的時候是這樣的：從讀寫頭在紙帶上讀出一個方格的資訊，並且根據它當前的內部狀態開始對程式進行查表，然後得出一個輸出動作，也就是是否往紙帶上寫資訊，還是移動讀寫頭到下一個方格。程式也會告訴它下一時刻內部狀態轉移到哪一個。
具體的程式就是一個清單，也叫做規則表，是這樣的：

當前內部狀態s	輸入數值i	輸出動作o	下一時刻的內部狀態s’
B	1	前移	C
A	0	往紙帶上寫1	B
C	0	後移	A
…	…	…	…

因此，圖靈機只要根據每一時刻讀寫頭讀到的資訊和當前的內部狀態進行查表就可以確定它下一時刻的內部狀態和輸出動作了。
圖靈機就是這麼簡單！不可思議吧？而只要你變化它的程式（也就是上面的規則表），那麼它就可能為你做任何電腦能夠完成的工作。因此可以說，圖靈機就是一個最簡單的電腦模型 ！
也許，你會覺得圖靈機模型太簡單，怎麼可能完成電腦的複雜任務呢？問題的關鍵是如何理解這個模型。

三、如何理解圖靈機？

1、小蟲的比喻
我們不妨考慮這樣一個問題。假設一個小蟲在地上爬，那麼我們應該怎樣從小蟲處理資訊的角度來建立它的模型？

首先，我們需要對小蟲所在的環境進行建模。我們不妨就假設小蟲所處的世界是一個二維的世界位於一個無限長的紙帶上(小蟲只能前進或後退)，這個紙帶有一些小的方格，而每個方格都僅僅只有黑和白兩種顏色。

很顯然，這個小蟲要有眼睛或者鼻子或者耳朵等等感覺器官來獲得世界的資訊，我們不妨把模型簡化，假設它僅僅具有一個感覺器官：眼睛，而且它的視力短得可憐，也就是說它僅僅能夠感受到它所處的方格的顏色。因而這個方格所在的位置的黑色或者白色的資訊就是小蟲的輸入資訊。
另外，我們當然還需要為小蟲建立輸出裝置，也就是說它能夠動起來。我們仍然考慮最簡單的情況：小蟲的輸出動作就是往紙帶上前爬一個方格或者後退一個方格。
僅僅有了輸入裝置以及輸出裝置，小蟲還不能動起來，原因很簡單，它並不知道該怎樣在各種情況下選擇它的輸出動作。於是我們就需要給它指定行動的規則，這就是程式！

假設我們記小蟲的輸入資訊集合為

I={黑色，白色}

它的輸出可能行動的集合就是：

O={前移，後移}

那麼程式就是要告訴它在給定了輸入比如黑色情況下，它應該選擇什麼輸出。

因而，一個程式就是一個從I集合到O集合的對應。我們也可以用清單的方式來表示程式，比如：

程式1：

輸入  輸出
黑色  前移
白色  後移

這個程式非常簡單，它告訴小蟲當讀到一個黑色方格的時候就往前走一個方格，當讀到一個白色方格的時候就後退一個格。
我們不妨假設，小蟲所處的世界的一個片斷是：黑黑黑白白黑白……，小蟲從左端開始。

那麼小蟲讀到這個片斷會怎樣行動呢？它先讀到黑色，然後根據程式前移一個方格，於是就會得到另外一個黑色資訊，這個時候它會根據程式再次前移一個方格，仍然是黑色，再前移，這個時候就讀到白色方格了，根據程式它應該後退一個格，這個時候輸入就是黑色了，前移，白色，後移……，可以預見小蟲會無限的迴圈下去。
然而，現實世界中的小蟲肯定不會這樣傻的在那裡無限迴圈下去。我們還需要改進這個最簡單的模型。首先，我們知道小蟲除了可以機械地在世界上移動以外，還會對世界本身造成影響，因而改變這個世界。比如蟲子看到旁邊有食物，它就會把那個東西吃掉了。在我們這個模型中，也就相當於我們必須假設小蟲可以改寫紙帶上的資訊。

因而，小蟲可能的輸出動作集合就變成了：

O={前移，後移，塗黑，塗白}。

這個時候，我們可以把程式1改為比如：

程式2：

輸入 輸出
黑   前移
白   塗黑

紙帶是黑黑白白黑……，小蟲會怎樣行動呢？下面的圖表示出了這個例子中每一步小蟲的位置（標有圓點的方格就是小蟲的當前位置），以及紙帶的狀況。
開始：小蟲在最左邊的方格，根據程式的第一行，讀入黑色應該前移。

第二步：仍然讀入黑，根據程式的第一行，前移。

第三步：這個時候讀入的是白色，根據程式的第二行，應該把這個方格塗黑，而沒有其他的動作。假設這張圖上方格仍然沒有塗黑，而在下一時刻才把它表示出來。

第四步：當前方格已經是黑色的，因此小蟲讀入黑色方格，前移。

第五步：讀入白色，塗黑方格，原地不動。

第六步：當前的方格已經被塗黑，繼續前移。

第七步：讀入黑色，前移

小蟲的動作還會持續下去……。我們看到，小蟲將會不停的重複上面的動作不斷往前走，並會把所有的紙帶塗黑。
顯然，你還可以設計出其他的程式來，然而無論你的程式怎麼複雜，也無論紙帶子的情況如何，小蟲的行為都會要麼停留在一個方格上，要麼朝一個方向永遠運動下去，或者就是在幾個方格上來回打轉。然而，無論怎樣，小蟲比起真實世界中的蟲子來說，還有一個致命的弱點：那就是如果你給它固定的輸入資訊，它都會給你固定的輸出資訊 ！因為我們知道程式是固死的，因此，每當黑色資訊輸入的時候，無論如何它都僅僅前移一個方格，而不會做出其他的反應。它似乎真的是機械的 ！
如果我們進一步更改小蟲模型，那麼它就會有所改進，至少在給定相同輸入的情況下，小蟲會有不同的輸出情況。這就是加入小蟲的內部狀態！我們可以作這樣的一個比喻：假設黑色方格是食物，蟲子可以吃掉它，而當吃到一個食物後，小蟲子就會感覺到飽了。當讀入的資訊是白色方格的時候，雖然沒有食物但它仍然吃飽了，只有當再次讀入黑色時候它才會感覺到自己饑餓了。因而，我們說小蟲具有兩個內部狀態，並把它內部狀態的集合記為：S={饑餓，吃飽}。這樣小蟲行為的時候就會不僅根據它的輸入資訊，而且也會根據它當前的內部狀態來決定它的輸出動作，並且還要更改它的內部狀態。而它的這一行動仍然要用程式控制，只不過跟上面的程式比起來，現在的程式就更複雜一些了，比如：
程式3：

前移，保持吃飽

輸入	當前內部狀態	輸出	下一時刻的內部狀態
黑	饑餓	塗白	吃飽
黑	吃飽	後移	饑餓
白	饑餓	塗黑	饑餓
白	吃飽	前移	吃飽

這個程式複雜多了，有四行，原因是你不僅需要指定每一種輸入情況下小蟲應該採取的動作，而且還要指定在每種輸入和內部狀態的組合情況下小蟲應該怎樣行動。看看我們的蟲子在讀入黑白白黑白……這樣的紙帶的時候，會怎樣？仍然用下面的一系列圖來表示，灰色的圓點表示饑餓的小蟲，白色的圓點表示它吃飽了。為了清晰，我們把小蟲將要變成的狀態寫到了圖的下一行。
假定它仍然從左端開始，而且開始的時候小蟲處於饑餓狀態。這樣讀入黑色，當前饑餓狀態，根據程式第一行，把方格塗白，並變成吃飽（這相當於把那個食物吃了，注意吃完後，小蟲並沒動）。
塗白方格，變成吃飽

第二步：當前的方格變成了白色，因而讀入白色，而當前的狀態是吃飽狀態，那麼根據程式中的第四條前移，仍然是吃飽狀態；
前移，仍然吃飽

第三步：讀入白色，當前狀態是吃飽，因而會重複第二步的動作。

第四步：仍然重複上次的動作。

前移，保持吃飽
第五步：讀入黑色，當前狀態是吃飽，這時候根據程式的第二行應該後移方格，並轉入饑餓狀態；
後移，變成饑餓

第六步：讀入白色，當前饑餓狀態，根據程式第三行應該塗黑，並保持饑餓狀態（各位注意，這位小蟲似乎自己吐出了食物 ！）；

塗黑，保持饑餓
第七步，讀入黑色，當前饑餓，於是把方格塗白，並轉入吃飽狀態（呵呵，小蟲把剛剛自己吐出來的東西又吃掉了 ！）。

塗白，轉入吃飽

第八步，讀入白色，當前吃飽，於是前移，保持吃保狀態。

前移，保持吃飽

這時候跟第四步的情況完全一樣了，因而小蟲會完全重複5、6、7、8步的動作，並永遠迴圈下去。似乎最後的黑色方格是一個門檻，小蟲無論如何也跨越不過去了。
小蟲的行為比以前的程式複雜了一些。儘管從長期來看，它最後仍然會落入機械的迴圈或者無休止的重複。然而這從本質上已經與前面的程式完全不同了，因為當你輸入給小蟲白色資訊的時候，它的反應是你不能預測的 ！它有可能塗黑方格也有可能前移一個。當然前提是你不能打開小蟲看到它的內部結構，也不能知道它的程式，那麼你所看到的就是一個不能預測的滿地亂爬的小蟲。如果小蟲的內部狀態數再增多呢，那麼它的行為會更加的不可預測 ！
好了，如果你已經徹底搞懂了我們的小蟲是怎麼工作的，那麼你已經明白了圖靈機的工作原理了 ！因為從本質上講，最後的小蟲模型就是一個圖靈機 ！
2、如何理解圖靈機模型*
剛才用小蟲說明瞭圖靈機的工作原理，相信你的第一個反映就是，這樣的模型太簡單了 ！他根本說明不了現實世界中的任何問題 ！下面，我就要試圖說服你，圖靈機這個模型是偉大的 ！
首先，我想說的是，其實我們每一個會決策、會思考的人就可以被抽象的看成一個圖靈機。
為什麼可以做這種抽象呢？首先我們可以考慮擴展剛才說的小蟲模型。因為小蟲模型是以一切都簡化的前提開始的，所以它的確是太太簡單了。然而，我們可以把小蟲的輸入集合、輸出行動集合、內部狀態集合進行擴大，這個模型就一下子實用多了。首先，小蟲完全可以處於一個三維的空間中而不是簡簡單單的紙帶。並且小蟲的視力很好，它一下子能讀到方圓500米的資訊，當然，小蟲也可以擁有其他的感覺器官，比如嗅覺、聽覺等等，而這些改變都僅僅是擴大了輸入集合的維數和範圍，並沒有其他更本質的改變。同樣道理，小蟲可能的輸出集合也是異常的豐富，它不僅僅能移動自己，還可以盡情的改造它所在的自然界。進一步的，小蟲的內部狀態可能非常的多，而且控制它行為的程式可能異常複雜，那麼小蟲會有什麼本事呢？這就很難說了，因為隨著小蟲內部的狀態數的增加，隨著它所處環境的複雜度的增加，我們正在逐漸失去對小蟲行為的預測能力。但是所有這些改變仍然沒有逃出圖靈機的模型：輸入集合、輸出集合、內部狀態、固定的程式！就是這四樣東西抓住了小蟲資訊處理的根本。
我們人能不能也被這樣的抽象呢？顯然，輸入狀態集合就是你所處的環境中能夠看到、聽到、聞到、感覺到的所有一起，可能的輸出集合就是你的每一言每一行，以及你能夠表達出來的所有表情動作。內部狀態集合則要複雜得多。因為我們可以把任意一個神經細胞的狀態組合看作是一個內部狀態，那麼所有可能的神經細胞的狀態組合將是天文數字 ！
似乎你會說，這個模型根本不對，還有很多思維本質的東西沒有概括進去。比如記憶問題，人有記憶，圖靈機有麼？其實，只要圖靈機具有了內部狀態，它就相應的具有了記憶。比如上面講到的具有饑餓和吃飽兩種狀態的小蟲就會記住它所經歷過的世界：如果吃到食物就用吃飽狀態來“記住”吃過的食物。什麼是記憶呢？假如你經歷了一件事情並記住了它，那麼只要你下一次的行動在相同條件下和你記住這件事情之前的行動不一樣了，就說明該事情對你造成了影響，也就說明你確實記住了它。
學習的問題反映在模型中了麼？學習是怎麼回事兒呢？似乎圖靈機模型中不包括學習，因為學習就意味著對程式的改變，而圖靈機是不能在運行過程中改變它的程式的。然而，我們不難假設，你實際上並不能打開一個人的腦袋來看，所以它的實際程式規則你是不知道的。很有可能一個圖靈機的規則沒有改變，只不過啟動了它的某些內部狀態，因而它的行為發生了本質上的變化，儘管給它相同的輸入，它給出了完全不同的輸出，因而在我們看來，它似乎會學習了 ！而實際上，這個圖靈機的程式一點都沒變。
還有很多很多現象似乎都能被圖靈機包括，什麼是人類的情緒、情感？你完全可以把它看作是某種內部狀態，因而處於心情好的情緒下，你的輸入輸出是一套規則，而心情不好的時候則完全是另一套。這仍然沒有逃出圖靈機的模型範圍。
接下來的問題就是我們人的思維究竟是不是和圖靈機一樣遵循固定的程式呢？這個問題似乎初看是不可能的，因為人的行為太不固定了 ！你不可預言它 ！然而我會爭辯道，無論如何神經元傳遞資訊、變化狀態的規律都是固定的，可以被程式化的，那麼作為神經元的整體：腦的運作必然也要遵循固定的規則也就是程式了。那麼，如果是這樣，正如圖靈相信的，人腦也不會超越圖靈機這個模型，所以，人工智慧也必然是可能的 ！然而，我認為針對這個問題的答案很有可能沒有這麼簡單，我們將在最後詳細討論這個問題。
無論如何，我相信你已經能夠體會到了，圖靈機模型實際上是非常強有力的 ！

四、計算

1、什麼是計算
說了這麼多，雖然也許你已經瞭解到了圖靈機的威力，也許還將信將疑，然而，你肯定仍然看不出來圖靈機和計算有什麼關係。而實際上，圖靈機是一個理論電腦模型，它最主
要的能耐還是在於計算上 ！所以，下面我們就來看看什麼是計算 ！
我可以先給出一個很摩登的對計算概念的理解：廣義上講，一個函數變化如把x變成了f(x)就是一個計算 ！如果我們把一切都看作是資訊，那麼更精確的講，計算就是對資訊的變換！如果採用這種觀點，你會發現，其實自然界充滿了計算 ！如果我們把一個小球扔到地上，小球又彈起來了，那麼大地就完成了一次對小球的計算。因為你完全可以把小球的運動都抽象成資訊，它無非是一些比如位置、速度、形狀等等能用資訊描述的東西嘛，而大地把小球彈起來就無非是對小球的這些資訊進行了某種變換，因而大地就完成了一次計算 ！你可以把整個大地看作是一個系統，而扔下去的小球是對這個系統的輸入，那麼彈回來的小球就是該系統的輸出，因而也可以說，計算就是某個系統完成了一次從輸入到輸出的變換！
這樣理解不要緊，你會發現，現實世界到處都是計算了 ！因為我們完全可以把所有的自然界存在的過程都抽象成這樣的輸入輸出系統，所有的大自然存在的變數都看作是資訊，因而計算無處不在 ！也的確，正是採取了這樣的觀點，國外才有可能發明什麼DNA電腦、生物電腦、量子電腦這些新鮮玩藝 ！因為人家把DNA的化學反應、量子世界的波函數變換都看作是計算了，自然就會人為地把這些計算組合起來構成電腦了。然而，似乎我們的理論家們還在力圖證明關於圖靈機的某個定理呢，卻完全沒有意識到計算其實就是這樣簡單！
下麵回到圖靈機 ！為什麼說圖靈機是一個計算的裝置呢？很簡單，圖靈機也是一個會對輸入資訊進行變換給出輸出資訊的系統。比如前面說的小蟲，紙帶上的一個方格一個方格的顏色資訊就是對小蟲的輸入，而小蟲所採取的行動就是它的輸出。不過這麼看，你會發現，似乎小蟲的輸出太簡單了。因為它僅僅就有那麼幾種簡單的輸出動作。然而，不要忘了，複雜性來源於組合 ！雖然每一次小蟲的輸出動作很簡單，然而當把所有這些輸出動作組合在一起，就有可能非常複雜！比如我們可以把初始時刻的紙帶看作是輸入資訊，那麼經過任意長的時間比如說100年後，小蟲通過不斷的塗抹紙帶最後留下的資訊就是輸出資訊了。那麼小蟲完成的過程就是一次計算。事實上，在圖靈機的正規定義中，存在一個所謂的停機狀態，當圖靈機一到停機狀態，我們就認為它計算完畢了，因而不用費勁的等上100年。
2、計算的組合
更有意思的是，我們可以把若干個計算系統進行合併構成更大的計算系統。比如還是那個小球吧，如果往地上放了一個蹺蹺板，這樣小球掉到地上會彈起這個蹺蹺板的另一端，而蹺蹺板的另一邊可能還是一個小球，於是這個彈起的小球又會砸向另一個蹺蹺板……。
我們自然可以通過組合若干圖靈機完成更大更多的計算，如果把一個圖靈機對紙帶資訊變換的結果又輸入給另一台圖靈機，然後再輸入給別的圖靈機……，這就是把計算進行了組合 ！也許你還在為前面說的無限多的內部狀態，無限複雜的程式而苦惱，那麼到現在，你不難明白，實際上我們並不需要寫出無限複雜的程式清單，而僅僅將這些圖靈機組合到一起就可以產生複雜的行為了。
有了圖靈機的組合，我們就能夠從最簡單的圖靈機開始構造複雜的圖靈機。那麼最簡單的圖靈機是什麼呢？我們知道最簡單的資訊就是0和1，而最簡單的計算就是對0或1進行布耳運算。而布耳運算本質上其實就三種：與、或、非。從最簡單的邏輯運算操作最簡單的二進位資訊出發我們其實可以構造任意的圖靈機 ！這點不難理解：任何圖靈機都可以把輸入、輸出資訊進行01的編碼，而任何一個變換也可以最終分解為對01編碼的變換，而對01編碼的所有計算都可分解成前面說的三種運算。也許，現在你明白了為什麼研究電腦的人都要去研究基本的布林電路。奧秘就在於，用布林電路可以組合出任意的圖靈機 ！

3、征服無限的方法 ！
回憶你小時候是如何學會加法運算的。剛開始的時候，你僅僅會死記硬背。比如你記住了1+1=2，記住了2+4=6，……。然而無論你記住多少固定數字的運算，你都不叫學會了加法。原因很簡單，假如你記住了n個數字的加法，那麼我總會拿出第n+1個數字是你沒有記住的，因此你還是不會計算。原則上。自然數的數字是無窮的，所以任何兩個數的加法可能結果也是無窮的，而如果採用死記硬背的方法，我們頭腦怎麼可能記住無窮數字的計演算法則呢？但是隨著年齡的增長，你畢竟還是最終學會了加法運算 ！說來奇怪，你肯定明白其實加法運算並不需要記住所有數位元的運算結果，而僅僅需要記住10以內的任意兩個數的和，並且懂得了進位法則就可以了。
你是怎麼做到的呢？假設要計算32+69的加法結果，你會把32寫到一行，把69寫到下一行，然後把他們對齊。於是你開始計算2+9=11，進一位，然後計算3+6=9，再計算9+1=10再進一位，最後，再把計算的這些每一位的結果都拼起來就是最終的答案101。這個簡單例子給我們的啟發就是：作加法的過程就是一個機械的計算過程，這裡輸入就是32和69這兩個數位，輸出就是101。而你的程式規則就是具體的把任意兩個10以內的數求和。這樣，根據固定的加法運算程式你可以計算任兩個數的加法了。

因而，可以說計算這個簡單的概念，是一種用有限來應付無限的方法 ！我們再看一個例子：假如給你4組數字：1,2 3,6 5,10 18,36，就這4組，這時問你102對應的數是多少？很顯然，如果僅僅根據你掌握的已知數對的知識，是不可能知道答案的，因為你的知識庫裡面沒有存放著102對應數字的知識。然而，如果你掌握了產生這組數對的程式法則，也就是看到如果第一個數是x，那麼第二個數就是2x的話，你肯定一下子就算出102對應的是204了。也就是說，你實際上運用2x這兩個字元就記住了無限的諸如1,2 3,6 102,204所有這樣的數對。
這看起來似乎很奇怪。我怎麼可能運用有限的字元來應對無限種可能呢？實際上，當沒有人問你問題的時候，你存儲的2x什麼也沒有，而當我問你102對應的是多少？我就相當於給你輸入了資訊：102，而你僅僅是根據這個輸入資訊102進行一系列的加工變換得到了輸出資訊204。因而輸入資訊就好比是原材料，而你的程式規則就是加工的方法，只有在原材料上進行加工，你才能輸出最終產品。
這讓我不禁想起了專家系統方法。其實專家系統就是一個大的規則庫。也就相當於存儲了很多很多的1,2 3,6 5,10這樣特殊的規則對。而無論它存儲的東西再多，總歸會是有限的，你只要找到一個它沒有存儲到的問題，它就無能為力了。因而專家系統就會在你問到102對應是多少的時候失敗 ！如何解決問題？人們想出了很多方法，就比如元規則的方法，其實元規則就相當於剛才所說的計算加法的程式，或者2x這樣的東西。運用元規則的確可以應對無限種情況了。所以，這就是為什麼你問電腦任何兩個數相加是多少，它總能給出你正確的答案的原因，雖然它不必記住所有這些加法對的資訊。
然而僅僅是元規則就能解決所有問題麼？
假如給你三組數對，排列成一個表：
1,2 3,6 4,8   100,200
3,9 2,6 8,24 100,300
1,4 2,8 3,12 100,400
…………………….
…………………….
…………………….
…………………….
…………………….
那麼問你在第6行上，3這個數字對應的是多少？我們先要找出第一行的規律是2x沒有疑問，第二行呢？是3x，第三行是4x，那麼第6行就應該是7x了，因而在第6行上3
應該對應的是21了 ！這裡跟前面不太一樣的是，雖然我們得到了每一行的規則比如第一行的2x，但是隨著行數的增加，這個規則本身也變化了，因而第2行是3x，第3行是4x等等，因而我們又得到了一個規則本身的規則，即如果行數是n的話，那麼這一行的規則就是(n+1)x。我們顯然能夠根據輸入的n和x計算出數值。把這個道理放到專家系統裡面，這種原理就是元規則的規則，元規則的元規則……，應該是無窮的 ！然而專家系統本身並不會自動的歸納這些規則，人必須事先把這些元規則寫到程式裡，這也就是專家系統最大的弊端。而我們人似乎總能在一些個別的事件中歸納出規則。進一步問，機器可以歸納麼？這就相當於說：可以為歸納方法編出程式麼？這也是一個很有趣的問題，下面將要詳細討論 ！可以設想，假如我們找到了真正歸納的方法，那麼編寫出這樣的程式，它就會一勞永逸的自己進行學習歸納了。我們完全再也不用給他編制程式和規則了。這正是人工智慧的終極目標 ！
4、歸納*
記得金庸在他的武俠小說《倚天屠龍記》中曾講述了這樣一段故事：武林泰斗張三豐在情急之下要把他新創的武功“太極拳”傳授給張無忌。張無忌除了有一身精湛的“內功修為”以外還對武學具有極高的悟性。因而當張三豐給他打過一趟太極拳以後，他就把所有的招式全部記下來了並且當場把所學的太極拳重新再打給張三豐看。在張無忌練拳的過程中，張三豐反復問他一個問題：“你還記得多少？”。他的回答令其他人異常不解，因為他越在那裡揣摩太極拳的奧秘，忘記的招數也越來越多。旁邊的人不明白，這樣的學法忘的這麼快，怎麼可能學會武功呢？然而，沒過多長時間，張無忌說已經忘掉了所有的招式。張三豐笑著說：“不錯，你終於學會了‘太極拳’”。
從這個例子中，我們看到了什麼？張無忌之所以能學會太極拳，正是因為他已經能夠從具體的一招一式之中抽象出了更高一層次的武學規律，因而，當他把所有的有形的武功招數都忘記的時候，已經掌握了太極拳的精髓。而太極武功講究的就是借力打力，以柔克剛。說白了就是事先並沒有固定招式存在，而等到敵人向我進攻的時候我再動態的生成破解的招術。
用到圖靈機模型中，我們不難發現，如果把具體的武功招術比喻成一些輸入，而對應招術比喻成圖靈機的輸出，那麼太極所講究的借力打力、以柔克剛的方法其實就是類似上節講過的2x這樣的圖靈程式！因而張無忌學太極拳的過程就是從特殊的輸入輸出提升到了一般的演算法的過程。也可以說，張無忌運用了歸納學習法 ！
然而，仔細觀察上一節的敘述，我們會發現。雖然圖靈機能夠將2x這樣的法則計算得出結果，但是抽象出2x本身並不是機器自動產生的，而是需要我們外在的人程式設計進去。那麼，面對這樣的問題，究竟圖靈機能不能像張無忌一樣進行歸納思維呢？
可以設想，如果電腦真有了張無忌那兩下子，我們人類可要省事兒多了 ！我們甚至不需要為電腦編程式，它就會自動的從若干個具體事例中歸納出一般的通用規律來。然而，究竟電腦能不能具有真正的歸納能力呢？讓我們來仔細考慮一下這個問題。
我們說如果電腦能自動歸納，也就意味著我們可以為歸納方法來編寫一段程式P。這個程式可以理解為輸入的是一些特殊的數對，輸出的是能夠生成這些數對的程式。也就是說輸入具體的“招術”，輸出的是這些“招術”的一般規律。如果說程式P真正可以自己歸納，那麼P就必然可以歸納出所有的規律。我們已經討論過了，其實任何一個程式都能夠被看作是對輸入的一個變換而得到輸出。那麼程式P自然也是。假設這些對子(a,b),(c,d),(e,f),……都是程式P的輸入輸出對，那麼我們挑選出前1000個（總而言之是足夠多的對子）。把這1000個特殊情況輸入到P中，那麼P就應該能夠產生這些對子的共性，也就是P自己這個程式了 ！換句話說，程式P產生了它自己，P自己把自己給歸納出來了 ！這似乎陷入了迴圈之中 ！另外，我們人類設計出來P，如果P可以歸納所有的規律，那麼P能否也能歸納出“人歸納P”本身這個規律呢？仍然是怪圈問題 ！這樣的問題似乎還有很多。反過來講，如果假設歸納出所有規律的程式P不存在，那麼為什麼我們人類總能歸納出規律呢？什麼樣的具體問題是可歸納的，什麼問題是不可歸納的？然而這些看起來非常重要的問題在目前還沒有統一的答案 ！
我們還將會看到很多問題都涉及到邏輯中的迴圈，而由於計算理論已經觸及了邏輯、資訊的根本，所以把一些問題引向邏輯迴圈並不奇怪。

五、模擬

1、什麼是模擬？
什麼是模擬？又是一個基本的問題，愛因斯坦說過，越是基本的概念就越是難以刻畫清楚。模擬這個概念就是一個很難說清的問題。
如果你站在一個朋友面前，衝著他做了一些鬼臉。那麼他也會學著你的動作衝你做鬼臉，那麼他就對你進行了模擬。
很明顯，在你和你朋友之間存在著一系列的對應關係：你的手對應他的手，你的眼睛對應它的眼睛，你的嘴巴對應他的嘴巴……。而且你的手、眼睛、嘴巴做出來的動作也會對應他的手、眼睛、嘴巴做出來的動作。因而，模擬的關鍵是對應！如果集合A中的元素可以完全對應B中的元素，那麼A就可以模擬B。
仍然用你對你的朋友做鬼臉的例子，假如這次你做出的鬼臉以及動作沒有被他立即模仿而是被他用某種符號語言記錄到了日記本上了。比如：“X年X月X日，瘋子XX對我做了一個鬼臉：他伸出了左手食指放到了右眼下面往下拉他臉上的肉，並且吐出了他長長的舌頭 ！”。過了N多天后，你的這位朋友掏出了日記本，按照上面的描述對著大家做了這個鬼臉。很顯然他仍然模擬了你當時的動作。那麼，你朋友日記本上的那段對話描述是不是對你鬼臉動作的模擬呢？似乎答案是否，因為這段文字跟你沒有半點相像。然而你的朋友正是根據這段描述才做出了對鬼臉動作的模擬。也就是說，他把那段文字翻譯成了他的動作，而他這個動作就是對你的模擬。這個翻譯的過程很顯然就是某種資訊的變換，我們完全可以把它理解為一個計算的過程，也就是可以用圖靈機來實現的演算法過程。所以，我們說日記本上的那段指令也構成了對你鬼臉動作的類比，原因是這些資訊也與你的鬼臉動作構成了對應。具體的，我們可以用下面的圖表示：
A <–模擬 B<–變換–> C
這裡A是你的鬼臉動作，B是你朋友做出來的鬼臉動作，C是日記本上的描述。你朋友的動作B模擬了你的動作A，而B的動作資訊是通過執行C上的描述得到的，也就是說存在著一個從C到B上資訊的變換。這樣我們認為C也對A進行了模擬。
2、圖靈機之間的模擬
下麵來考慮圖靈機之間的模擬。按照前面的定義，一台圖靈機包括：輸入集合I，輸出集合O，內部狀態集合S，程式規則表T四個要素。那麼如果兩個圖靈機之間的這些元素都存在剛才說的對應關係，就認為兩個圖靈機可以相互模擬了。然而圖靈機的功能是完成對輸入資訊進行變換得到輸出資訊的計算。我們關心的也僅僅是輸入輸出之間的對應關係。因而一台圖靈機A如果要模擬B並不一定要模擬B中的所有輸入、輸出、內部狀態、程式規則
表這些元素，而只要在給定輸入資訊的時候，能夠類比B的輸出資訊就可以了。
因此，我們可以用下面的圖來表示圖靈機之間的模擬：

也就是說在給定相同輸入資訊的情況下，只要輸出資訊o’能夠類比資訊o就可以，也就認為B模擬了A。而資訊o’對資訊o的類比又符合我們上面對一般集合之間類比的定義。也就是說如果存在另外一台圖靈機能夠把資訊o’計算並對到資訊o，就認為o’模擬了o。說白了也就是o’可以與o不一樣，但是只要你能用一個圖靈機把o’經過一系列運算變換到相同的o，就認為o’模擬了o。因而也就是圖靈機B模擬了圖靈機A。
進一步，我們可以假設A和B輸入的資訊也不一樣，一個i，另一個是i’，那麼如果i和i’之間也存在著模擬對應關係的話，我們仍然認為B可以模擬A。也就是下面的圖：

有一點需要注意，如果A圖靈機模擬了B圖靈機，那麼並不一定B圖靈機可以模擬A圖靈機。因為有可能A圖靈機比B圖靈機處理的資訊更多。也就是說假如B能處理的資訊就是1，2，34，而A處理的資訊除了這四個數之外，還有5,6,7,8，那麼顯然當輸入1234的時候A能夠模擬B，而當輸入5678的時候B沒定義了，不能完成任何操作。在這個時候B顯然不能模擬A了。
3、計算等價性
講了這麼多關於模擬的知識有什麼用呢？模擬的一個關鍵作用就是闡明什麼是等價的。比如為了完成加法運算，你寫了一段程式，而我也寫了另一段程式，雖然我們兩個的程式可能完全不一樣，然而只要我們兩個程式之間能夠相互類比，也就是說只要給定兩個數，我們都能正確的一模一樣的算出它們的和，那麼我們兩個程式就是等價的 ！
具體地說，如果A能夠模擬B，並且B也能模擬A，那麼A和B就是計算等價的。計算等價性是非常強有力的，因為它揭示了在我們這個宇宙中某種非常普遍的規律。我們仍然用剛才說的加法演算法為例子來說明。雖然計算兩個數的加法的方法可能有無窮多種，也有可能用各種各樣的電腦語言，什麼C,Basic,JAVA等等來實現，更有可能用在不同的電腦上，然而所有這些程式，這些計算的結果意義都是相同的。也就是說所有與加法運算演算法計算等價的電腦程式都是一回事兒，因而加法演算法這個東西是某種永恆而獨立的 ！
看 ！我們在宇宙中找到了某種永恆性了，這種永恆性反映了宇宙規律中某種本質上的美！計算等價性就和能量守恆定律一樣具有這種高級的對稱性，我甚至覺得計算等價性要比能量守恆定律更加深刻 ！因為無論如何能量守恆定律仍然是刻畫了物理系統的某種屬性，而計算等價性則刻畫的是非常廣泛的資訊系統之間的某種守恆和對稱性，而一切系統都可以被轉換為資訊系統，甚至是物質世界，所以，計算等價性是跨越所有系統之間的某種高級對稱的、永恆的、美的東西。
為了進一步理解計算等價性的威力所在，我們不妨想像一下。假設我們能夠用電腦類比某個人，比如說張三的思維程式化了（也就是假設真正的人工智慧可以實現了）。那也就是說我們可以用一個電腦軟體X來完成對張三思維的模擬。這樣，這個軟體就會在一切與它具有計算等價性的程式甚至系統上實現張三這個人的思維過程 ！比如我們完全有可能讓一大堆分子的碰撞來實現X這個軟體，那麼就會在這大堆分子碰撞的過程中完成對張三思維的模擬，也就是說張三這個人的意志蹦到了這一大堆分子系統中去了 ！更進一步，我們還可以找來足夠多的人比如這個星球上所有的人來模擬那大堆分子的碰撞，從而完成軟體X的計算。這意味著什麼？意味著張三這個人的思維或者說意識在那群人的整體上突現了 ！很有可能，這些構成軟體X的人都並沒有意識到在他們上層的張三的意識的出現。更有趣的是，張三自己很有可能就在那一群人之中呢 ！
相信你已經能夠參悟到了什麼是計算等價性的威力了，那麼我也相信你能夠理解為什麼說任何一台我們使用的電腦都不過是圖靈機的翻版了。

4、意義*

考慮下麵三句話：“請把窗戶關上 ！”，“Please close the window ！”，“01001110111”。這三句話分別說給不同房間中的三個人。第一句話告訴給一個中國人，於是他關上了窗戶；第二句話告訴了一個外國人，他也關上了窗戶；第三句話告訴的是一個機器人，他也關上了窗戶。這三句話從表面看顯然是完全不一樣的，然而當它們讓不同的對象來聽的時候，卻達到了相同的最終結果：窗戶被關上了。那麼，我們自然會想，這三句話有何相同呢？顯然，答案是他們的意義相同。然而什麼又是意義呢？
真正回答意義的本質是一個很困難的問題，現在人們正在努力理解語義是什麼。雖然我們仍沒有完全回答這個問題，但是，不妨從圖靈機、計算以及計算等價性的觀點來考慮該問題。如果把中國人、外國人、機器人看作是三個不同的圖靈機，而這三句話就是對這個三個圖靈機的輸入資訊，那麼最終的結果就是圖靈機計算的輸出。這個時候我們看到三種結果是相同的。也就是說這些圖靈機之間是可以相互模擬的。
考慮這三句話，顯然它們都具有相同的意義。而根據前面的敘述，能夠相互模擬的圖靈機是具有相同的計算等價性的。因而描述聽到關窗指令後並按照指令行事的圖靈機具有相同的計算等價性。而這種計算等價性就好像是前面說到的加法規則一樣是獨立於計算系統、執行機構的。因而，我們能得到下面的圖：

通過這個對比圖，我們不難得出結論：所謂語言的意義，就是執行這個語言系統的計算等價性 ！
我們如何知道不同的語言表達了相同的意義呢？顯然，我們只要有了翻譯就可以明白“請把窗戶關上”與“Please close the window”具有相同意義，而翻譯所作的工作無非就是輸入中文資訊輸出英文資訊這樣的資訊轉換工作，因而，也就是一個計算過程 ！
然而當不存在從一個語言到另外一個語言的翻譯的時候，我們也並不能斷定某一個符號序列對於固定的圖靈機是否有意義。這就是說，我們雖然不能明白鳥叫是什麼含義，但並不能否認它們的叫聲可能有意義，因為只有鳥自己才能明白叫聲的含義。

六、萬能圖靈機

1、編碼
其實我說的這個“萬能圖靈機”就是電腦術語中的“通用圖靈機”，英文是Universal Turing Machine。而我之所以稱之為“萬能圖靈機”完全是因為這個名字似乎聽起來更加貼切。
前面已經講述了模擬的概念，那麼自然會產生這樣一個問題：是否存在一台圖靈機能夠模擬所有其他的圖靈機呢？答案是存在的。這種能夠模擬其他所有圖靈機的圖靈機就叫做通用圖靈機，也就是我們所說的“萬能圖靈機”。這種機器在圖靈計算這個範疇內，是萬能的 ！
“萬能圖靈機”會怎樣工作呢？假如我把資訊X輸入到了圖靈機M中，M就能計算出一個結果O。那麼如果我把X和M的資訊都輸入給萬能圖靈機，那麼萬能圖靈機也會輸出O，也就是萬能圖靈機可以模擬任何一台特殊的圖靈機。這樣的話我們就可僅僅通過改變輸入X和M的值就能“改變”萬能圖靈機的程式規則了。因而也可以認為萬能圖靈機就是可以任意程式設計的。這裡的改變兩個字加上了引號，是因為事實上任何圖靈機在誕生之後規則就不能改變了，因而我們能夠改變“萬能圖靈機”的規則，僅僅是因為看上去是這樣的，其實根本沒有改變。
要說明為什麼“萬能圖靈機”是存在的，以及它是怎樣模擬其他任何圖靈機的動作的，我們必須先要理解究竟怎樣把任何一台圖靈機輸入到“萬能圖靈機中”，這就需要理解編碼的概念。什麼是編碼呢？你可以理解為對某一堆事物進行“編號”就是編碼。
其實我們每人每天都在跟編碼打交道。每個人都有一個身份證，而這個身份證都有一個ID號碼吧？那麼這個號碼就是你的編號。上學的時候老師給我們每個人都分配一個學號也是編碼。
26個字母能夠被編碼，比如a對應1，b對應2，……，這是顯而易見的。然而任意一個英文單字都是可以被編碼的則不那麼容易一眼看出來。事實上，我們可以按照字典順序把所有的單字都列出來。也就是說字母順序越靠前，字元長度越短的單字排在前面，其次長單字，字母順序靠後的單字就排在後面。比如一種可能的字典順序：
a, about, an…, bad, be, behave…..
只要這樣一排好序，我們就能給每個單字賦予一個數字，最簡單的方法是，給第一個字母分配1，第二個分配2，……，因而我們就給所有的單字都編碼了。
下麵討論任意一個圖靈機能不能被編碼。我們假設討論的所有圖靈機的輸入集合都是僅有0,1兩種，而它的輸出也僅僅有0,1,2,3四個動作分別表示前移，後移，塗寫0，塗寫1。而內部狀態數最多為10000個（總之足夠多就可以了）。下麵考慮程式。
假設圖靈機的程式表為：

當前內部狀態s	輸入數值i	輸出動作o	下一時刻的內部狀態s’
2	1	0	3
1	0	3	2
3	0	1	1
…	…	…	…

那麼我們可以把它寫到一行中，這就是2,1,0,3; 1,0,3,2; 3,0,1,1，注意用“，”分開了內部狀態，輸入數值，輸出動作和下一時刻的狀態，而用“；”分開了一行一行具體的程式。這樣無論這個表有多大，我們都可以把它寫成這樣的一個字串。這個字串就相當於一個英文單詞，這就是對該圖靈機程式的一個描述。同理，其他的圖靈機也能夠得到這樣的一個單詞描述，那麼我們再用字典序的方法對這些描述進行編碼，也就得到了對所有圖靈機的編碼。
如果一台圖靈機的編碼是M，它讀入的資訊是x，這樣只要把M和x用“.”號隔開的方法分開作為資料登錄到“萬能圖靈機”中，運用特殊的演算法，這個萬能的機器就能得出對M計算x的模擬結果了。事實上可以由定理證明萬能圖靈機對於任意的編碼都是存在的，在這裡我們就不敘述證明過程了。
2、自食其尾
既然“萬能圖靈機”能夠模擬任何一台圖靈機的動作，那麼它能不能模擬它自己的動作呢？答案是肯定的。我們首先看到“萬能圖靈機”也是圖靈機，也有固定的輸入、輸出、狀態的集合、固定的程式，因而它也能被編碼。於是我們就可以把它自己的編碼資訊輸入給它自己了。這就好像一條蛇咬到了自己的尾巴。會發生什麼呢？自食其尾就會產生迴圈，雖然我們現在還沒有看到任何不好的徵兆，然而在下一節裡面，我們將看到這種迴圈會誕生什麼樣的結論。而且我們也會看到，其實這個迴圈是和康得對角線法則、哥德爾定理有關的。
圖靈機一旦能夠把程式作為資料來讀寫，就會誕生很多有趣的情況。首先，存在某種圖靈機可以完成自我複製 ！事實上，電腦病毒就是這樣幹的！

我們簡單說明一下，這個特殊的圖靈機是如何構造的。我們假定，如果一台圖靈機X，那麼它的編碼就記為<X>，這樣能夠自我複製的圖靈機T的功能是，把T的編碼<T>寫到紙帶上輸入到“萬能圖靈機”，
那麼“萬能圖靈機”就能根據讀入的<T>，在紙帶上再次輸出<T>的一份拷貝<T>’，並且<T>=<T>’。下面就來大概解釋如何構造這樣的T。首先T由兩部分構成AB。第一部分A的功能是指導“萬能圖靈機”把B的編碼<B>原封不動的列印到紙帶上，這個時候紙帶上就有了<B>，如果這個時候你想用同樣的方法列印<A>到紙帶上是不行的，因為會出現迴圈定義。然而B可以這樣做，讀入紙帶上的資訊X，生成能夠列印X的圖靈機：p(X)的編碼<p(X)>列印到紙帶上，並把X和<p(X)>的內容前後調換，有定理保證這樣的圖靈機是存在的。這樣當B讀到紙帶上的資訊<B>之後就會列印出能夠列印<B>的圖靈機的編碼也就是<A>了，然後把<A>和<B>位置對換就構成了<AB>也就是<P>，所以P把自己進行了一次拷貝。初看起來，這種自我複製的程式是不可能的，因為這包含了無窮無盡的迴圈。P要能產生它自己<P>就意味著P中至少包含了一個<P>，而這個<P>中又包含了至少一個<P>……，最後P必然是一個無限大的程式，然而我們卻能夠證明P是可能的。
有了“萬能圖靈機”還能得到很多有趣的結論，比如假設有一大群圖靈機，讓它們彼此之間隨機的相互碰撞，當碰到一塊的時候，一個圖靈機可以讀入另一個圖靈機的編碼，並且修改這台圖靈機的編碼。那麼這樣一個圖靈機“湯”中會產生什麼呢？聖塔菲研究院的Fontana她已經研究了這個實驗，並得出了驚人的結論：在這樣的系統中會誕生自我繁殖的、自我維護的類似生命的複雜組織，而且這些組織能進一步聯合起來構成更大的組織 ！

七、停機問題

1、無限迴圈
在進行正式討論之前，我們先來看看一個非常簡單的猜硬幣遊戲。
假如我的兩個拳頭中一個藏著一枚硬幣，另一個沒有，然後讓你猜是哪一個？於是你告訴我左手中有。這時候我不會把手張開，而是背過身去做一番手腳，然後把拳頭伸過來，張開手 ！哈，你錯了吧，硬幣在右手中 ！大概傻子都能看出來我的伎倆之所在 ！不用說，採用這種方法我保證百戰百勝。因為我總是等你說出來哪個手有硬幣之後再動態的改變我的策略。所以，這改變之後的狀態就已經不是你猜的了。
大概你會覺得不可思議：其實圖靈停機問題就是一個類似該遊戲的原理 ！
下面我們來看看圖靈停機問題是怎麼回事兒。讓我們考慮這樣一個問題：是否存在一個程式比如說P，能夠判斷出任意一個程式X是否會在輸入Y的情況下陷入無限迴圈？我們不妨設P(X,Y)表示P判斷程式是X，資料是Y的結果。如果存在無限迴圈，那麼P(X,Y)就輸出一個yes。如果不存在無限迴圈，那麼P(X,Y)就輸出一個no。我們的問題就是這樣的P(X,Y)存在麼？這就是停機問題。所謂的某個程式X在輸入Y上停機就是說X不存在著無限迴圈，反過來如果不停機就是存在著無限迴圈，因而這裡停機和無限迴圈是一回事兒。那麼，這種判斷停機問題的程式P存在麼？
答案是不存在的。下面我可以證明我的這個結論。
我們不妨假設程式P存在。那麼我們可以根據P設計一個新的程式Q如下：
X是任何一段程式的編碼：
Program Q(X){
m=P(X,X)
do while (m=no)
…
…
end do
if m=yes
then return
}

這段程式通俗來講就是：輸入任何一段程式X，調用函數P(X,X)並得到返回值m，如果m=no，也就是說根據P的定義，P判斷出程式X作用到它自己身上X不存在無限迴圈。那麼Q就不停的做do while和end do之間的語句。如果m=yes，我們知道這表示P判斷出程式X在X上存在無限迴圈。就返回，結束該函數。
我們可以看到，這樣定義的函數Q(X)是沒有問題的。下麵就進入關鍵時刻了：我們問Q這個程式作用到Q自身的編碼上也就是Q(Q)會不會形成無限迴圈呢？當然我們可以運用強有力的函數P(Q,Q)來計算這個問題。
假設Q(Q)會發生無限迴圈，那麼P(Q,Q)就會返回yes。然而根據Q函數的定義，把X=Q代入其中會發現由於P(Q,Q)返回的是yes，也就是m=yes，因此Q函數會馬上結束，也就是程式Q(Q)沒有發生無限迴圈。然而如果假設Q(Q)不發生無限迴圈，那麼P(Q,Q)應該返回no，這樣根據Q函數的定義，把X=Q代入Q(Q)之中會得到m=no，這樣程式就會進入do while迴圈，而這個迴圈顯然是一個無限迴圈。因而Q(Q)發生了無限迴圈！這又導致了矛盾。
無論Q(Q)會不會發生無限迴圈，都會產生矛盾，然而哪裡錯了呢？答案只能是最開始的前提就錯了，也就是說我們最開始的假設P(X,Y)能夠判斷任意程式X在輸入Y的時候是不可判定的 ！也就是說這樣的程式P(X,Y)不存在 ！
2、如何理解
也許你會感覺整個論證過程有些怪異，為什麼不存在這種P(X,Y)程式呢？而上面的論證過程中僅僅說P(X,Y)當作用到P(Q,Q)上時會產生矛盾。似乎並不能說明P作用到其他程式上不能判斷是否迴圈。比如你可以考慮編寫這樣一段程式，當一發現某個程式do while(T)，這裡T總是為真，這樣的語句出現的時候就判斷這個程式有迴圈。這顯然是可能的。但問題的關鍵是，你假設了P(X,Y)能夠判斷任意的一個程式是否可判定！最關鍵的就是這“任意程式”上了。因為假如你已經按照剛才提到的判斷是否有do while(T)語句的方法寫出了一個程式P來判斷某程式是否可判定，那麼我就會根據你這個程式P再構造出一個程式Q，就是利用上面提到的論證方法，我們不妨寫成QP(這裡下標P的含義表示根據你的程式P而構造的Q)。這樣你的P在遇到了P(Q,Q)這樣的怪東西的時候無能為力了 ！
可能你還不服輸，於是你又改進了你的程式變成了P’，這個時候P’能夠判斷包含了QP這個程式時候的所有程式情況了。那麼我又會根據你的新程式P’來構造出一個更新的QP’，你的程式P’仍然不能判斷，當然你還可以構造P’’,P’’’,……，我也會跟著構造QP’’,QP’’’，……，總而言之這個過程是無窮的 ！因為我總在你之後構造程式，所以你是水我是船，水漲船高，我總能比你高一級吧 ！
這很像剛開始敘述的那個猜硬幣的遊戲。你想猜對我的硬幣，就必須告訴我一個答案是左手還是右手，然而關鍵問題是我總能根據你做出的答案動態調整，使得你永遠也猜不對 ！停機問題也是如此，我總能根據你的程式P來構造你的P判定不出來的問題Q，我總會贏 ！很簡單，因為你總要在我之前構造好P呀，就相當於你總要先說出硬幣在哪個手 ！
3、對角線刪除方法
我在開始的時候就提到了圖靈停機問題、哥德爾定理等等都來源於康得的對角線刪除法則。那麼，下面我們就來看看，如果運用對角線刪除法則如何證明圖靈停機問題呢？採用這種全新的證明思路，也許你會更加清楚地認識到停機問題的本質。
問題沒有變，是否存在一個程式P(X,Y)判斷出來任意一個程式X當輸入Y的時候是否可判定，或者說是否停機。
我們仍然用反證法，假設這樣的P(X,Y)是存在的。而我們知道程式X本身是可以被編碼的。也就是可以為所有的程式進行編號：1，2，3，……，而資料Y本身也是這樣的編號1，2，3，……，因而我們就可以把每一對X和Y的組合排列在一張表上。比如列表示的是資料Y，而行表示的是程式X，這樣P(X,Y)的值也就是yes或no就可以寫在第X行第Y列的對應位置上，表示P(X,Y)判斷出的X作用在輸入Y上是否會可判定的結果。比如下面的表就是一個示例：

	1	2	3	4	5	6	……	Y	……
1	yes	no	no	yes	yes	no	……	yes	……
2	no	no	yes	yes	no	yes	……	no	……
…	……	……	……	……	……	……	……	……	……
X	no	yes	no	yes	yes	no	……	yes	……
…	……	……	……	……	……	……	……	……	……

到這裡沒有發生什麼毛病。我們知道上表中的每一個行都表示一個確定的程式作用到不同的資料上所得到的結果。比如程式X作用在1,2,3,4,……,Y，……上給出的結果就是一個序列：
no,yes,no,yes,yes,no,……,yes，……
下麵我們把上表對角線上的元素挑出來。也就是專門找那些第1行第1列，第2行第2列，……這樣的元素就可以得到一個序列：
yes,no,no,yes, ……
而根據這個序列我們完全可以構造這樣一個反序列：
no,yes,yes,no, ……
也就是說這個序列在每一個位上都與前一個對角線序列不同 ！這個序列就稱為對角線刪除序列。那麼我問，這個對角線刪除序列是否在我們表中的某一行上呢？答案是否定的 ！這個序列必然不在表中。因為它的第一個元素與1,1不同，第二個元素與2,2不同，……。換一種方式解釋就是：假設對角線刪除序列的確在表上列出來了，那麼這個序列必然在某個固定的行，比如說就在第1001行。那麼我們就要考慮這第1001行，第1001列的元素，我們知道根據序列的構造方法(1001,1001)對應表中的元素是與序列上的這個元素不同的 ！因而產生了矛盾，也就是說這個構造出來的對角線刪除序列不在表中 ！
我們完全可以設計出一段程式Q使得Q作用在1,2,……,X,……上產生的序列就是對角線刪除序列，而對角線刪除序列不在表中就意味著程式P並不能判斷出程式Q作用在任意輸入上是否停機。其實這裡的程式Q就是前一節論證停機問題的程式Q。
用對角線刪除的證明方法來看究竟哪裡出錯了呢？錯誤就出在我們假設程式P能夠得到這樣一張完整的表，這張表對所有的程式計算得到是否停機的答案。
4、意味著什麼？
我們已經看到了，的確存在著一類問題我們人類能構造出來，而圖靈機是不能解的，。我們知道圖靈機不能解的問題也就是一切電腦不能解的問題，因而這類問題也叫做不可判定的(Undecidable problem)。因此，必然存在著電腦的極限。實際上，運用我們前面敘述的計算等價性原理，有很多問題都可以被歸結為圖靈停機問題，也就是說圖靈停機問題揭示了宇宙中某種共同性的東西，所有那些電腦不能解決的問題從本質上講都和圖靈停機問題是計算等價的。比如在最開始我們就提到的希爾伯特第10問題就是一個典型的不可判定的問題 ！還有很多問題是不可判定的，尤其是那些涉及到計算所有程式的程式。比如是否存在一個程式能夠檢查所有的電腦程式會不會出錯？這是一個非常實際的問題。我們都知道電腦程式特別容易犯錯，為了檢查出某段電腦程式的錯誤，我們人類就必須對這個程式進行人工的檢查。那麼能不能發明一種聰明的電腦軟體，輸進去任何一段程式碼，這個軟體就會自動幫你檢查輸入的程式是否有錯誤？答案仍然是不存在的，其實這個問題可以被證明和圖靈停機問題實質上是一樣的 ！於是我們的夢想又破滅了 ！
圖靈停機問題也和複雜系統的不可預測性有關。我們總希望能夠預測出複雜系統的運行結果。那麼能不能發明一種聰明的程式，輸入進去某個複雜系統的規則，輸出的是這些規則運行的結果呢？從原則上講，這種事情是不可能的。它也是和圖靈停機問題相同的。因而，我們得出來的結論就是：要想弄清楚某個複雜系統運行的結果，唯一的辦法就是讓這樣的系統實際運作，沒有任何一種電腦演算法能夠事先給出這個系統的運行結果。你很有可能不同意我的觀點，因為畢竟我們能夠發明很多預測複雜系統的方法，而不需要一定要讓複雜系統去真實的運作。但是，你還是沒有理解我這裡的問題。我們強調的是不存在一個通用的程式能夠預測所有複雜系統的運行結果，但並沒有說不存在一個特定的程式能夠預測某個或者某類複雜系統的結果。那麼這種特定的程式怎麼得到呢？顯然需要我們人為地程式設計得到 ！也就是說存在著某些機器做不了的事情，而人能做。這似乎為人工智慧的崇拜者給以了沉重的打擊 ！
人工智慧真的是不可能的麼？彭羅斯曾經寫過一本著名書藉：《皇帝新腦》來論證人工智慧的不可能性。它所運用的方法就是我們上面的邏輯。因為對於任何一個人工智慧程式來說，總存在著它解決不了的問題 ！但是似乎我們人類卻不受這種限制，我們總是能夠發現一個程式是否有迴圈，總是能夠找到對某類複雜系統預測的方法，並且我們還能構造出來圖靈停機問題這樣的問題。然而事實並沒有那麼簡單，反對者馬上就會論證到，其實針對某一個具體的人，比如說就是彭羅斯，我們也能夠運用前面的方法構造出一個彭羅斯自己不能解的問題 ！然而事實情況下要構造彭羅斯不可解的問題太麻煩了，而我們只是說原則上講這種問題是存在的 ！因而電腦超越不了的問題，人自己也超越不了，所以說人工智慧是可能的 ！
看看上面提到的兩方面論證似乎都很有道理，究竟哪個正確呢？真的會存在某個人不可解的類似圖靈停機的問題麼？其實要想徹底回答這個問題就相當於問超越圖靈計算的限制是否可能？如何超越圖靈機停機問題呢？下面我們將詳細討論一下這個問題。
5、超越圖靈計算的限制
我們仍然用那個猜硬幣的遊戲為例來說明。
在進行了幾輪猜硬幣的遊戲之後，你已經很惱火了，認為這樣的遊戲不公平。於是你想了一個妙招來對付我：每當我讓你說硬幣在哪個手中時，你先胡亂的說一個答案，比如左手。這個時候我會根據你的答案動態調整把硬幣放到了右手中。這個時候你趕緊搶著說，不對，我猜你的硬幣在右手 ！我沒辦法只能再次調整策略把硬幣放到了左手。你又趕快說：是在左手 ！……。就是這樣，你也學會了我的方法，根據我的策略不斷調整你的策略從而讓我不可能贏你。能不能把這種方法用到超越圖靈停機問題呢？
前面我們已經看到了類似這樣的過程。假如你寫出了一個程式P能夠判斷所有程式是否停機，那麼我就能夠構造一個程式Q是你的程式判斷不了的。這個時候還沒有結束，你又根據我的Q構造了新的程式P’，然而我又能構造一個程式Q’仍然讓你的程式P’解決不了。但是你沒有結束，又構造了新的程式P’’，我於是又構造Q’’……。
乍一看，似乎這個過程並不能說明任何問題。原因很簡單，我要求的是構造一個固定的程式P判斷出所有程式是否停機。而你給我的並不是一個具體的實實在在的程式，而是一個不斷變化的捉摸不定而虛無飄渺的程式序列 ！並且你的這些總在變化的程式序列總是要根據我構造的程式才會確定改變！
首先一點值得肯定的是，運用這種方法，我們的確能夠超越圖靈計算了，只要反復不停的變換我們的程式就不可能找出它不能解的問題。然而，另一方面又會讓我們很失望：這樣的變換過程並不能給出一個實實在在的程式來 ！我們擁有的僅僅是不斷改變的程式序列，而不是一個實際存在的程式！
這正是問題的關鍵所在：要想徹底超越圖靈計算的限制，我們必須要放棄程式的實在性。也就是說程式在每時每刻都要變化成不是它自己了。那麼這樣的一個不斷變化得不是它自己的怪東西存在麼？
幾千年的人類科學一直在研究實實在在的東西。無論是原子、分子還是電腦程式，它們必須是一個實實在在存在的個體，在這種前提下科學才能夠對它進行研究 ！如果當我們研究它的時候，它已經變得不是它自己了，那麼科學就對它無能為力了。然而，我不禁要提出這樣的問題：真的一切都是固定不變的存在著麼，有沒有某種東西在每一時刻都在變得不是它自己了呢？
這個問題似乎是一個古老的哲學問題，有人問赫拉克利特斯(Heraclitus)，過去的事情能否更改？赫拉克利特回答說：「人不能兩次踏進同一條河流。」(No man ever steps in the same river twice, for it’s not the same river and he’s not the same man.)。我想他說的正是這樣的問題：因為河流在每時每刻都不再是它自己了。河流是一大群流動的水滴構成的整體，在每時每刻這些水滴都在不停的運動、流逝，因而當你兩次踏入這條河的時候，所有的水滴可能都不一樣了，那麼我們怎麼能說這些水滴構成的整體還是同一條河呢？
再考慮我們人自己。你很可能拿著一個你3歲時候的照片興奮的對你的朋友說：“看，我3歲的時候多可愛呀 ！”。然而你這句話意味著什麼呢？意味著照片反映的3歲時的你和現在的你是同一個個體 ！然而，3歲的你和現在的你是多麼不同呀 ！我們知道，你無疑就是一大堆細胞構成的一個整體。而基本生理學知識告訴我們，實際上人體的所有細胞每隔大約8至10年就會因為新陳代謝的作用全部更新一遍。也就是說，你的細胞全被掉了包了，更何況3歲時候的你和現在的你差了多少個8年呀？那憑什麼說那個3歲時候的你就是現在的你呢？
這個問題看似玄學，不過我認為現在我們的確應該認真對待該問題了。儘管從分析的角度來說3歲的你和現在的你的確不是一個個體，然而常識告訴我們，這兩個你的確都是同一個人 ！那就意味著，你這個個體並不是一成不變的一些固定的細胞，而是一個每時每刻都在變化，都在更新的一個一大堆細胞組成的構形。這個構形在每時每刻都要利用更新的一大堆細胞去維持自己的存在 ！我們得到了什麼？和我們前面敘述的超越圖靈機的討論結合起來，就發現，原來人還有赫拉克裡特斯的河流這種東西剛好就滿足那種超越圖靈計算的要求。也就是說人還有赫拉克裡特斯的河流在每時每刻都在不停的更新它自己從而變得不是它自己了。那麼很有可能，某一種做類似變化的個體的變化規律就是不停超越它自己的圖靈停機程式，這樣的虛幻的個體就真的能夠超越圖靈計算了 ！
總結前面的討論，我們不難給出結論，一個固死的能夠被寫出就不再變化的程式不可能超越圖靈計算的限制，然而如果一個程式每時每刻都已經變化得不是它自己了，這個程式就能夠超越圖靈計算。聯繫到人這個個體，我們能得到：因為每時每刻的人都已經由於細胞的變化而變得不再是它自己了，所以人是超越圖靈計算的 ！還記得在前面我提到的一個問題麼：“人腦的資訊處理過程能不能被表示成固定的程式呢？”。我這裡的答案就是否定的 ！也就是說人腦資訊處理的過程並不是一個固定的程式！如何製造真正的人工智慧呢？很顯然，我們不能用一個簡單的程式來構造，而必須是利用其它的方法，這個方法是什麼呢？現在還沒有結論 ！

八、懸而未決

到此，我已經把全部的有關圖靈機、可計算理論、停機問題的一些重要概念介紹完了。然而在計算理論這個領域裡還有很多重要的問題沒有介紹。但我不得不根據我的興趣進行取捨。在整個介紹過程中，一方面我介紹了人們已經得出來的結論，另一方面，我儘量把一些沒有解決的問題展現給大家。回憶起來，這些懸而未決的問題包括下面幾個：

1、是否一切的資訊處理過程都具有固定的程式呢？人腦有固定的程式麼？
2、如何用電腦程式進行歸納？能對所有事物進行一勞永逸的歸納演算法是否存在？
3、什麼是意義，更確切的，什麼是語義？
4、圖靈停機問題是不可超越的麼？
5、人工智慧是否可能？
還有很多問題在本文中沒有提出來，然而我認為也是相當重要的。例如，我們如何用圖靈機模型表示一般的學習過程？若干小的圖靈機是如何自動的構造出更大的圖靈機的（這就是萬事萬物自組織的過程）？生命的目的性如何用圖靈機模型表示？
另外最近的計算主義已經把宇宙中一切的過程都歸結為計算過程了，也就是說到處都是圖靈機正在做運算。那麼我們能不能從圖靈機的角度探討時間和空間的本質呢？我們知道，計算理論另外一大類問題就是探討計算的時間和空間複雜度，那麼這種計算的時間和空間與我們這個宇宙的時間和空間有什麼關係呢？
我希望，在新世紀的今天，人們會最終解答這些問題 ！
書寫之中不免有很多錯誤和遺漏，也可能會存在很多邏輯漏洞，如果你發現了，請告訴我，我不禁感謝 ！
2004年7月

作者：張江(email: jakezj@163.com) 集智俱樂部：http://www.swarmagents.com

---

國外的劉子千

發表： 2011 年 07 月 20 日 | 作者： vinsonhsieh | Filed under: 科技趣味創意 | 發表留言

不懂的人可以先看一下…

http://www.youtube.com/watch?v=9TWGMFTORfw&feature=related

原版Friday.

http://www.youtube.com/watch?v=fuIydD_9jBE&feature=player_embedded#at=190

Kuso版的Saturday ，但Quality很不錯喔

http://www.youtube.com/watch?v=4Chelwb6Ja4&NR=1&feature=fvwp

新聞提到的Autotune…

http://halcyonstar.blogs.com/middle_of_nowhere/2010/09/%E8%87%AA%E5%8B%95%E8%AA%BF%E9%9F%B3autotune.html

I am T-Pain

http://iamtpain.smule.com/

把對話轉歌曲…

http://www.youtube.com/watch?v=OOVeDwFLd1Q&feature=player_embedded#at=87

Another example.

http://www.youtube.com/watch?v=hMtZfW2z9dw&feature=player_embedded

Original source.

http://www.youtube.com/watch?feature=player_embedded&v=_oYnM9DD04s#at=100

---

傅立葉轉換演算法專論

發表： 2011 年 07 月 18 日 | 作者： vinsonhsieh | Filed under: 演算法 | 發表留言

傅立葉轉換演算法專論 – for dragon    2004/11/9 下午 09:52:27

第一節　傅立葉轉換的定義

什麼是傅立葉轉換（Fourier Transform）？在早期的定義中，傅立葉序列（Fourier Series）被看成是一堆波的合成波，其公式約略如下（假設週期為T）：

∞
F(u) = Σ　 (a(k)*cos(ku*2π/T) + b(k)*sin(ku*2π/T))
k=-∞

由於傅立葉序列許多的合成波高低頻特性，在訊號處理（Signal Processing）上有很多應用價值，因此傅立葉序列所在的領域一般又稱為頻譜領域（Frequency Domain），而傅立葉轉換前的領域便稱為時序領域（Time Domain）。然而有很多函數並不屬於週期函數，但我們仍然可以將之寫成類似上面的公式。考慮任何一個函數f(x)，令：

1　 ∞
A(u) = ─ ∫　 f(x) cos(ux) dx
π　-∞

1　 ∞
B(u) = ─ ∫　 f(x) sin(ux) dx
π　-∞

如此便可將f(x)改寫下面類似傅立葉序列的公式，稱為傅立葉積分（Fourier Integral）：

∞
f(x) = ∫　 (A(u)*cos(ux) + B(u)*sin(ux)) du
-∞

再將之化解成複數形態（Complex Form），可得：

1　　 ∞　　 ∞　　　 iuy　　　 -ixu
f(x) = ── ∫　〔 ∫　 f(y) e　 dy 〕 e　　 du
2π　 -∞　　-∞

然後即可定義出傅立葉轉換公式及反轉公式如下：

∞　　　iux
F(u) = ∫　f(x) e　 dx
-∞

1　　 ∞　　　-ixu
f(x) = ── ∫　F(u) e　　du
2π　 -∞

有些書則將傅立葉轉換公式定義為：

1　　　　∞　　　iux
F(u) = ───── ∫　f(x) e　　dx
sqrt(2π)　 -∞

1　　　　∞　　　-ixu
f(x) = ───── ∫　F(u) e　　du
sqrt(2π)　 -∞

或

∞　　　i2πux
F(u) = ∫　f(x) e　　dx
-∞

∞　　　-i2πxu
f(x) = ∫　F(u) e　　du
-∞

不過在本文中，將都以第一式為準。必須注意的是，上面的i=sqrt(-1)，表示
複數（Complex Number）值，而非變數。而式中e^(iux)經尤拉（Euler）公式
化解可得：

iux
e　　= cos(ux) + i sin(ux)

當然，這類數學式子的化解較為繁複，因此本文中並不將之全數列出。有興趣的人不妨自行導導看，或是參考其他書籍了解。

由於電腦內的資料都是離散（Discrete）的形式，因此上面這種連續性積分的公式，並不太適用於電腦資料的運算上，於是便有離散式傅立葉轉換（簡稱DFT）公式的需要。至於如何導出這個公式呢？假設給定一個長度為N的向量A，欲轉換成傅立葉頻譜中的向量B，則我們可視向量A的週期為N，且希望向量A可寫成下列式子：

1　N-1　　　 -ikm*2π/N
A(k) = ─ Σ　B(m) e　　　　　 ，k = 0～N-1
N　m=0

可算出：

N-1　　　 imk*2π/N
B(m) = Σ　A(k) e　　　　　　　，m = 0～N-1
k=0

這便是離散式傅立葉轉換的公式，而前一個式子便是其反轉公式。但是式子中的e^(i2π/N)看起來和寫起來均很不方便，因此我們定義

i (2π/N)
ω = e

則前述轉換公式可寫成：

N-1　　　　mk
B(m) = Σ　A(k) ω　　　　　，m = 0～N-1
k=0

1　N-1　　　　-km
A(k) = ─ Σ　B(m) ω　　　 ，k = 0～N-1
N　m=0

這個ω便是乘法單位元素的主要N次方根（Principle Nth Root），它具有下列特性：

2004/11/9 下午 09:55:44

　　　0　　　 1　　　 2　　　　　　N-1　　　 N
ω ＝1, ω ≠1, ω ≠1, …, ω　 ≠1, ω ＝1

1　　2　　　　N-1
1 + ω + ω + … ω　　= 0

N/2
ω　 ＝ -1　　，當 N = 2k 時

這些特性對於後面要說明的快速傅立葉轉換演算法（簡稱FFT）是很重要的，請各位要記下來。其實這樣的式子可以更進一步擴展至通用的Commutative Ring(R,+,*,0,1)數學領域上，但為方便說明，以下我們還是只局限在複數領域上。所謂Commutative Ring，其定義為：

1.具有分配律 (a+b)+c = a+(b+c)、 (a*b)*c = a*(b*c)
2.具有結合律 (a+b)*c = a*c + b*c、a*(b+c) = a*b + a*c
3.0為加法單位元素 a + 0 = 0 + a = a
4.1為乘法單位元素 a * 1 = 1 * a = a
5.加法可交換 a + b = b + a
6.具有加法反元素 a + b = 0
7.乘法可交換 a * b = b * a

若去除第7項特性，則稱為Ring；若加上具有乘法反元素特性，則稱為Field。通常傅立葉轉換後的函數，可表示成下式：

F(u) = ｜F(u)｜ e^iφ(u)

｜F(u)｜為傅立葉頻譜大小（Fourier Spectrum），φ(u)則為相位角度（Phase Angle），而｜F(u)｜^2則通常稱為次方頻譜（Power Spectrum）或頻譜密度（Spectral Density）。其定義如下：

2　　　　　　 2
｜F(u)｜= sqrt(F(u).Real　+ F(u).Image )

-1
φ(u) = tan　 (F(u).Image/F(u).Real)

上式之Real、Image分別表示該數值的實數及虛數部份。這些定義，在傅立葉頻譜的顯示上是相當常用的。而對於二維影像而言，相位角度包含了大部份影像的邊緣（Edge）資訊，這點對於影像處理上是很有用的。

至於傅立葉轉換有那些特性呢？下面我們便列出幾個一維傅立葉轉換的特性。

1.線性（Linearity）－對於任兩個時序領域向量A(x)、B(x)，則

FT(c1*A(x)+c2*B(x)) = c1*FT(A(x)) + c2*FT(B(x))

其中FT表示傅立葉轉換，c1、c2為任何兩個複數值。

2.平移（Shifting）與調變（Modulation）－對於任一個時序領域向量A(x)，傅立葉轉換後的向量為B(x)，則

FT(A(x-k)) = B(x) * ω^xk
FT(A(x)*ω^-xk) = B(x-k)

3.共輒對稱性（Conjudate Symmetry）－對於任一個時序領域實數向量A(x)，若其傅立葉轉換後的向量為B(y)，則B(y)與B(N-y)互為共軛數，其中N為向量長度。由於這點在傅立葉轉換計算中亦很有用，因此我們特別再加以說明。假設有兩個實數向量f(n)、g(n)，經傅立葉轉換後的向量為F(n)、G(n)，若我們將f、g組合成一複數向量x(n)=f(n)+ig(n)，經傅立葉轉換後的向量為X(n)，則依線性化特性，可得

X(n) = F(n) + i G(n)

也就是說

X(n).Real　= F(n).Real　- G(n).Image
X(n).Image = F(n).Image + G(n).Real

因f(n)、g(n)均為實數，故F(n)、F(N-n)及G(n)、G(N-n)均互為補數，可得：

X(N-n).Real　=　F(n).Real　+ G(n).Image
X(N-n).Image = -F(n).Image + G(n).Real

解這兩組式子便可得到

F(n).Real　= ( X(n).Real　+ X(N-n).Real ) / 2
F(n).Image = ( X(n).Image – X(N-n).Image) / 2
G(n).Real　= ( X(n).Image + X(N-n).Image) / 2
G(n).Image = (-X(n).Real　+ X(N-n).Real ) / 2

由這個結果可知，一次傅立葉轉換計算，其實是可以算出兩個實數向量的傅立葉轉換值的（記得F、G因有共輒數關係，上式真正的計算只需一半）。

2004/11/9 下午 09:57:13

4.乘積（Product）與迴旋（Convolution）－對於任兩個時序領域向量a(x)、b(x)，則

FT(a(x)‧b(x)) = FT(a(x)) × FT(b(x)) / N
FT(a(x)×b(x)) = FT(a(x)) ‧ FT(b(x))

其中‧表向量乘積、×表向量迴旋。

5.滯後乘積（Lagged Products）－滯後乘積為關連運算（Correlation）的主要部份。對於任兩個時序領域向量f(n)、g(n)，經傅立葉轉換後的向量為F(n)、G(n)，則

N-1
Σ　f(n)g(n-k) = F(n)‧G*(n)
n=0

其中G*(n)表G(n)的共輒數。

6.時序差（Cyclic Difference）－對於任一個時序領域向量A(x)，若B(x)=A(x)-A((x-1) mod N)，則FT(B(x)) = FT(A(x)) * (1 – ω^x)。

當然，傅立葉轉換還有其他特性，如旋轉、放大、積分、微分、Parseval理論等等，這裡只是列出其最基本的特性而已。

至於傅立葉轉換到底有什麼用處？傅立葉轉換除了可用做一般數位訊號或影像在傅立葉頻譜上的處理，如影像銳化、平滑、特徵擷取、資料壓縮、濾波等等之外，也可以用來解決一些難解的數學問題。例如線性微分方程式（Linear Differential Equation），經傅立葉轉換會成為一般的代數方程式（Algebraic Equation），要解決這樣的方程式便容易多了。此外，也可加速如多項式迴旋、矩陣乘法、代數運算、高精度乘法等數學計算上的速度。其他的應用例子還很多，這些便由讀者自行看書並親自試驗去體會了。在本文的最後，我們也會列出一些實際應用的例子，供讀者參考。

第二節　快速傅立葉轉換演算法

在提到快速傅立葉轉換演算法之前，我們要先算出傳統傅立葉轉換所需的時間（Time Complexity），以便和快速傅立葉轉換做一番比較，並進行驗證的工作。傳統計算法所需的時間其實很簡單求的，因要向量長度為n時需要計算n個元素值，而每個元素值均需cn個乘法和cn個加法，其中c為一常數，因此它的時間便是O(n^2)。傳統的傅立葉轉換計算，其程式如下：

#define PI 3.141592653589793
typedef struct { float real, image; } complex;

void DFT(int data_no,complex *in_data,complex *out_data)
{
/* ————————————————————
作用：直接式的離散傅立葉轉換
參數：in_data = 輸入資料
data_no = 資料數目
輸出：out_data = 輸出資料
———————————————————— */
int m, k;
float w_cos, w_sin, angle_step, angle;

angle_step = 2*PI/data_no;　　/* 乘法單位元素根ω之間隔角度 */
for (m=0; m<data_no; m++)　　 /* 計算B(m) */
{
out_data[m].real = in_data[0].real;
out_data[m].image = in_data[0].image;
for (k=1; k<data_no; k++)　/* 計算A(k)*ω^(mk) */
{
angle = (float) m*k*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
out_data[m].real += in_data[k].real * w_cos – in_data[k].image * w_sin;
out_data[m].image += in_data[k].real * w_sin + in_data[k].image * w_cos;
}
}
}
至於反向傅立葉轉換程式要如何寫呢？當然我們也可以直接從定義來寫，但其實正向與反向傅立葉轉換程式是相當類似的，因此我們提出了兩種做法：

1.將乘法單位元素根ω之間隔角度2π/N由正變負，然後進行同樣的傅立葉轉換，最後將回傳值全部除以N即可（若一開始將輸入值全部除以N，其效果相同）。下面便是這個方法的反轉程式：

void IDFT(int data_no,complex *in_data,complex *out_data)
{
/* ————————————————————
作用：直接式的反向離散傅立葉轉換
參數：in_data = 輸入資料
data_no = 資料數目
輸出：out_data = 輸出資料
———————————————————— */
int m, k;
float w_cos, w_sin, angle_step, angle;

angle_step = -2*PI/data_no;　/* 乘法單位元素根ω之間隔角度改負 */
for (m=0; m<data_no; m++)　　/* 計算B(m) */
{
out_data[m].real = in_data[0].real;
out_data[m].image = in_data[0].image;
for (k=1; k<data_no; k++)　/* 計算A(k)*ω^(-mk) */
{
angle = (float) m*k*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
out_data[m].real += in_data[k].real * w_cos – in_data[k].image * w_sin;
out_data[m].image += in_data[k].real * w_sin + in_data[k].image * w_cos;
}
}
/* 傳回值全部除以N */
for (m=0; m<data_no; m++)
{
out_data[m].real /= data_no;
out_data[m].image /= data_no;
}
}

2.將輸入資料改成共輒數，呼叫傅立葉轉換程式後，再將輸出資料全部除以N，並轉成共輒數。以下為這個方法的反轉程式：

void IDFT(int data_no,complex *in_data,complex *out_data)
{
/* ————————————————————
作用：直接式的反向離散傅立葉轉換
參數：in_data = 輸入資料
data_no = 資料數目
輸出：out_data = 輸出資料
———————————————————— */
int m;

/* 將輸入資料全部改成共輒數 */
for (m=0; m<data_no; m++) in_data[m].image = -in_data[m].image;
/* 進行傅立葉轉換 */
DFT(data_no,in_data,out_data);
/* 將輸入資料改回原樣 */
for (m=0; m<data_no; m++) in_data[m].image = -in_data[m].image;
/* 傳回值全部除以N, 並改成共輒數 */
for (m=0; m<data_no; m++)
{
out_data[m].real /= data_no;
out_data[m].image = -out_data[m].image/data_no;
}
}
上面這兩種方式以第一種速度較快些，若能使用預先計算ω^k值方式，將此值均除以N，那麼傳回值便不必再全部除以N，速度會更快些。但重複程式碼也比較多些，這是無法兩全的事。而第二者卻有不重複程式碼的好處。為了能適應各種不同的傅立葉轉換演算法，且不必在反向傅立葉轉換程式中再寫同樣的傅立葉轉換程式碼，因此後面各演算法的反轉程式，我們將都採用第二種方法來撰寫。如果確實很在意執行速度的話，讀者不妨自行依傅立葉轉換程式碼，將之修改成第一種的反轉方式。

相對於一般直接運算的O(n^2)時間，快速傅立葉轉換所需的時間則只要O(nlogn)即可。至於快速傅立葉轉換的計算方式是怎麼來的呢？最早的時候是有人將傅立葉轉換公式視為矩陣運算：

N-1　　　　mk
B(m) = Σ　A(k) ω　　　　，m = 0～N-1
k=0

B　　= W * A

其中W為一二維矩陣，且W[m,k]=ω^mk。當N=2^L時，由ω之特性可將矩陣W分解成L個Sparse矩陣。所謂Sparse矩陣，指的一個矩陣T，其元素T[i,j]除了i=j時有值外，其餘均為0。如此每個二維Sparse矩陣對一維矩陣相乘，只需O(n)次運算，共有logn層，因此全部的運算時間只需O(nlogn)。然而分解W成為Sparse矩陣並不太容易寫成程式，因此後來便有人利用分割解決法來解決此一問題（請參閱拙作”法則分析”）。首先將傅立葉轉換公式寫成下列的多項式格式（注意，在這邊我們均假設N=2^k，如果項次不足時只要補係數為0之項次即可）：

N-1　　　 k　　　　　　 0　　N-1
B(x) = Σ　A(k) x　　　，x = ω ～ω
k=0

我們將奇數項放一堆，偶數項放一堆，成為：

B(x) = C(y)*x + D(y)　 ，y = x^2

其中

C(y) = A(N-1)*y^(N/2-1) + A(N-3)*y^(N/2-2) + … + A(1)
D(y) = A(N-2)*y^(N/2-1) + A(N-4)*y^(N/2-2) + … + A(0)

由於ω^2亦為乘法單位元素的主要N/2次方根，因此C(y)和D(y)均可視為A的奇次項和偶數項之傅立葉轉換後的值。又由於ω是乘法單位元素的主要N次方根，依其特性可知ω^(N/2)=-1，因此前式可綜整為：

j　　　　　　 2j　　j　　　 2j
B(ω　　 ) =　C(ω　) ω　+ D(ω　)

j+N/2　　　　 2j　　j　　　 2j
B(ω　　 ) = -C(ω　) ω　+ D(ω　)　　，j = 0～N/2-1

假設T(n)為計算長度n向量的傅立葉轉換時間，則依前面公式可得時間計算公式為：

T(n) = 2*T(n/2) + cn　　 ，c為一常數

由此公式可解得T(n)=O(nlogn)。因此利用分割解決法，便可以得到一個更有效率的演算法，我們將分割解決法的遞迴呼叫計算方式寫成下列程式：

void RecursiveFFT(int data_no,complex *data)
{
/* ————————————————————
作用：遞迴式的快速離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：資料數目必須為2的次方方可呼叫本函數
———————————————————— */
int i, j, new_pos;
complex temp;
float angle_step, angle, w_cos, w_sin;

if (data_no <= 2)
{
temp = data[1];
data[1].real = data[0].real – temp.real;
data[1].image = data[0].image – temp.image;
data[0].real += temp.real;
data[0].image += temp.image;
return;
}
angle_step = 2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */

/* 奇偶分堆, 偶在左, 奇在右 */
for (i=2; i<data_no; i+=2)
{
/* 將位置i資料調至i/2 */
new_pos = i/2;
temp = data[i];
for (j=i-1; j>=new_pos; j–) data[j+1] = data[j];
data[new_pos] = temp;
}
/* 呼叫自己計算子向量 */
RecursiveFFT(data_no/=2,data);
RecursiveFFT(data_no,data+data_no);
/* 合併出答案 */
for (i=0; i<data_no; i++)
{
angle = (float) i*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
new_pos = i+data_no;
temp.real = data[new_pos].real * w_cos – data[new_pos].image * w_sin;
temp.image = data[new_pos].real * w_sin + data[new_pos].image * w_cos;
data[new_pos].real = data[i].real – temp.real;
data[new_pos].image = data[i].image – temp.image;
data[i].real += temp.real;
data[i].image += temp.image;
}
}

void RecursiveIFFT(int data_no,complex *data)
{
/* ————————————————————
作用：遞迴式的快速反向離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：資料數目必須為2的次方才可呼叫本函數
———————————————————— */
int i, j, new_pos;
complex temp;
float angle_step, angle, w_cos, w_sin;

if (data_no <= 2)
{
temp = data[1];
data[1].real = (data[0].real – temp.real) / 2;
data[1].image = (data[0].image – temp.image) / 2;
data[0].real = (data[0].real + temp.real) / 2;
data[0].image = (data[0].image + temp.image) / 2;
return;
}
angle_step = -2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 奇偶分堆, 偶在左, 奇在右 */
for (i=2; i<data_no; i+=2)
{
/* 將位置i資料調至i/2 */
new_pos = i/2;
temp = data[i];
for (j=i-1; j>=new_pos; j–) data[j+1] = data[j];
data[new_pos] = temp;
}
/* 呼叫自己計算子向量 */
RecursiveIFFT(data_no/=2,data);
RecursiveIFFT(data_no,data+data_no);

/* 合併出答案 */
for (i=0; i<data_no; i++)
{
angle = (float) i*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
new_pos = i+data_no;
temp.real = data[new_pos].real * w_cos – data[new_pos].image * w_sin;
temp.image = data[new_pos].real * w_sin + data[new_pos].image * w_cos;
data[new_pos].real = (data[i].real – temp.real) / 2;
data[new_pos].image = (data[i].image – temp.image) / 2;
data[i].real = (data[i].real + temp.real) / 2;
data[i].image = (data[i].image + temp.image) / 2;
}
}

雖然我們已經寫出了一個快速的傅立葉轉換計算程式，而且經實際試驗結果，確實比前面直接計算的方式快了相當多，但由於它是遞迴呼叫方式，在效率上仍未能達到最佳。而且遞迴呼叫之前的分堆資料搬移動作也佔了相當多的時間，這點在資料數增多時便開始逐漸浮現。因此我們必須找出其非遞迴呼叫的程式，並設法去除分堆的資料搬移動作，直接併入資料計算中，如此便可再提高快速傅立葉轉換的計算速度。但這點必須展開遞迴呼叫的程序，觀察其規律性，才能得到結論。不過在展開遞迴呼叫分析前，我們先說明另一種利用多項式除法的觀念來解決成非遞迴呼叫的方式，以便後面進行遞迴呼叫展開分析時，會比較容易懂些。以下我們便說明一下這個多項式除法的觀念。

依照前面所提及之多項式：

N-1　　　 k　　　　　　 0　　N-1
B(x) = Σ　A(k) x　　　，x = ω ～ω
k=0

要求B(ω^m) 的值，即是求B(x)除以(x-ω^m)之餘數值，但m=0～N-1，全部要除上N次也是很花時間的。由於其中有相當多的重複計算，因此我們可以取適當的分式積：

m1　　　 m2　　　　　 mk
(x-ω　) (x-ω　) … (x-ω　)

求出該式除B(x)的餘式R(x)後，再以各分式除R(x)來得到所要的餘數，這樣的計算速度將會快些。當N=2^k次方時（若不足可補係數0，前面已提過），由於所有(x-ω^m)的乘積將為(x^N-1)，因此我們可以利用二分的方式持續進行分式除法（有點類似分割解決法）。例如N=8時，各層所使用的除式可列舉如下：

x^8-1
┌───────┴──────┐
x^4-1　　　　　　　　　　　　x^4-ω^4
┌───┴───┐　　　　　　┌───┴───┐
x^2-1　　　　 x^2-ω^4　　　 x^2-ω^2　　　 x^2-ω^6
┌─┴─┐　　　┌─┴─┐　　┌─┴─┐　　　┌─┴─┐
x-ω^0　x-ω^4　x-ω^2　x-ω^6　x-ω　x-ω^5　x-ω^3　x-ω^7

由其結果可發現，利用此種二分分式除法方式，不論如何調換分式，最後得到的結果都不會是我們所要的順序，但我們卻可排出一種位元反轉的順序出來。所謂位元反轉，指的是數字中的位元順序完全調換，例如3=011b，其位元反轉後的值便是110b=6。上例得到的結果便是：

B[0] B[4] B[2] B[6] B[1] B[5] B[3] B[7]

B[m]即B(ω^m)。讀者可以檢驗看看這個順序便是位元反轉順序。利用這個位元反轉順序的關係，我們便可以很容易將這些資料重排回我們所要的順序。但若要有效率的重排，不妨使用Cooley、Lewis和Welch所提出的演算法：

j = 0;
for (i=1; i<N-1; i++)　/* N為資料數 */
{
k = N/2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
if (i < j) 資料i和資料j對調
}
接下來我們要說明一下上述這些分項式的除法要如何計算。對於一個2r-1次多項式除以(x^r-c)的餘式，其實便是：

r-1　　　　　　　　　　　　　　r-1
A(r-1)x　 + … + A(0) + c * (A(2r-1)x　 + … + A(r))

以上述例子來說，假設各項次係數為[A7 A6 A5 A4 A3 A2 A1 A0]，則除以(x^4-ω^4)的餘數便是：

[A3 A2 A1 A0] + ω^4 [A7 A6 A5 A4]
= [A3 A2 A1 A0] – [A7 A6 A5 A4]
= [A3-A7 A2-A6 A1-A5 A0-A4]

至於各層除式所應乘上的ω^k，k值應為何呢？上面已提過，各層的分式均保持著位元反轉的順序，因此當分式在該層第i個位置時（由0計算起），k值便是i的位元反轉值。不過在實際計算時，是兩個分式為一組計算的，除了儲存空間上的考量外，也因可使用同一個ω^k值來加快速度。

依照上述的這些特性，我們已經得到一種非遞迴呼叫方式的快速傅立葉轉換演算法。至於其時間Order如何？由於上述之二分結果共有logn層，每層只需cn次運算，因此其時間也是O(nlogn)的，有興趣的讀者不妨做一下較正式的證明。下面便是我們依照上述的多項式除法觀念所寫成的程式：

void PolynomialFFT(int data_no,complex *data)
{
/* ————————————————————
作用：多項式除法的快速離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int i, j, k, level, depth, depart_no, sub_data_no, windex;
int pre_index, post_index;
float w_cos, w_sin, angle_step, angle;
complex temp;

if (data_no < 2) return;
angle_step = 2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 計算層次數 */
depth = 15;
sub_data_no = data_no;
while ((sub_data_no & (int)0x8000) == 0)
{
sub_data_no <<= 1;
depth–;
}
data_no = 1 << depth;
/* 開始轉換 */
depart_no = sub_data_no = data_no;
for (level=0; level<depth; level++)
{
depart_no /= 2;
windex = 0;
/* 處理每個多項式除法 */
for (j=0; j<data_no; j+=sub_data_no)
{
/* 計算該多項式所應乘上的ω^k值 */
i = 0;
for (k=0; k<depth; k++)
i = (i << 1) | ((windex >> k) & 1);
angle = i*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
windex += 2;

　　 /* 處理該多項式的除法 */
for (i=0; i<depart_no; i++)
{
pre_index = i + j;
post_index = pre_index + depart_no;
temp.real = data[post_index].real * w_cos – data[post_index].image * w_sin;
temp.image = data[post_index].real * w_sin + data[post_index].image * w_cos;
data[post_index].real = data[pre_index].real – temp.real;
data[post_index].image = data[pre_index].image – temp.image;
data[pre_index].real += temp.real;
data[pre_index].image += temp.image;
}
}
sub_data_no = depart_no;
}
/* 利用位元反轉重排資料 */
j = 0;
for (i=1; i<data_no-1; i++)
{
k = data_no/2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
if (i < j)
{
temp = data[j];
data[j] = data[i];
data[i] = temp;
}
}
}

void PolynomialIFFT(int data_no,complex *data)
{
/* ————————————————————
作用：多項式除法的快速反向離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int m;

/* 將輸入資料全部改成共輒數 */
for (m=0; m<data_no; m++) data[m].image = -data[m].image;
/* 進行傅立葉轉換 */
PolynomialFFT(data_no,data);
/* 傳回值全部除以N, 並改成共輒數 */
for (m=0; m<data_no; m++)
{
data[m].real /= data_no;
data[m].image = -data[m].image/data_no;
}
}

讀者可在本節最後面所列的各演算法程式測試比較表發現，這個程式又比遞迴呼叫的方式快了相當多，只不過每次要算ω^k值時都必須做一次位元反轉，比較稍費時些。

現在我們再回到原來的遞迴呼叫演算法上。原來的遞迴呼叫演算法的遞迴呼叫部份，主要來自於下式（前面已列過）：

j　　　　　　 2j　　j　　　 2j
B(ω　　 ) =　C(ω　) ω　+ D(ω　)

j+N/2　　　　 2j　　j　　　 2j
B(ω　　 ) = -C(ω　) ω　+ D(ω　)　　，j = 0～N/2-1

由於遞迴呼叫計算時，計算資料來源和計算後資料存回的位置都不一樣，因此有人便想辦法調換輸入資料的順序，使得計算資料來源和計算後資料存回的位置能夠一致，如此只要用迴圈的方式一層一層計算，不必再經由遞迴呼叫的方式計算，也不用進行資料分堆搬移的動作。經過展開分析的結果，發現輸入資料的順序應該改成上面我們提到的位元反轉順序，而計算公式也改成下式（在第k層時，k由0算起）：
B[i　　] = B[i] + B[i+2^k] ω^(i*N/2^(k+1))
B[i+2^k] = B[i] – B[i+2^k] ω^(i*N/2^(k+1))　，i=0～2^k-1

由於上面公式這種交叉運算存回的特性，因此一般我們便稱之為Butterfly法，而又由於上面的公式是經由時序領域的資料分析所得來的，因此又稱之為DIT（Dicimation in Time）Butterfly法（亦稱Cooley-Tukey法）。以下便是這個方法的傅立葉轉換程式：

void DITFFT(int data_no,complex *data)
{
/* ————————————————————
作用：DIT式的快速離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int i, j, k, level, depth, depart_no, sub_data_no, wstep;
int pre_index, post_index;
float w_cos, w_sin, angle_step, angle;
complex temp;

if (data_no < 2) return;
angle_step = 2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 計算層次數 */
depth = 15;
sub_data_no = data_no;
while ((sub_data_no & (int)0x8000) == 0)
{
sub_data_no <<= 1;
depth–;
}
data_no = 1 << depth;
/* 利用位元反轉重排輸入資料 */
j = 0;
for (i=1; i<data_no-1; i++)
{
k = data_no/2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
if (i < j)
{
temp = data[j];
data[j] = data[i];
data[i] = temp;
}
}
/* 開始轉換 */
depart_no = 1;
wstep = data_no;
for (level=0; level<depth; level++)
{
sub_data_no = depart_no * 2;
wstep /= 2;
for (i=0; i<depart_no; i++)
{
/* 計算該butterfly所應乘上的ω^k值 */
angle = i*wstep*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
for (j=0; j<data_no; j+=sub_data_no)
{
pre_index = i + j;
post_index = pre_index + depart_no;
temp.real = data[post_index].real * w_cos –
data[post_index].image * w_sin;
temp.image = data[post_index].real * w_sin +
data[post_index].image * w_cos;
data[post_index].real = data[pre_index].real – temp.real;
data[post_index].image = data[pre_index].image – temp.image;
data[pre_index].real += temp.real;
data[pre_index].image += temp.image;
}
}
depart_no *= 2;
}
}
void DITIFFT(int data_no,complex *data)
{
/* ————————————————————
作用：DIT式的快速反向離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int m;

/* 將輸入資料全部改成共輒數 */
for (m=0; m<data_no; m++) data[m].image = -data[m].image;
/* 進行傅立葉轉換 */
DITFFT(data_no,data);
/* 傳回值全部除以N, 並改成共輒數 */
for (m=0; m<data_no; m++)
{
data[m].real /= data_no;
data[m].image = -data[m].image/data_no;
}
}

這個程式只比多項式除法的方式稍快一點點而已，主要是因為ω^k在計算時，不必每次都要位元反轉。

除了對輸入的時序領域資料進行分析以外，也有人直接對輸出的頻譜領域資料進行分析，而得到類似於上述的演算法，這種方法稱之為DIF（Dicimation in Frequency）Butterfly法（或稱Sande-Tukey法）。其分析過程是這樣子的。對於離散式傅立葉轉換的公式：

N-1　　　　mk
B(m) = Σ　A(k) ω　　　　　，m = 0～N-1
k=0

將之分解如下（同樣假設N=2^k的情況）：

N/2-1　　　 mk　　N-1　　　　mk
B(m) =　Σ　A(k) ω　 +　Σ　A(k) ω
k=0　　　　　　 k=N/2

N/2-1　　　　　　　　　　 mN/2　　 mk
=　Σ　〔 A(k) + A(k+N/2) ω　　 〕ω
k=0

於是

N/2-1　　　　　　　　　　 2mN/2　　　 2mk
B(2m)　 =　Σ　〔 A(k) + A(k+N/2) ω　　　〕　ω
k=0

N/2-1　　　　　　　　　　 (2m+1)N/2　 (2m+1)k
B(2m+1) =　Σ　〔 A(k) + A(k+N/2) ω　　　〕　ω
k=0

由於

2mN/2　　　　　(2m+1)N/2
ω　　　= 1，　ω　　　　　= -1
因此

　　　　　　N/2-1　　　　　　　　　　　 2mk
B(2m)　 =　Σ　〔 A(k) + A(k+N/2) 〕ω
k=0

N/2-1　　　　　　　　　　　 k　　 2mk
B(2m+1) =　Σ　〔 (A(k) – A(k+N/2)) ω　〕ω
k=0

這個公式便是DIF Butterfly法的由來。將其計算過程逐一列出，便可發現其運算過程恰與DIT Butterfly相反（在第k層時，k由0算起）：

B[i　　　　　] = B[i] + B[i+N/2^(k+1)]
B[i+N/2^(k+1)] = (B[i] – B[i+N/2^(k+1)]) ω^(i*2^k)，i=0～N/2^(k+1)

而所得結果便是位元反轉順序。以下便是DIF Butterfly演算法的程式：

void DIFFFT(int data_no,complex *data)
{
/* ————————————————————
作用：DIF式的快速離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int i, j, k, level, depth, depart_no, sub_data_no, wstep;
int pre_index, post_index;
float w_cos, w_sin, angle_step, angle;
complex temp, temp2;

if (data_no < 2) return;
angle_step = 2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 計算層次數 */
depth = 15;
sub_data_no = data_no;
while ((sub_data_no & (int)0x8000) == 0)
{
sub_data_no <<= 1;
depth–;
}
data_no = 1 << depth;
/* 開始轉換 */
depart_no = sub_data_no = data_no;
wstep = 1;
for (level=0; level<depth; level++)
{
depart_no /= 2;
for (i=0; i<depart_no; i++)
{
/* 計算該butterfly所應乘上的ω^k值 */
angle = i*wstep*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
for (j=0; j<data_no; j+=sub_data_no)
{
pre_index = i + j;
post_index = pre_index + depart_no;
temp.real = data[pre_index].real + data[post_index].real;
temp.image = data[pre_index].image + data[post_index].image;
temp2.real = data[pre_index].real – data[post_index].real;
temp2.image = data[pre_index].image – data[post_index].image;
data[post_index].real = temp2.real*w_cos – temp2.image*w_sin;
data[post_index].image = temp2.real*w_sin + temp2.image*w_cos;
data[pre_index] = temp;
}
}

sub_data_no = depart_no;
wstep *= 2;
}
/* 利用位元反轉重排輸入資料 */
j = 0;
for (i=1; i<data_no-1; i++)
{
k = data_no/2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
if (i < j)
{
temp = data[j];
data[j] = data[i];
data[i] = temp;
}
}
}

void DIFIFFT(int data_no,complex *data)
{
/* ————————————————————
作用：DIF式的快速反向離散傅立葉轉換
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int m;

/* 將輸入資料全部改成共輒數 */
for (m=0; m<data_no; m++) data[m].image = -data[m].image;
/* 進行傅立葉轉換 */
DIFFFT(data_no,data);
/* 傳回值全部除以N, 並改成共輒數 */
for (m=0; m<data_no; m++)
{
data[m].real /= data_no;
data[m].image = -data[m].image/data_no;
}
}

由於和DIT Butterfly法計算方式相當類似，因此這兩種方法的執行時間也相差不多（雖然DIF Butterfly法略快一點點）。那麼為什麼要再找出DIFButterfly法呢？其實這兩種方法在一般的情況下的傅立葉轉換是相同的，但在某些條件下的傅立葉轉換卻是有不同的優點的。例如傅立葉轉換修剪（Pruning）問題，當輸入訊號或輸出訊號只有一部份有值，其餘均為0的情況下，某些運算是可以省略的。而在輸入訊號較少時，可由DIT Butterfly法修改出更有效率的演算法，在輸出訊訊號較少時，可由DIF Butterfly法修改出更有效率的演算法，因此這兩種不同的方法，仍然是有其不同的應用價值的。但作者對這些特殊條件下的傅立葉轉換並不想多談（主要也是因為太多了），因此有興趣的讀者不妨多找這方面的書來看，以看完本文後所打下來的基礎，相信看起來應不吃力才對。

其實上述的DIT Butterfly法和DIF Butterfly法還有很多的形式，主要是在計算次序的安排上。例如將兩兩交叉運算存回的位置安排成每層都相同，較有利於硬體上的製作；或是利用額外的空間暫存，使得輸出入資料的順序都不必做位元反轉順序重排等等。不過這些改變對於程式執行效率上來講，改進都不大。另外必須提的一點是，資料的位元反轉順序重排，可直接加在輸入資料，亦可直接加在輸出資料，兩者的結果是相同的。因此不管DIT Butterfly法或DIF Butterfly法，我們均可將位元反轉順序重排的程式調至前頭或後頭。

雖然我們已經列出各種不同現有的演算法及程式，但實際寫作時卻還有許多變化來加速，例如改良指標處理方式、或將ω^k的cos、sin計算提到迴圈外面等等，但減少的時間均相當有限。對於後者，我們也以DIF Butterfly為例修改程式如下，供讀者參考。其中將記錄ω^k值的變數做成公用變數，以便將來能在二維的傅立葉轉換中使用。其實DIT Butterfly或多項式除法等演算法，將ω^k的cos、sin計算提到迴圈外面之後，這幾種演算法的計算速度便相差無幾了，只是要多使用些記憶體罷了。當然，我們也考慮到了記憶體不夠的情況。只不過增加N/2的記憶體只減少N/2的cos、sin計算，是否划算，便取決於讀者了。但和一維不同的是，通常在二維的傅立葉轉換計算上都會先算好並記錄ω^k的值，主要因使用O(n)的記憶體可節省O(n^2)的時間，是非常划算的事。

2004/11/9 下午 10:10:05

static Keep_No=0;
static complex *W_Root;
void FFT(int data_no,complex *data)
{
/* ————————————————————
作用：加強版的快速離散傅立葉轉換（採DIF Butterfly法）
參數：data_no = 資料數目，若為0表示要清除函數內使用之記憶體
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int i, j, k, level, depth, depart_no, sub_data_no, wstep;
int pre_index, post_index;
float w_cos, w_sin, angle_step, angle;
complex temp, temp2;

if (data_no == 0)
{
if (Keep_No != 0) free(W_Root);
Keep_No = 0;
}
if (data_no < 2) return;
angle_step = 2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 計算層次數 */
depth = 15;
sub_data_no = data_no;
while ((sub_data_no & (int)0x8000) == 0)
{
sub_data_no <<= 1;
depth–;
}
data_no = 1 << depth;
/* 取得並計算ω^k值 */
if (data_no != Keep_No)
{
if (Keep_No != 0) free(W_Root);
Keep_No = 0;
/* 只需算一半即可 */
sub_data_no = data_no / 2;
W_Root = (complex*) malloc(sub_data_no*sizeof(complex));
if (W_Root != NULL)
{
Keep_No = data_no;
for (i=0; i<sub_data_no; i++)
{
angle = i*angle_step;
W_Root[i].real = cos(angle);
W_Root[i].image = sin(angle);
}
}
}
/* 開始轉換 */
depart_no = sub_data_no = data_no;
wstep = 1;
for (level=0; level<depth; level++)
{
depart_no /= 2;
for (i=0; i<depart_no; i++)
{
/* 計算該butterfly所應乘上的ω^k值 */
j = i*wstep;
if (Keep_No != 0)
{
w_cos = W_Root[j].real;
w_sin = W_Root[j].image;
}
else
{
angle = j*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
}

for (j=0; j<data_no; j+=sub_data_no)
{
pre_index = i + j;
post_index = pre_index + depart_no;
temp.real = data[pre_index].real + data[post_index].real;
temp.image = data[pre_index].image + data[post_index].image;
temp2.real = data[pre_index].real – data[post_index].real;
temp2.image = data[pre_index].image – data[post_index].image;
data[post_index].real = temp2.real*w_cos – temp2.image*w_sin;
data[post_index].image = temp2.real*w_sin + temp2.image*w_cos;
data[pre_index] = temp;
}
}
sub_data_no = depart_no;
wstep *= 2;
}
/* 利用位元反轉重排輸入資料 */
j = 0;
for (i=1; i<data_no-1; i++)
{
k = data_no/2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
if (i < j)
{
temp = data[j];
data[j] = data[i];
data[i] = temp;
}
}
}

void IFFT(int data_no,complex *data)
{
/* ————————————————————
作用：加強版的快速反向離散傅立葉轉換（採DIF Butterfly法）
參數：data_no = 資料數目
data = 輸入資料
輸出：data = 輸出資料
限制：1.資料數目必須為2的次方方可呼叫本函數
2.資料數目不可超過16384個
———————————————————— */
int m;

/* 將輸入資料全部改成共輒數 */
for (m=0; m<data_no; m++) data[m].image = -data[m].image;
/* 進行傅立葉轉換 */
FFT(data_no,data);
/* 傳回值全部除以N, 並改成共輒數 */
for (m=0; m<data_no; m++)
{
data[m].real /= data_no;
data[m].image = -data[m].image/data_no;
}
}

也許讀者有個問題：如果給定的向量長度不是2的次方，而且也無法預留補0係數的空間，那麼要如何運用上面的程式呢？在這個情況下，便無法使用快速傅立葉轉換的演算法了，只好使用最原始的計算方式了。不過由於前面所列的DFT程式，其ω^k值的cos、sin運算相當頻繁，嚴重影響了執行速度，因此我們便將此部份提到迴圈外，而程式便修改如下，但執行效率仍不夠好，這是沒有辦法的事（其實仍是有方法的，但由於頗為複雜，故不在此討論）。
void EnhancedDFT(int data_no,complex *in_data,complex *out_data)
{
/* ————————————————————
作用：加強版的直接式離散傅立葉轉換
參數：data_no = 資料數目，若為0表示要清除函數內使用之記憶體
in_data = 輸入資料
輸出：out_data = 輸出資料
———————————————————— */
static keep_no=0;
static complex *w_root;
int m, k, i;
float w_cos, w_sin, angle_step, angle;

if (data_no == 0)
{
if (keep_no != 0) free(w_root);
keep_no = 0;
}
if (data_no < 2) return;
angle_step = 2*PI/data_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 取得並計算ω^k值 */
if (data_no != keep_no)
{
if (keep_no != 0) free(w_root);
keep_no = 0;
w_root = (complex*) malloc(data_no*sizeof(complex));
if (w_root != NULL)
{
keep_no = data_no;
for (i=0; i<data_no; i++)
{
angle = i*angle_step;
w_root[i].real = cos(angle);
w_root[i].image = sin(angle);
}
}
}
for (m=0; m<data_no; m++)　　 /* 計算B(m) */
{
out_data[m].real = in_data[0].real;
out_data[m].image = in_data[0].image;
for (k=1; k<data_no; k++)　/* 計算A(k)*ω^(mk) */
{
if (keep_no != 0)
{
i = (int) ((long) m*k % data_no);
w_cos = w_root[i].real;
w_sin = w_root[i].image;
}
else
{
angle = (float) m*k*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
}
out_data[m].real += in_data[k].real * w_cos – in_data[k].image * w_sin;
out_data[m].image += in_data[k].real * w_sin + in_data[k].image * w_cos;
}
}
}

void EnhancedIDFT(int data_no,complex *in_data,complex *out_data)
{
/* ————————————————————
作用：加強版的直接式反向離散傅立葉轉換
參數：data_no = 資料數目
in_data = 輸入資料
輸出：out_data = 輸出資料
———————————————————— */
int m;

/* 將輸入資料全部改成共輒數 */
for (m=0; m<data_no; m++) in_data[m].image = -in_data[m].image;
/* 進行傅立葉轉換 */
EnhancedDFT(data_no,in_data,out_data);
/* 將輸入資料改回原樣 */
for (m=0; m<data_no; m++) in_data[m].image = -in_data[m].image;
/* 傳回值全部除以N, 並改成共輒數 */
for (m=0; m<data_no; m++)
{
out_data[m].real /= data_no;
out_data[m].image = -out_data[m].image/data_no;
}
}

現在我們實際對以上七種傅立葉轉換程式的執行速度做一番測試，以便讓讀者能直接從這些測試數據裡了解到各個程式的優劣。下表便是上述幾種傅立葉轉換計算程式的平均執行時間比較表（取100次平均），內部的數據為秒。

┌───────┬───┬───┬───┬───┐
│資料數目　　　│　 64 │　256 │ 1024 │ 4096 │
├───────┼───┼───┼───┼───┤
│DFT　　　　　 │　0.13│　2.07│ 33.62│541.43│
├───────┼───┼───┼───┼───┤
│RecursiveFFT　│　0.01│　0.09│　1.13│ 15.54│
├───────┼───┼───┼───┼───┤
│PolynomialFFT │　0.00│　0.02│　0.09│　0.43│
├───────┼───┼───┼───┼───┤
│DITFFT　　　　│　0.00│　0.02│　0.09│　0.41│
├───────┼───┼───┼───┼───┤
│DIFFFT　　　　│　0.00│　0.02│　0.09│　0.40│
├───────┼───┼───┼───┼───┤
│FFT　　　　　 │　0.00│　0.01│　0.07│　0.34│
├───────┼───┼───┼───┼───┤
│EnhancedDFT　 │　0.04│　0.60│　9.83│160.16│
└───────┴───┴───┴───┴───┘
表一　各傅立葉轉換計算程式執行時間比較表

最後要再提的一點便是，根據第一節所提的傅立葉轉換特性三，我們可利用FFT一次找出兩個實數向量的傅立葉轉換值。由於在時間領域中，我們通常是使用實數向量，因此這個方法可使我們在進行多個向量的傅立葉轉換時，可省下近一半的時間。至於如何寫成程式，只是呼叫上述傅立葉轉換程式後的一些簡單計算而已，我們留給讀者當一個習題好了。但是要注意的是，這個特性在反轉過程中便無法適用，而對於二維的傅立葉轉換，也只能節省1/4的時間而已。
第三節　二維傅立葉轉換

關於二維傅立葉轉換，其定義如下：

∞　　　　　i2π(ux+vy)
F(u,v) = ∫∫ f(x,y) e　　　　　 dxdy
-∞

∞　　　　　i2π(xu+yv)
f(x,y) = ∫∫ F(u,v) e　　　　　 dudv
-∞

亦可定義成：

∞　　　　　i(ux+vy)
F(u,v) = ∫∫ f(x,y) e　　　 dxdy
-∞

1　　∞　　　　 i(xu+yv)
f(x,y) = ───∫∫ F(u,v) e　　　 dudv
4π^2　-∞

至於二維離散傅立葉轉換之定義則為下式：

M-1 N-1　　　　 i2π(sm/M+tn/N)
B(s,t) = Σ　Σ　A(m,n) e　　　　　　　 ，s = 0～M-1、t = 0～N-1
m=0 n=0

1　M-1 N-1　　　　 -i2π(ms/M+nt/N)
A(m,n) = ─　Σ　Σ　B(s,t) e　　　　　 ，m = 0～M-1、n = 0～N-1
MN　s=0 t=0

由於上式可寫成：

M-1　 N-1　　　　 i2πtn/N　 i2πsm/M
B(s,t) = Σ 〔 Σ　A(m,n) e　　　　〕e
m=0　 n=0

1 M-1　 1 N-1　　　　 -i2πnt/N　 -i2πms/M
A(m,n) = ─ Σ〔 ─ Σ　B(s,t) e　　　　 〕e
M s=0　 N t=0

因此我們可視二維傅立葉轉換為先對行或列進行一維傅立葉轉換後，再對列或行進行一維傅立葉轉換的結果。利用第二節所述之快速傅立葉轉換，我們便可以得到一個O(n^2logn)的演算法。下面便是這個演算法的程式：

void Decomposed2DFFT(int x_no,int y_no,complex **data)
{
/* ————————————————————
作用：行列分解法快速二維離散傅立葉轉換
參數：x_no = 列數，若為0表示要清除函數內使用之記憶體
y_no = 行數
data = 輸入資料
輸出：data = 輸出資料
限制：1.行數及列數均必須為2的次方方可呼叫本函數
2.行數或列數不可超過16384個
———————————————————— */
int i, j, k, level, depth, depart_no, sub_data_no, wstep;
int pre_index, post_index;
float w_cos, w_sin, angle_step, angle;
complex temp, temp2;

if (x_no == 0)
{
if (Keep_No != 0) free(W_Root);
Keep_No = 0;
}
/* 先對各行進行一維傅立葉轉換 */
for (line=0; line<y_no; line++) FFT(x_no,data[line]);
/* 再對各列進行一維傅立葉轉換 */
if (y_no < 2) return;
angle_step = 2*PI/y_no;　 /* 乘法單位元素根ω之間隔角度 */
/* 計算層次數 */
depth = 15;
sub_data_no = y_no;
while ((sub_data_no & (int)0x8000) == 0)
{
sub_data_no <<= 1;
depth–;
}
y_no = 1 << depth;

/* 取得並計算ω^k值 */
if (y_no != Keep_No)
{
if (Keep_No != 0) free(W_Root);
Keep_No = 0;
sub_data_no = y_no / 2;
W_Root = (complex*) malloc(sub_data_no*sizeof(complex));
if (W_Root != NULL)
{
Keep_No = y_no;
for (i=0; i<sub_data_no; i++)
{
angle = i*angle_step;
W_Root[i].real = cos(angle);
W_Root[i].image = sin(angle);
}
}
}
/* 開始轉換各列 */
for (line=0; line<x_no; line++)
{
depart_no = sub_data_no = y_no;
wstep = 1;
for (level=0; level<depth; level++)
{
depart_no /= 2;
for (i=0; i<depart_no; i++)
{
/* 計算該butterfly所應乘上的ω^k值 */
j = i*wstep;
if (Keep_No != 0)
{
w_cos = W_Root[j].real;
w_sin = W_Root[j].image;
}
else
{
angle = j*angle_step;
w_cos = cos(angle);
w_sin = sin(angle);
}
for (j=0; j<y_no; j+=sub_data_no)
{
pre_index = i + j;
post_index = pre_index + depart_no;
temp.real = data[pre_index][line].real +
data[post_index][line].real;
temp.image = data[pre_index][line].image +
data[post_index][line].image;
temp2.real = data[pre_index][line].real –
data[post_index][line].real;
temp2.image = data[pre_index][line].image –
data[post_index][line].image;
data[post_index][line].real = temp2.real*w_cos –
temp2.image*w_sin;
data[post_index][line].image = temp2.real*w_sin +
temp2.image*w_cos;
data[pre_index][line] = temp;
}
}
sub_data_no = depart_no;
wstep *= 2;
}

　/* 利用位元反轉重排輸入資料 */
j = 0;
for (i=1; i<y_no-1; i++)
{
k = y_no/2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
if (i < j)
{
temp = data[j][line];
data[j][line] = data[i][line];
data[i][line] = temp;
}
}
}
}

void Decomposed2DIFFT(int x_no,int y_no,complex **data)
{
/* ————————————————————
作用：行列分解法快速二維離散傅立葉轉換
參數：x_no = 列數
y_no = 行數
data = 輸入資料
輸出：data = 輸出資料
限制：1.行數及列數均必須為2的次方方可呼叫本函數
2.行數或列數不可超過16384個
———————————————————— */
float divide;
complex *sub_data;
int x, y;

/* 將輸入資料全部改成共輒數 */
for (y=0; y<y_no; y++)
{
sub_data = data[y];
for (x=0; x<x_no; x++) sub_data[x].image = -sub_data[x].image;
}
/* 進行二維傅立葉轉換 */
Decomposed2DFFT(data_no,data);
/* 傳回值全部除以M*N, 並改成共輒數 */
divide = (float) x_no * y_no;
for (y=0; y<y_no; y++)
{
sub_data = data[y];
for (x=0; x<x_no; x++)
{
sub_data[x].real /= divide;
sub_data[x].image = -sub_data[x].image/divide;
}
}
}

上面這個程式是假設使用動態記憶體來存放資料，若使用的是靜態記憶體，那麼上面的函數宣告上便必須做一些更動，但程式實體並不需改變。為了方便讀者使用，作者也寫了一個通用型的二維動態記憶體取得及釋放函數如下：

void** Alloc2DArray(int xsize,int ysize,int item_size)
{
/* ————————————————————
作用：取得二維陣列動態記憶體
參數：xsize, ysize = 二維陣列大小
item_size = 陣列元素大小
傳回：記憶體空間指標，若為NULL表記憶體不夠
限制：xsize*item_size必須小於65536
———————————————————— */
int i, j;
void** ptr;

ptr = (void**) malloc(ysize*sizeof(void*));
if (ptr == NULL) return NULL;
for (i=0; i<ysize; i++)
{
ptr[i] = (void*) malloc(xsize*item_size);
if (ptr[i] == NULL)
{
/* 記憶體不夠, 釋放原取得之記憶體 */
for (j=0; j<i; j++) free(ptr[j]);
free(ptr);
return NULL;
}
}
return ptr;
}
void Free2DArray(int xsize,int ysize,void **ptr)
{
/* ————————————————————
作用：釋放二維陣列動態記憶體
參數：xsize, ysize = 二維陣列大小
ptr = 記憶體空間指標
備註：xsize實際上並沒用到，但為容易辨識使用，故加上此一參數
———————————————————— */
int i;

if (ptr == NULL) return;
for (i=0; i<ysize; i++) free(ptr[i]);
free(ptr);
}

例如要取得256×128的complex二維陣列記憶體空間，只需呼叫下列函數：

data = (complex**) Alloc2DArray(256,128,sizeof(complex));
if (data == NULL) 記憶體不夠…

要釋放時則呼叫下列函數：

Free2DArray(256,128,data);

這樣的函數對於二維資料的記憶體空間處理上是頗為方便的。

現在我們已經得到一個快速的二維傅立葉轉換演算法，但是否有更快的演算法？答案是有的，我們仍然可以利用第二節所述的分割解決法來得到更快的演算法，此種方法稱為Vector-Radix法。對於二維傅立葉轉換的公式：

M-1 N-1　　　　 i2π(sm/M+tn/N)
B(s,t) = Σ　Σ　A(m,n) e　　　　　　　 ，s = 0～M-1、t = 0～N-1
m=0 n=0

……. (未完) …….
……. 以下暫略, 可能不再寫下去 ………　青衫 ’97 3/29
2004/11/9 下午 10:15:42

第四節　傅立葉轉換的應用

一．傅立葉頻譜顯示

二．地程高度量測誤差模式

三．矩陣乘法

習題

1.試由傅立葉積分式，導出傅立葉轉換公式。
2.證明第一節所述之各種傅立葉轉換特性。
3.證明多項式除法、DIT Butterfly法、DIF Butterfly法之傅立葉轉換計算時間均為O(nlogn)。
4.證明2r-1次多項式除以(x^r-c)的餘式為

r-1　　　　　　　　　　　　　　r-1
A(r-1)x　 + … + A(0) + c * (A(2r-1)x　 + … + A(r))

其中A(n)為x^n的係數。
5.請以N=16為例，展開DIT Butterfly法和DIF Butterfly法的運算方式，並驗證第二節所列的運算公式無誤。
6.試寫出傅立葉轉換修剪的程式（包括輸入資料修剪和輸出資料修剪兩種）。
7.試利用傅立葉轉換的共輒對稱性，寫出一個同時計算兩個實數向量的傅立葉轉換程式。

參考文獻

1."Image Processing – Theroy, Algorithms, and Architectures"，Maher A. Sid-Ahmed，1995，ISBN 0-07-113697-5。
2."Digital Image Processing"，Rafael C. Gonzalez、Richard E. Woods，1992，開發。
3."Fractal Image Compression"，Progress Report，July 1992。
4."C Language Algorithms for Digital Signal Processing"，Paul M. Embree、Bruce Kimble，1991，開發。
5."Principles of Pictorial Information Systems Design"，Shi-Kuo Chang，1989，ISBN 0-13-710393-X。
6."The Fast Fourier Transform and Its Applications"，E. Oran Brigham，1988，ISBN 0-13-307547-8。
7."One-Dimensional Digital Signal Processing"，Chi-Tsong Chen，1985，歐亞。
8."The Art of Computer Programming"，Vol.2，Donald E. Knuth，1985，台北。
9."Discrete Fourier Transformation and Its Application to Power Spectra Estimation"，Nezih C. Geckinli、Davras Yavusz，1985，巨擘。
10."Two-Dimensional Fast Fourier Transforms in Image Processing"，Donald E. Lewis，Systems Research Laboratories, Inc.，1984。
11."Advanced Engineering Mathematics"，Erwin Kreyszig，1983，淡江。
12."Fast Transforms – Algorithms, Analyses, Applications"，Douglas F. Elliott、K. Ramamohan Rao，1982，ISBN 0-12-237080-5。
13."Fast Fourier Transform and Convolution Algorithms"，Henri J. Nussbaumer，1982，ISBN 0-387-11825-X。
14."2-D FFTs of Large Images with AP-120B"，Richard E. Twogood，Floating Point Systems Users Group Meeting，Feb. 1980。
15."Computer Image Processing and Recognition"，Ernest L. Hall，1979，儒林。
16."Digital Image Processing"，William K. Pratt，1978，儒林。
17."The Design and Analysis of Computer Algorithms"，Aho、Hopcroft、Ullman，1974，開發。
18."Terrain Correlation Suitability"，Wang Tang、Robert L. McClintock，SPIE Vol.2220。
19."數位影像處理"，余松煜、周源華、吳時光編著，1993，儒林。
20."數位信號處理"，王啟川、郭弘政編譯，1992，全華。
21."數位影像處理"，衛祖賞編著，1989，全華。
 

 

---

x264

發表： 2011 年 07 月 14 日 | 作者： vinsonhsieh | Filed under: 影片音樂格式 | 發表留言

http://www.hd.club.tw/thread-15891-1-1.html

http://zh.wikipedia.org/wiki/X264%E7%B7%A8%E7%A2%BC%E5%99%A8

我想首先應該先區分x264和H264有何不同

我們不要用太嚴謹的定義來講的話
H264/AVC,VC1,MPEG2這些可以說是一種標準(Standard)或是一種格式(format)
或是一種codec(encoding編碼和decoding解碼)
而x264是一種類似程式program以library的型式提供出來的一種編碼encoding的程式
而且它是免費的,你可以自由修改它的介面和程式內容
相對於要錢的nero,quicktime等
所以你播放h264影片需要一些解碼器例如coreavc,ffdshow等等而不是x264

所以利用x264可以把原本的H264的影片作處理
你可以保留原本的VIDEO ADUIO不動(這樣就沒有意義了)
或者次之VIDEO不動,把某些AUDIO取消或是重編碼讓檔案小一點
或者把VIDEO,RE-ENCODING,RESIZE(很多X264的1080P都不是1920X1080,
720P也不是1280X720)
所以出來的檔案和原本的檔案大小差別由很小到很大都可能
這些多少都有可能會產生和原本影片的差異
至於差異的多大,是否容易分辨,應該要依各影片而論
不果我個人是下載很多x264,主要是usenet就有一個x264群組
而且我現在的設備也不到1080p和次世代聲道
除非是有些大片不然我也懶的下三五十G的完整影片

---

陳鍾誠，於金門大學資訊工程系

發表： 2011 年 07 月 13 日 | 作者： vinsonhsieh | Filed under: 有趣跟有用的網站 | 發表留言

http://ccckmit.wikidot.com/

---

品管機遇原因跟非機遇原因

發表： 2011 年 07 月 13 日 | 作者： vinsonhsieh | Filed under: 管理相關轉貼 | 發表留言

http://hanchingchung.myweb.hinet.net/deming/2/21998065.htm

共同原因和特殊原因

(Common cause VS. Special cause )

劉振 教授

凡曾修習過品質管制圖的人，全都知道在管制界限以內的變異原因，稱為機遇原因(Chance Cause)，超出管制界限以外的變異原因，稱為非機遇原因(Assignable Cause)。這是Walter A. Shewhart博士在 1924年發明第一張品質管制圖時，就這樣命名的。當時的學者為崇敬Shewhart博士把統計學應用在工業上所作的貢獻，故亦有稱為Shew- hart圖者。

直到民國59年(1970年)，美國的品管大師 W. Edwards Deming博士首次應邀訪台，在他所發表的學術演講中，他提出了兩個嶄新的名詞，就是在管制界限以內變異的原因，稱為共同原因(Common Cause)，超出管制界限以外的變異原因，稱為特殊原因(Special Cause)。許多美英等國的品管書籍上，也逐漸採用這個名詞。尤其近幾年來，Deming博士時有新著問世，聲望陡增，已成為品管界中之牛耳。

我們當時都以為Deming博士把共同原因來代替機遇原因，特殊原因來代替非機遇原因，無非是取其淺顯達意，容易明白而已。而且一直以為這兩個名詞是相等的(即：共同原因 = 機遇原因，特殊原因 = 非機遇原因)，包括我自己在內，也認為是如此，沒有再去作進一步的研究。

近幾年來，美國新出版的品管書籍很多，也曾討論過這些名詞的問題，才知道這兩個名詞並不完全相等，因為Shewhart和Deming兩人的出發點，並不一樣。茲就讀書所得，說明於次。

Shewhart博士的觀念是著眼於製程中的變異，凡是在管制下的變異(Controlled Variation)，可以認為是穩定(Stable)而經常存在的，所以稱之為機遇原因，亦有人稱之為無法避免之原因(Unavoidable Cause)者。這種變異是製程中固有的(Inherent)的變異，如果要把它減少(或減小)，製程亦必隨之而變動。

至於製程中不能管制的變異(Uncontrolled Variation)，是隨時在變動，既不穩定，亦非經常存在，並不是製程中的一部分。這種不穩定的性質，使製程不能按照預定的目標操作，致產生過多的變異(Excess Variation)。這種變異原因稱為非機遇原因，應該可以把這種變異原因找出來，並加以消除掉，使製程回復正常。故這種原因亦有稱之為可避免之原因(Avoidable Cause)者。

上述兩種方法，都可以用來改進製程，但基本上是不同的。前者是去修改(Modify)一個經常而穩定的製程，而後者是去創造(Create)一個經常而穩定的製程。到底要採用那一種方法？就要看製程中的變異而定。因此，改進製程的首要步驟，是先要決定製程中有無不能管制的變異原因存在。而Shewhart的品質管制圖是根據機率原理和統計學而來的，正是判斷製程中有無不能管制的變異原因的有效利器。

至於Deming博士，他的著眼點是放在：誰該對這種變異負責？於是創造了共同原因和特殊原因這兩個名詞。

共同原因的變異是存在於製造系統(System)或管制良好的系統中。因為這種變異是屬於系統內的，所以應該由管制這一系統的人員來負責，那即是：經理人員─特別是高階層的經理人員。共同原因的變異僅祇能由經理階層人員採取行動，才能把它移除掉。

特殊原因的變異在本質上是局部的。Deming博士曾說過，他喜歡用〝特殊〞這個形容詞來指特定的一群作業員，或特定的某一位作業員，或特定的機器，或特定的局部環境等所造成的特殊原因。他並說，名詞並不重要，重要的是觀念。一般都能從發生問題的那一特定階層人員，就可自行找出原因，採取行動，並把困擾的問題消除掉。

現在把兩位大師的區分方法，用圖解來說明，作一比較如下。

Shewhart博士：

非機遇原因 (不能管制的變異)	機遇原因 (管制下的變異)
每一個都有明顯的影響	每一個的影響都很小

Deming博士

特殊原因	共同原因
本質上是局部的	系統的一部分

從上圖可知，特殊原因始終是不能管制的非機遇變異原因。這些非機遇原因可以由作業員或那一部分的領班來採取行動。

而共同原因卻不然，它可以是機遇原因，也可以是非機遇原因。無論在那一種情況，都是系統上發生的錯誤(Faults)，僅祇能由經理階層人員採取行動，才能改正過來。

Deming博士是著重在誰該對這種不同型式的變異負責，現在可以用一句話來說明他的觀念：減少任何一個品質特性(如厚度，或績效等)的變異，不管這個品質特性是否在管制狀態之下，甚至祇有少數幾個或沒有不良品產生，都認為是良好的管理。

〝零缺點〞並不能算夠好。工業界必須做到比符合規格還要來得好。經理人員要去研究製程，並且應該去找尋變異的來源，進而把它消除掉，以經常改進產品。品質管制圖正是用來找出這些來源的有效利器。當製程中的變異減小後，則零件將更為相似，產品亦將愈佳。這些都是Deming博士在1950年帶給日本人的金玉良言，日本人懂得這些話，經過了石川一郎(按：即石川馨氏之尊翁)及小柳賢一等繼續不斷地努力，才奠定了實施SQC成功的基礎。

現在，將再介紹另一個觀念(Concept)，那就是：工程上(Engineering)對變異的觀念，和Shewhart博士(包括Deming博士)對變異的觀念，兩者完全不同。工程上對變異的觀念，其目的是要求產品能符合規格，不管產品中的變異有多大，祇要在規格範圍以內，就儘可能地讓它去變化(Vary)。如果結果是在規格範圍以內，就認為是〝夠好〞了。而Shewhart博士觀念的目的，是要製程經常穩定一致(Consistency)，結果自然是產品儘可能地穩定一致了。因為它們的目的不同，結果也就隨之而不同。我們沒有必要把這兩種觀念加以協調。經理人員必須採用其中之一，作為生產製程的指導原則：僅祇要求符合規格，還是繼續不斷地改進製程。

經理人員在工業革命時期就開始採用第一種觀念了。大約經過了200年以後，這個目的仍沒有達到。因為我們僅僅把目標放在符合規格上，結果卻缺乏進步。這就沒有理由使人相信將來會有什麼進步的緣故。

另一方面，根據日本的經驗，應用了Shewhart博士的觀念之後，證明了持續不斷地改進製程，使得日本的工業產品提高了它們的品質，也提高了它們的生產力。結果，不再侷限於以符合規格為滿足，而是在持續不斷地改進途程中前進。

要產品完全符合規格，唯有持續不斷地改進製程，而且僅祇有當經理人員用言詞和行動來支持這個目標時，才會獲得品質和生產力的增加。

參考資料：

1. W. A. Shewhart,The Economic Control of Quality of Manufactured Product (Van Nostrand,1931,ASQC,1980)
2. E. L. Grant and R. S. Leavenworth,Statistical Quality Control (McGraw-Hill,5th,ed.,1980)
3. W. E. Deming,Quality,Productivity,and Competitive Position (MIT,1982)
4. D. J. Wheeler and D. S. Chambers,Understanding Statistical Process Control (Statistical Process Controls,Inc,1986)

 

---