本文首發于i春秋,本人為本文原作者,轉載請标明出處
序言
佛語中有箴言：坐亦禅,行亦禅,一花一世界,一葉一如來
世間萬物是複雜的，但是又是純粹的簡單的，從宏觀的花花世界到微觀的原子電子，萬物都在按照它的規律運作，而我們的先輩前人，一直都在用自己的方式與經驗，總結着萬物運作的規律。
不要誤會，本文并不是哲學論題的讨論，一個偶然的機會拜讀了知乎上的一篇關于頻域水印的問答。（
阿裡巴巴公司根據截圖查到洩露資訊的具體員工的技術是什麼？https://www.zhihu.com/question/50735753/answer/122593277）當中有對頻域數字水印的實作與讨論，在群裡有不少的朋友對此頗感興趣，于是我就想以後機會寫一個“從零開始的頻域變換到水印的完整解答”
如果閱讀到這裡你還不明白本文的主旨與目标是什麼，你可以直接跳轉到本文的最後一章節“魯棒盲水印”來檢視頻域水印是如何在版權保護中起作用并對抗各類版權偷盜者的攻擊的。
當然，本人并非全才，本文中多有疏漏錯誤還望各位讀者多多指正，文章将從最基礎的三角函數開始，一步一步地推到并演化到傅裡葉變換直到頻域簽名。
我相信學習如流水，希望讀者能在本文的循序漸進的推導過程中滿足對相關技術的求知欲，并對文中的不足與錯誤不吝賜教。
最後再次引用箴言的後半句作為導言的結束語：
春來花自青,秋至葉飄零,無窮般若心自在,語默動靜體自然
從三角形開始
有人說，上帝使用三角形創造了這個世界，這點我完全同意，三角形擁有如此之多的特性，足夠讓每一個探究者為之所着迷，它仿佛維系着幾何與世界的基石，誕生出數學中衆多的定律并在今天成就了我們的學術與技術大廈。
當然，本文并不是抒情散文，我并不打算也沒有這個能力去探究那些更深層次的數學理論關系，但了解三角形并不是制作變形金剛，你可以在一張紙上畫三個點然後用直線把他們連起來，那就是一個三角形了如圖（a.1）
（圖a.1）
三角形有非常多的特性，首先，确定它是一個穩定的結構。另外确定一個平面也僅僅隻需要一個三角形就足夠了，三角形的所有内角角度之和是180°（a.2）。
（圖a.2）
但最為有意思的是被稱之為直角三角形的東西，在直角三角形中有一個角的角度是90°（a.3），例如圖（a.3）就是一個直角三角形。
圖（a.3）
實際上在這個直角三角形中，三角形的邊長比例，将随着角度θ的變化而變化，一個角度θ決定了abc三邊長度的比例關系，于是在這裡，我們引入了一個叫三角函數的東西，其中
，當θ小于90°的時候，我們可以在直角三角形中非常直覺地看到三角函數的變化關系，但在直角三角形中，θ的取值範圍，也被限制到了0到90°，為了能讓θ表示更大的範圍，我們就需要引入新的表達方式了。
我們首先先建立一個直角坐标系（a.4），然後以原點為院校，畫一個圓。同時，我們在圓上任意取一點，并将該點的坐标設為（x，y）
（圖a.4）
通過勾股定理我們知道，圓上的點到原點的距離，
那麼，對三角函數我們重新定義為
，這樣一來，我們就可以表示θ為任意實數時三角函數對應的值了。
通過這個圖我們同時可以知道，當
時，實際上相當于在一個圓上繞了若幹個個圈，你可以想象看着你家裡的時鐘，秒針每過60秒就轉了一個圈，你現在看到秒針的位置，在60秒後，120秒後，180秒後…它仍然會指向同樣的位置，這個特性在三角函數中同樣有效，我們管它叫做三角函數的周期性。
現在，讓我們引入一個新的符号π，π在角度上代表180°，除了角度外，我們還引入弧度，它的定義是弧長等于半徑的弧，其所對的圓心角為1弧度。實際上半圓弧長和半徑的比例恰好是一個定值，是以，π除了在角度上表示180°，在弧度上則是一個定值，這個值是一個無限不循環小數，我們常常将它約等于為3.14。
是以
，同時因為周期性，
其中n是一個整數。
那麼，在第一章節，三角函數sin,cos的幾何意義，就顯而易見了。
讓三角函數動起來
假如我們把時間引入進來，那麼，三角函數就開始變得生動了，最直覺的比喻就是現在挂在大廳牆上的時鐘，秒針每分鐘都會轉動360°，現在，讓我們假設秒針指向“12”，也就是垂直于水準面時，它的角度是0°，那麼經過15秒後，秒針指向“3”，也就是轉動了90°，30秒後經過了半分鐘，指向“6”也就是轉動了180°，那麼，秒針每秒轉動的角度是
,我們用t來表示時間，那麼時間與秒針轉動角度我們可以用θ=6t來表示，在實體上，我們常常使用ω來表示角速度，那麼，時間内轉過的角度就是
θ=ωt
現在，讓我們畫一個二維坐标系，并設橫坐标為時間，角速度ω=2π，那麼，
的坐标如下圖（b.1）所示
（圖b.1）
可以看到，在一秒0-1秒，
已經是一個完整的周期了，因為它在1秒内擁有一個完整的周期，我們就稱他的頻率為1hz，顯然的當ω=4π時1秒内有兩個完整的周期，是以，它的頻率就為2hz（圖b.2）。
（圖b.2）
可以看到，頻率實際上和角速度是相關聯的，我們使用字母f來表示頻率，函數公式如下
角速度和頻率是一個正比關系，角速度越大，頻率也就越大
對正餘弦函數而言，我們也常常寫作
三角函數的一些常用公式
我們可以非常直覺的從幾何圖像中推導出三角函數的一些性質，畫一個坐标系，同時以原點為圓心繪制一個半徑為1的圓（圖c.1）
（圖c.1）
那麼
可知
由周期性可知
同時，三角函數間是可以互相轉換的，我們很容易得出這樣的公式：
因為正餘弦函數存在這種互相轉換的關系，是以之後我們統稱它們為正弦函數
我們還可以進一步證明二角的更多公式（證明過程略），這些公式我們都可以直接拿來使用
二角和差：
和差化積
通過聯立二角和差公式，我們還可以得到積化和差公式
和二角和差推導出的二倍角公式
三角函數的正交性
在讨論正交性之前，讓我們用最簡單的了解下微積分，我們取一個最簡單的一進制函數做比方
現在，讓我們畫一個坐标軸，那麼，這個函數的圖像是這樣的(d.1)：
(d.1)
現在，做一條經過（2,0）并垂直于x軸的直線，交
于點A，（2,0）為B，原點為C如圖（d.2）
圖d.2
那麼，我們很容易求出三角形ABC的面積
實際上這個求面積的過程，就是微積分于該函數上的幾何意義，這個求面積的過程，我們用微積分來表示，就是
這個積分方程實際上是求四邊形ABCD的面積（d.3）
圖d.3
在二維平面上我們很容易将微積分了解為函數在一定範圍内和x軸圍成的面積，在但實際上面積和微積分稍有不同，因為面積肯定是一個非負數，但積分卻可以是負數，
他在坐标系中是三角形HCI的面積，但是它的x軸下方，我們可以了解為三角形的高度是一個“負數”，是以，這個區域的面積也是一個“負的面積”，是以，這段的積分為-2（d.3）
圖d.3
實際上連續的中心對稱函數圖形h(x)，在
的積分都是0
顯然的，正弦函數圖像滿足這個性質（d.4），不僅如此，我們可以非常直覺的看出來，正餘弦函數在函數的一個周期内的積分都是0（d.4 d.5）
圖d.4
圖d.5
現在，我們假設有兩個自然數m，n且
，并假設兩個函數
然後我們将這兩個正餘弦函數任意進行組合，并對他們在-π到π進行積分
通過上述的常用公式，我們可以推導出下面的三個式子
可以看到，m不等于n時，積分結果都是0
從正交性得到的啟發
還記得之前我們使用的将時間作為變量的三角函數麼
也就是正弦波了，我們用更加通用的公式來表示一個正弦波
或者，餘弦波也同樣是這個道理，畢竟它與正弦波的不同僅僅隻是“移動了一下”
其中，t表示時間，ω表示角速度，如果我們将ω除以一個2π，那麼結果就是頻率也就是f，
我們稱為相位，
表示初相位，A表示振幅，當然，振幅和相位并不影響其正交性，從上章節的正交證明推導出的
我們很容易證明這點。
現在我們可以想象一個正弦波在二維坐标上的樣子，或者是我們在科幻電影或科研實驗室中，看到儀器儀表上那一段段的波形。或者更加直覺的，我們将一塊石投進水池裡，蕩起的波浪也像極了正弦波。
當然，在自然界一般不會出現标準的正弦波，各種波的疊加，你可以想象在一個房間裡一群朋友聊天，或者是在KTV中歌聲和談話聲混雜的場景，它就是聲波的各種疊加
盡管聊天中大家都在說話，但是我們仍然不會把朋友們說的話混淆起來，因為不同的人發出的聲音頻率也不一樣，這也就是我們常說的未見其人先聞其聲，我們天生具備有分别不同頻率聲波的能力，是以盡管環境嘈雜，我們仍然能夠取得我們想要的資訊。
但是現在我們如何使用數學的手段，檢視波形函數f(x)中是否有我們想要頻率的波形呢。
顯而易見的，答案當然就是本章節的标題，應用波的相關性，例如我們使用一個正弦波
去乘以一個由多個正餘弦波疊加而成的函數f(t)，然後對它們進行積分
從相關性可以知道，當角速度（上一章節的m，n）不同時，也就是頻率不同時，積分結果是0，隻有相同時，積分結果才不為0，是以，如果f(t)中包含1HZ的正弦波，那麼積分的結果就不為0，否者，結果就是0
那麼，一個檢波手段也就誕生了。
歐拉公式複變函數
正弦函數和複數表面上并沒有什麼關聯，但正所謂千裡姻緣一線牽，在我們考慮着那根看不見的紅線另一端系的是誰的時候，歐拉早在幾個世紀前就将正弦函數與複數牽了一條紅線，進而撮合了數學上這一段流傳千古的因緣，不過如果再在紅線上扯下去，這篇文章就要變愛情小說了，但是我們仍然需要回到這個渣作的現實三次元，現在讓我們來介紹一下大名鼎鼎的一個複變函數，它的知名程度基本上在是數學上的安徒生童話，人人皆知
它的公式如下
其中，i表示複數的虛部，它是
，也就是說
，當然我們并不能深究他在現實生活中的意義，但它在數學多個方面，起着舉足輕重的意義。
那麼，複數如何了解呢
我們來看看複數的标準形式
a + bi
其中，a，b為任意實數，a在複數中表示實部，b在複數中表示虛部
假設我們現在畫一條二維坐标系，x軸為實部，y表示虛部，那麼，1+2i實際上表示的是坐标上的（1,2）
圖d.6
那麼
在坐标系上，實際上是一個半徑為1的圓（d.7）
圖d.7
複數在信号的表示
現在讓我們來考慮下面這個複數
它的實部為
，如果我們希望用
的形式來表示它，那麼，它就變成了
畫在坐标系中，如圖d.8
圖d.8
按照極坐标的角度而言，即一個長度為2的變繞X軸逆時針旋轉
現在讓我們更進一步放大我們的腦洞，以原點為圓心，2為半徑畫一個圓（圖d.9）
圖d.9
那麼,假如我們将
看作是與原點距離為1，繞x軸逆
時針旋轉的點，将
與之相乘
就變成了，初始位置為
以角速度以原點為圓心逆時針旋轉的點
那麼，假如我們将複數用在信号中，就從複數坐标系上的點拓展到了線上信号的幅度與初相了。
還記得之前我們寫的通用的正弦/餘弦函數麼，
這就是上面的通用公式，因為這個公式太重要了（或者說避免你翻回去找這個公式），我将它再寫了一遍
可以看到，對于一個正選函數，我們僅需要知道它的幅度，角頻率，初相，就可以确定這個正弦函數是什麼樣子的了
但正弦函數比善變的女人還有能耐，通過一些數學變形，你還會發現它有時比哄妹紙有意思多了，通過三角函數的公式，我們可以得到
設
那麼就有
那麼，這個函數就變成了正餘弦函數的組合，同時我們也得出了
那麼就有
=
=
傅裡葉的故事
這裡我們抛開那些繁瑣的數學公式，然後把時間拉回到17世紀，當時有一位法國的男爵名叫巴普蒂斯·約瑟夫·傅裡葉（Baron Jean Baptiste Joseph Fourier）當然，并不是因為它當了一官半職我們才介紹它，但他的成就，卻從17世紀一直影響到我們今天，毫不誇張的說他奠定了信号系統的基礎(或者說給出了指導方向？)。
具體時間還是要回到18世紀，伯努利（D.Bernoulli）就曾經提出：一個弦的實際運動都可以用标準張模的線性組合來表示，但是當時另一個數學大神拉格朗日是不茲磁這個說法的，拉格朗日認為，不可能用三角級數來表示一個具有間斷點的函數，就在這個環境下，傅裡葉仍然堅信：“任何”周期信号都能夠成諧波關系的正弦級數來表示，他還專門寫了一篇論文，但迫于拉格朗日當時在學術界的威望，傅裡葉理論的論文一直沒能被發表，直到傅裡葉的晚年，他才得到他應有的承認。
當然，以我們現在的角度來說，傅裡葉當時的理論是有缺陷的，在後人不懈努力下，傅裡葉變換才趨于完善，實際上傅裡葉變換并不是因為傅裡葉對這個理論在數學上做了多大的貢獻，但傅裡葉當時對問題的前瞻性與指導性，卻奠定了這個數字信号舉足輕重的公式基礎。
檢波與傅裡葉變換
現在我們讨論的就是傅裡葉大神“任何連續周期信号都可以由一組适當的正弦曲線組合而成”這一命題，在信号系統或者是工學，傅裡葉變換絕對是入門的加減乘除，實際上傅裡葉變換并不是一個公式，而是多個分别應對不同類型信号的多個公式，但萬變不離其宗，不管變體如何變換，其核心的思想是不會改變的。
那麼傅裡葉變換的意義是什麼呢，回想一下我們之前說的正弦函數的相關性，通過這個相關性，我們可以檢測某個波裡是否包含某個頻率的正弦波，進一步的，假設我們用無限多個不同頻率的正弦波與它正交組合，我們就能知道，這個波是由哪些頻率的正弦波組合而成的了，我們管波在時間軸上的表達，叫做時域，圖b.1其實就是
的時域圖，通過傅裡葉變換，我們能夠将
是由哪些頻率組成的圖像表示出來（當然隻有一個1HZ頻率），也就是橫坐标由時間變成了頻率，也就是我們說的頻域了。
簡單來說，傅裡葉變換就是将時域信号轉換為頻域信号的一個變換過程
這裡，我們先來看看連續傅裡葉變換，因為傅裡葉英文的開頭是F，是以我們使用大寫的F來表示傅裡葉變換，f(t)用來表示某連續非周期性的時域信号函數
看看我們的歐拉公式
然後把歐拉公式代入傅裡葉變換
角速度乘以時間，不就是
了麼，我們将這個變換函數寫成這種形式
結果變得顯而易見了，簡單來看，傅裡葉變換與檢波的手段多多少少的相似性
傅裡葉級數的三角函數表達
但僅從上面的檢波手法來推斷傅裡葉變換，是不嚴謹的，
設一個函數由一個直流分量（簡單來說就是一個常數）和多個正弦函數組成，那麼它可以寫成這種形式
其中
表示一個累加符号，
相當于0+1+2+3…..+N-1。因為定積分隻能用線上性可導的函數中，對于離散的信号，我們隻能用累加。如果你學過C或java，C#之類的語言，那麼它很像這段代碼
int sum;
for(n=0;n<N-1;n++)
{
Sum+=n;
}
并且由上式可知，
表示這個三角函數的角速度或者叫角頻率，當n=1時，我們管
叫基頻，管
叫傅裡葉級數（餘弦信号形式）
利用三角函數的變換公式，上式可變形為
設
那麼，上式變為
現在，讓我們正式的引入正交性的性質，還記得檢波手段麼，這裡，我們假設對f(x)用
進行檢波，那麼就有
假設f(x)中含有
角頻率的正弦波系數為，那麼根據三角函數的正交性，上式就有
進一步計算，可得
同樣，
也可以使用相同的方式進行推導
是以，通過
也就是其傅裡葉級數，我們可以知道這個波的波幅與相位：
周期連續時間傅裡葉級數
現在，讓我們來想象一個函數f(x)它是一個周期
函數，那麼根據傅裡葉的理論，它能夠表示成若幹個（無窮個）一組“适當”的正弦曲線組合而成，在前面幾個章節，我們通過歐拉公式，得到
顯然，它是一個周期性函數，當然i是一個複數，顯然他是周期複指數信号，同時因為那麼，這個公式也可以寫成
或者是
顯然的，這個複指數信号的頻率是
，現在我們假設有另一組“适當”的正弦函數，它們的頻率剛好是
的整數倍并且幅度也不同，那麼，這組信号我們可以使用
來表示（k為自然數）。那麼從上面的根據歐拉公式，我們也很容易得出下面的推導
現在将這兩個公式帶入
得到
化簡後得
是以，上式最終變為了
設
上式就寫成了,n=k
顯然，一個信号如果可以寫成這種形式，那麼
就是對應不同頻率正弦波的系數（複數形式），我們稱之為頻譜系數，那麼，如何通過f(x)求的
呢
首先，我們将這個等式的兩邊分别乘以
并對兩邊同時在一個周期内積分，那麼我們就得到公式
進一步變形為
于是，得到
也就是
實際上$$a_{n}$$就是我們所說的傅裡葉級數，或者說是頻域系數。通過這個系數的值，我們可以知道這個頻率的波對原始波的能量貢獻值，在之前的《信号的複數表示》章節中，我們可以了解到這個系數确定了該頻率波的幅度，初相，進而完成信号時域到頻域的分解，并且我們還知道了
也就是說，通過
的模
我們知道對應角頻率波的波幅（等于波幅的一半）
通過
可以得出其相角。
周期離散時間傅裡葉級數
在上一個章節中，讨論了周期連續時間信号的傅裡葉級數求解方式，那麼，連續信号可以求其級數，離散的是否也有這樣一個公式呢
但在介紹離散變換變換前，我們先來了解連續和離散是什麼，。其實顧名思義。比如給你一段長度100米的繩子，當然，這個100米的繩子是連續的，如果你在50米的地方剪短這根繩子，那麼它就是不連續的了，假設我們把這根繩子切成若幹個段。直到每個段都變成一個“點”，我們可以直接用數字編号每一個點，那麼他就變成離散的了。
用圖像來繼續說明，如圖e.1，這是一個sinx的函數圖像，當然，它是周期無限長且連續的
圖e.1
現在，我們每隔
就取一個點，那麼這些點在坐标系上就是周期離散的（e.2）
圖e.2
說到這裡，我們現在要介紹的變換，就是處理這些點的變換方式
首先，因為是等間距對周期信号采樣，是以采樣出的點也是周期性的，假設采樣的點用數組x[n]來表示，也就是說假設周期為N，x[a]與x[a+Nk]是完全相等的值，也就是說，整個離散樣本中取任意周期N内他們的累加和是一樣的，同樣的，x[n]仍然可以用“恰當”的正弦函數組合而成（或者說可以由基波頻率為
的一系列波形組合而成）基于這點，我們就可以将x[n]寫成這種形式（k=<N>表示取任意連續N點即一個周期内點，不管如何取，結果都是一樣的）
那麼，
就是我們要求的周期離散時間傅裡葉級數了，它的推導與上一章節周期連續時間傅裡葉級數的求解方式出奇的相似
首先，因為周期為N，是以
，那麼，式子就變形為
我們取
可以知道，當K不等于N時的結果為0，僅當K等于N時的結果為N。同樣的，我們再次将兩邊同時乘以
得到
然後再同時對兩邊進行N項上求和
顯然的，僅當k=r或n為N的整數倍時，右式才不為0，那麼，上式就變為
那麼
頻譜系數求得！。
非周期離散時間傅裡葉變換
如果嚴格來說，自然界大多是沒有理想的周期信号的，那麼，是否有辦法處理非周期離散時間的信号麼。
我們來看看這樣一個離散信号（圖e.3），它隻有一個脈沖取樣
圖e.3
這個脈沖信号僅在-3，-2，-1上是有值的，其餘的值都是0
那麼我們是否可以把它當成一個周期無窮大的周期信号呢。
我們先來看看上一節中周期離散時間傅裡葉級數的分析公式
假如把這個分析公式套用在上述的信号中，那麼因為僅在三點有值（其它都是0）并且周期是無窮大那麼我們就得到了公式
當然，它和下面這個式子是等價的
這樣我們就将公式推廣到了更加通用的公式類型
現在定義函數
那麼
我們現在再将$$a_{r}$$代入原式
中，得到
因為
是以
從式中看出，随着N趨近于無窮大，
趨近于無窮小，那麼，上式就從累加變成了積分，且因為
的周期為2π，且其僅在周期内有值，于是，上式也随之變為了
非周期離散有限長度傅裡葉變換
最後，我們将上述的變換公式進行進一步的推廣，就是非周期離散有限長度傅裡葉變換了，實際上它與周期離散傅裡葉變換已經非常的接近了
如果我們将有限長的信号推廣到無限長的信号，那麼我們先假設信号的樣本點數為N個，那麼，信号
的n取值範圍就可以定義在[0,N-1]
我們假設将這個有限長的區間補到無限長，除了[0,N-1]區間，我們僅僅需要在其他區間再補上N個同樣的離散信号就行了，這并不影響其結果，那麼我們就可以将有限長非周期離散信号變為周期離散信号了，這樣我們就可以直接套用周期離散時間傅裡葉變換的分析公式：
為了友善計算，我們周期取[0,N-1]，那麼，公式就變成了
當然在實際應用中，我們常常設：
那麼就有了有限長非周期傅裡葉變換的的分析公式：
又因為
将
代入回周期離散時間傅裡葉變換的綜合公式中，得到：
為了友善計算，我們仍然取[0,N-1]作為計算周期那麼就得到了有限長離散時間傅裡葉變換的綜合公式
非周期離散時間傅裡葉變換的應用
在上面幾個章節中，我們從周期連續時間的傅裡葉級數逐漸的推廣到了有限長離散時間傅裡葉變換。雖然内容不少，但是真正實際在日常用的多的，是有限長非周期離散時間傅裡葉變換，為何？當然我們幸運的不是在一個老牛拉破車的年代計算隻能靠人腦算盤，現在的信号處理除非一些數學推導應用，大多數實際生活應用都是靠計算機來完成的，而這也決定了我們的公式必須與計算機的硬體相關，我們的計算機目前存儲空間是有限的，而連續的信号（不管是頻域還是時域），就相當于由無窮多個點組成，很遺憾，現在仍然沒有存儲容量無窮大的記憶體，就算有，也沒有能夠處理無窮大資料的CPU，是以，我們無法直接處理連續的信号，觀察上面幾個變換，也就隻有時域和頻域都是有限長度且離散的有限長離散時域傅裡葉變換能夠滿足我們的需求了。
下面為了說明友善，我們将有限長度且離散的有限長離散時域傅裡葉變換的分析公式統一簡稱為離散傅裡葉變換（Discrete Fourier Transform）或者其英文縮寫DFT，将它的逆變換（Inverse Discrete Fourier Transform）也就是綜合公式簡稱為IDFT。
現在讓我們來看看離散傅裡葉變換對
其中，
為采樣率,N表示信号的采樣點數，n的範圍是0到N-1，
的第一個點，
的第二個點，以此類推，k的取值範圍是0到N-1，由
可知，其分辨率等于
為了更加友善淺顯地了解其中的計算，我們繼續觀察公式中的
仍然是我們熟悉的味道，現在我們使用歐拉公式替換它于是離散傅裡葉變換的公式變成了
可以看到，它實際上是所有樣本點的累加和，那麼它意味着什麼呢，我想象一個波形函數
假設當中的正弦波最大頻率為F并且頻率都是自然數, 現在在紙上你可以随手畫一個正弦波，你可以看看你至少需要采樣幾個點才能夠還原出這個波出來，顯然的，我們可以在波峰波谷分别取樣一個點，這樣我們就能靠這兩個點将它還原了，這也就是說，我們至少需要2F的采樣頻率，才不至于讓它失真，實際上，采樣頻率至少要被采樣頻率的2倍才不至于失真，這個采樣定律我們管它叫香農采樣定律，現在假設我們就使用2F的頻率對g(x)進行等距采樣，并采樣了1秒，顯然的，采樣後的結果每個波都進行了周期整數倍的采樣，因為正弦波在一個周期内的累加是0，是以累加的結果也會是0，是以，将所有樣本點加起來的和，實際上就是$$beta$$的N倍，我們将累加和除以N就能得到
，實際上我們管
叫做直流分量。
當k=1時，也就是
或者寫作
它表示基頻成分，因為我們從公式中可以看到，頻率都是其整數倍，至于多少倍，就是由k值決定的。
現在，讓我們用一個實際的範例來驗證離散傅裡葉變換
假設正弦函數
我們假設對這個函數進行采樣，采樣頻率是4HZ，那麼，實際上我們将在0.25s，0.5s，0.75s，1s處取得其樣本點，那麼，對應的值應該是
現在對其進行離散傅裡葉變換那麼
因為有四個樣本點，是以，N的值是4，k的取值範圍是0-N-1，也就是0到3了
我們先來計算
的值
這點沒錯，因為
沒有直流分量，是以它理所當然是0
繼續是
的值
可以看到，結果的實部是有值的，那麼它對應的頻率是多少6呢，以一個不大準确的運算來講，首先讓我們回到變換公式中的
，這個函數的頻率是多少呢，實際上，n就是我們在時域的t，那麼
，基波頻率
，因為k和采樣頻率是對應的，是以，基波頻率f實際上就是采樣頻率除以采樣點數。
也确實是1HZ，那麼，結果中實部的2是怎麼回事，2與這個波的振幅有什麼關系。
從運算的過程我們可以看到，
進行了相關操作，實際就是
，我們采樣的點是波峰與波谷，然後将所有值累加，波峰是1，波谷是-1，于是平方後都變成了正數，累加後就是0+1+0+1=2了
這是2HZ的值，顯然，它為0
最後是
的值是2，但是我們知道，并不包含3HZ的波，實際上根據香農采樣定律，4HZ的采樣隻能表達2HZ的波，是以這個點實際上是不準确的。但是，出現2的結果并不是偶然，我們接着往下看。
巧妙的對稱性，共轭
如果說對稱是世界上共有的一種美的表達，那麼，在幾何平面上，無窮多的函數就擁有這種對稱性，在對稱美這一個方面，出色的數學家也許并不遜色于畢加索。
當然深究的話就是後話了，在這裡我們來簡單看看幾個擁有對稱性的簡單函數：
圖f.2
這是一個非常簡單的二進制一次方程，可以從圖f.2中看到，他是關于y軸對稱的，這句話如果用數學的語言來将，非常簡單且直覺的，我們可以用下面的公式來表達這個“對稱”的思想
對于這類關于滿足上述公式的函數，我們管它叫偶函數。
現在再讓我們看另外一個函數
它的函數圖像如圖f.3
圖f.3
可以看到，函數的圖像是原點對稱的，相對的，我們用如下公式表示這種“原點對稱”關系
我們管它叫奇對稱
當然，對稱未必一定要是y軸或者是原點，現在我們将
的函數圖像向右平移兩個機關，變成
這樣，它就關于x=2這條豎線對稱了（f.4）
圖f.4
這回我們在數學上用這個方程式來表示它關于x=2對稱
最後，我們用更加通用的方程式來表示某個函數f(x)關于x=N/2這條垂線對稱
現在，讓我們來看看共轭是什麼：
你看過幾米的《向左走，向右走》麼，共轭是數學中一個文藝而浪漫的代名詞，一句相當文藝詩來總結這個關系，就是：
向左走，向右走，縱使背道而馳，相隔萬裡，但隻要心連接配接在一起，也終會在大洋的彼岸，迎來相會的交點。
文藝的詩歌藝術生看到後也許就開始感歎，多麼美的詩句，世間萬物芸芸衆生，千裡有緣一線牽，感謝我能遇見你，但對于地理大神也許不削一顧：這不就是說地球是圓的麼，數學系的牛人毅然站出來，别做夢了少年，你們隻會越離越遠，就算到了世界末日，你們也碰不到一塊兒。
結果，文藝青年赢的了女神的芳心。理工大神則斬獲了“屌絲”的稱号
實際上，共轭我們可以了解成相關聯，也就是所謂的“緣分”，就是這一對數存在某種聯系的意思，這裡我們不深究更深層次的含義，我們主要說說共轭複數，首先，複數的形式是：
a + bi
它的共轭是
a - bi
非常簡單的一句概括：實部相等，虛部取反，這兩個複數互為共轭。
那麼，和我們之前說的對稱性相結合的話什麼是共轭對稱性呢
我們這樣給出定義：
當一個函數f其實部為偶函數，虛部為奇函數時，此函數就為共轭對稱函數，即f(x)的共轭等于f(-x)，舉一個非常簡單的例子：
顯然的，實部是個
偶函數，虛部xi是個奇函數，是以它是一個共轭對稱函數，現在，我們來看看離散的範例，假設有以下複數序列
我們分别提取實部與虛部那麼分别是
實部：1,2,3,3,2,1
虛部：1,2,3，-3，-2，-1
我們很清晰的看到了序列的對稱性，實部偶對稱，虛部奇對稱。在這裡，因為他是離散的序列，是以我們管它叫共轭對稱序列。
離散傅裡葉變換的共轭對稱性
對稱之美存在世間萬物中的每一個角落，作為信号與系統中最美的變換函數之一，她也存在這種對稱性。
我們回到離散傅裡葉變換的公式：
我們假設：
依據歐拉公式，可以得出
那麼，上式可以寫成
我們設k=N-k，并将它帶入
當中，那麼就有
進行對比，發現，實部相同虛部相反，假如輸入的信号為實信号時，剛好呈共轭對稱性，将它代入離散傅裡葉變換方程後
其中*的意思就是共轭。也就是說離散傅裡葉變換具有這種共轭對稱性（輸入為實信号時）。
也就是之前所說的直流分量了。
的序列離散傅裡葉變換的結果
根據共轭對稱性
也就是2+0i=2-0i，這也解釋了為什麼
快速傅裡葉變換
在上文的範例中，我們僅僅是計算了序列長度隻有4的離散傅裡葉變換，也可以看到非常麻煩，倘若序列再長點的話，工作量就會呈指數級增長，我們使用矩陣運算來表達離散傅裡葉變換的運算過程
可以看到需要16次乘法運算，按照時間複雜度來算的話，它的複雜度是
，當然，這還沒算上sin cos與複數加減法帶來的性能開銷。是以，如果你想将一段2000長度的離散信号進行DFT運算，意味着你至少需要做400萬次的運算，巨大的性能花銷，導緻了DFT在以前并不被看好，畢竟除非一些非常重要的信号，誰會花大量的人力物力去做這幾百萬次的運算呢，即使是在今天，幾百萬次的計算開銷在有了計算機的幫助之後，也并不算一個小數目，倘若對于那些實時的頻域分析，這些計算開銷都是非常昂貴的，不過幸運的是，一個DFT的優化算法很快被開發出來，我們稱之為快速傅裡葉變換（Fast Fourier Transformation），當然，其本質上仍然是DFT隻是算法上進行了優化使它更加适合于計算機處理，不過話說歸說，不給出理論證明的廣告都是瞎扯淡，那麼，接下來就繼續看看，DFT是如何被優化的：
首先第一點，我們仍然先把離散傅裡葉變換對貼上來
現在，我們假設對離散序列x[n]的離散傅裡葉變換寫作
其中，k，n為範圍為[0,N-1]的整數，現在，我們将離散序列x[n]分成兩組，例如x[0],x[2],x[4],x[8]……x[2i]為偶數序列，記為x[2i],意思是偶（even）序列，而将x[1],x[3],x[5]….x[2i+1]記為x[2i+1],意思是奇(odd)序列
那麼，不管是偶序列還是奇序列，他們的個數都變為了
個，是以，i的取值範圍是[0，
],同時為了剛好平分，也要求N必須是一個偶數。
那麼
就可以寫成
進一步變形，得到
還記得
麼
那麼
就是
于是，上式又可以寫成
後面兩個式子是不是非常眼熟，沒錯，一個DFT變換變成了2個DFT變換，不同的是，後面序列的長度隻有前面的一半，為了保證後面兩個DFT變換成立，k的範圍也應随之變為了
既然，前半部分可以變為兩個DFT變換，那麼後半部分呢，我們使用
來表示離散序列後半部分，那麼，DFT就變為了
進一步變形得到
根據歐拉公式，因為
是以
可以看到，僅僅是多了一個
，其它都與之前的推導相同，那麼這也就是為什麼我們要将離散序列分為奇數列與偶數列的原因了，偶數列不變，奇數列多個負号，于是，公式變為了
那麼，公式總結為
那麼，計算的複雜度就從
變為了
隻要n的值大于2，計算的時間複雜度無疑是降低了，我們進一步對多項式進行分解，直到最後僅剩下2個離散點的傅裡葉變換，那麼，它的複雜度将會降低至
，同時，這也就要求我們的離散序列項的個數必須是2的整數次幂，是以，這種快速傅裡葉變換又叫做基2快速傅裡葉變換，當然同樣的，也有其它基底的快速傅裡葉變換，但這裡就不做更多的讨論了。
快速傅裡葉逆變換
快速傅裡葉逆變換完全可按照正變換的過程進行推導，實際上就是換湯不換藥，根據IDFT公式
同樣的，我們将離散序列X[n]分成兩組，例如X[0],X[2],X[4],x[6]……x[2i]為偶數序列，記為x[2i],，而将X[1],X[3],X[5]….X[2i+1]記為X[2i+1],意思是奇(odd)序列
那麼，公式可以寫成
進一步變形，得到
後式與DFT正變換，僅僅隻是W指數正負的不同，我們隻需要修改下符号就可以了。
按照上面步驟同樣推導出（因推導步驟一樣，推導過程略）：
使用C語言編寫DFT/IDFT代碼
回到離散傅裡葉變換對
我們需要先将它們變成容易編碼的格式，首先是正變換
利用歐拉公式，變為
然後是逆變換
利用歐拉公式，變形為
首先因為頻域涉及到複數運算（輸入信号一般為實信号）是以我們需要複數結構體complex,現定義結構體
typedef struct __complex
{
float re;// really
float im;// imaginary
}complex;
及複數的加乘運算函數
complex complexAdd(complex a,complex b)
{
complex ret;
ret.re=a.re+b.re;
ret.im=a.im+b.im;
return ret;
}
complex complexMult(complex a,complex b)
{
complex ret;
ret.re=a.re*b.re-a.im*b.im;
ret.im=a.im*b.re+a.re*b.im;
return ret;
}
之後定義DFT運算函數
void DFT(complex x[],complex X[],int N);
其中，x[]為輸入時域信号，X[]為輸出頻域信号，N為時域信号的樣本點數
void DFT(_IN complex x[],_OUT complex X[],int N)
{
int k,n;
complex Wnk;
for (k=0;k<N;k++)
{
X[k].re=0;
X[k].im=0;
for (n=0;n<N;n++)
{
Wnk.re=(float)cos(2*__PI*k*n/N);
Wnk.im=(float)-sin(2*__PI*k*n/N);
X[k]=complexAdd(X[k],complexMult(x[n],Wnk));
}
}
}
同時定義DFT逆變換函數
void IDFT(complex x[],complex X[],int N);
其中，X[]為輸入頻域信号，x[]為輸出時域信号，N為頻域信号的樣本點數
void IDFT(_IN complex X[],_OUT complex x[],int N)
{
int k,n;
float im=0;
complex ejw;
for (k=0;k<N;k++)
{
x[k].re=0;
x[k].im=0;
for (n=0;n<N;n++)
{
ejw.re=(float)cos(2*__PI*k*n/N);
ejw.im=(float)sin(2*__PI*k*n/N);
x[k]=complexAdd(x[k],complexMult(X[n],ejw));
}
x[k].re/=N;x[k].im/=N;
}
}
實際上通過離散傅裡葉變換後頻域的共轭對稱性，我們僅僅需要計算前半部分就可以了，這個優化操作若讀者有興趣可以自己完成。
使用C語言編寫FFT/IFFT代碼
根據FFT的公式有
其中，$$W_{N}^{k}$$寫成易于編碼的格式為
那麼，我們就可以采用疊代的方式，将輸入序列不斷地拆分重組，直到它變為2點的DFT，代碼如下：
void FFT_Base2(_IN _OUT complex x[],int N)
{
int exbase,exrang,i,j,k;
complex excomplex,Wnk,cx0,cx1;
if (N>>2)
{
// x[] 4 base odd/even Sort
exbase=1;
exrang=0;
while (exrang<N)
{
exrang=exbase<<2;
for (i=0;i<N/exrang;i++)//for each token
{
for (j=0;j<exbase;j++)//for each atom in token
{
excomplex=x[exrang*i+exbase+j];
x[exrang*i+exbase+j]=x[exrang*i+exbase*2+j];
x[exrang*i+exbase*2+j]=excomplex;
}
}
exbase<<=1;
}
FFT_Base2(x,N>>1);
FFT_Base2(x+(N>>1),N>>1);
for(k=0;k<N>>1;k++)
{
Wnk.re=(float)cos(-2*__PI*k/N);
Wnk.im=(float)sin(-2*__PI*k/N);
cx0=x[k];
cx1=x[k+(N>>1)];
x[k]=complexAdd(cx0,complexMult(Wnk,cx1));
Wnk.re=-Wnk.re;
Wnk.im=-Wnk.im;
x[k+(N>>1)]=complexAdd(cx0,complexMult(Wnk,cx1));
}
}
else
{
//2 dot DFT
cx0=x[0];
cx1=x[1];
x[0]=complexAdd(cx0,cx1);
cx1.im=-cx1.im;
cx1.re=-cx1.re;
x[1]=complexAdd(cx0,cx1);
}
}
void FFT(_IN complex x[],_OUT complex X[],int N)
{
int size=1;
complex *i_px;
while((size<<=1)<N);
i_px=(complex *)malloc(sizeof(complex)*size);
memset(i_px,0,sizeof(complex)*N);
memcpy(i_px,x,sizeof(complex)*N);
FFT_Base2(i_px,size);
memcpy(X,i_px,sizeof(complex)*N);
free(i_px);
}
IFFT的理論則依據下面兩個公式，就不再複述了：
C語言代碼如下：
void IFFT_Base2(_IN _OUT complex X[],int N)
{
int exbase,exrang,i,j,n;
complex excomplex,Wnnk,cx0,cx1;
if (N>>2)
{
// x[] 4 base odd/even Sort
exbase=1;
exrang=0;
while (exrang<N)
{
exrang=exbase<<2;
for (i=0;i<N/exrang;i++)//for each token
{
for (j=0;j<exbase;j++)//for each atom in token
{
excomplex=X[exrang*i+exbase+j];
X[exrang*i+exbase+j]=X[exrang*i+exbase*2+j];
X[exrang*i+exbase*2+j]=excomplex;
}
}
exbase<<=1;
}
IFFT_Base2(X,N>>1);
IFFT_Base2(X+(N>>1),N>>1);
for(n=0;n<N>>1;n++)
{
Wnnk.re=(float)cos(2*__PI*n/N);
Wnnk.im=(float)sin(2*__PI*n/N);
cx0=X[n];
cx1=X[n+(N>>1)];
X[n]=complexAdd(cx0,complexMult(Wnnk,cx1));
Wnnk.re=-Wnnk.re;
Wnnk.im=-Wnnk.im;
X[n+(N>>1)]=complexAdd(cx0,complexMult(Wnnk,cx1));
}
}
else
{
//2 dot IDFT
cx0=X[0];
cx1=X[1];
X[0]=complexAdd(cx0,cx1);
cx1.im=-cx1.im;
cx1.re=-cx1.re;
X[1]=complexAdd(cx0,cx1);
}
}
void IFFT(_IN complex X[],_OUT complex x[],int N)
{
int size=1,i;
complex *i_px;
while((size<<=1)<N);
i_px=(complex *)malloc(sizeof(complex)*size);
memset(i_px,0,sizeof(complex)*N);
memcpy(i_px,X,sizeof(complex)*N);
IFFT_Base2(i_px,size);
// 1/N operate
for (i=0;i<N;i++)
{
i_px[i].re/=N;
i_px[i].im/=N;
}
memcpy(x,i_px,sizeof(complex)*N);
free(i_px);
}
信号與采樣
《淮南子·說山訓》中有“見一葉落而知歲之将暮。”, 宋朝《文錄》則引曰：“山僧不解數甲子，一葉落知天下秋。”直到今天的“一沙一世界,一花一天堂.雙手握無限,刹那是永恒”都在傳達着局部概括全局，一刻即是永恒的思想。
如果是一個充滿文藝氣息的理科青年，應該很容易就能發出這樣的感慨：
我們都是宇宙的一份子，我們是宇宙的縮影，即便微不足道，但冥冥之中我們也是世界不可或缺的一環。
但如果是一個理科大神，這會兒可要費點腦子了，大神抱着心愛的四路泰坦32G記憶體1T的PCIE SSD和ryzen 18x的電腦，琢磨着如何把
的波形圖像完整地存入電腦，很快的大神發現即使它再把電腦容量更新一倍，他也無法存儲sinx的完整波形，哪怕是一段他也辦不到:
首先sin(x)是周期函數，它的區間是
,很遺憾，就算把整個地球的沙子都做成記憶體顆粒，也沒有辦法存儲一個無窮大的量。
那麼存儲一個周期如何，很遺憾，sin(x)是連續的函數，就算是一個周期也包含着無窮多個幅度資訊，就算把整個銀河系的沙子都做成記憶體顆粒，也無法存儲無窮多的幅度資訊。
不過很快的，大神找到了門道，既然無法存儲完整的波形，那麼存儲當中一些關鍵的點總能辦得到
很快，大神就找到了一個周期内的波峰和波谷,如圖x.1
圖x.1
很快，波峰和波谷的水準間距剛好是周期的一半，而波峰與x軸的垂直距離剛好就是波幅。那麼最終大神用2個離散的點表示了sin(x)，并且我們也知道了，假設sin(x)的頻率是n，那麼我們至少要以2n的頻率進行采樣才能夠還原出原始的波形。
那麼，波形函數如何以數學的形式對一個時域點進行采樣呢，實際上前輩們早就定義了一個理想的函數
而非無窮大以便于我們的處理
是以我們就可以得到一個采樣的方程，例如對sin(x)的x=0點進行采樣那麼我們就可以用
進行右移位處理就可以了
當然假設我們需要對多個點進行采樣，那麼我們完全可以寫成這種形式(對sinx整數值采樣)
在奧本海姆的《信号與系統》7.1.1中有如下的推論
設對一信号以T間隔取樣，那麼就有如下數學式
根據卷積性質時域内的相乘對應于頻域内的卷積，那麼就有
為沖擊串頻域函數
根據公式（《信号與系統》例4.8推論,
為基波頻率或叫頻域間距）
及（*為卷積符号）
于是有
上式說明
是頻率為
的周期函數，由一組移位的X(jw)疊加而成，在幅度标以
的變化，當
大于頻帶寬度的一半時
頻帶不會發生堆疊，反之會發生堆疊（引用信号與系統 7.1.1 沖擊串采樣，詳細推論請查閱書籍）。
實際上上面的推論可以用一個簡單的話來說：
表示一個持續期為T且最高頻率為W的時間函數，有2TW的樣本個數就足夠了。
實際上這句話說的内容就是著名的香農采樣定律，看！那些實體學家，數學家和詩人都是雖然表達方式不一樣，但他們表達的東西常常都是一個意思，不管怎麼說，他們都是會玩的。
鋼琴與按鍵識别
錄音頻率與采樣
“未見其人先聞其聲”說的就是靠聲音識别人的道理，從聲學的角度來說，各個人發出的嗓音也是各有不同的，我們大腦進化出了自動篩分頻率的功能，是以盡管我們見不到真人，仍然可以從他說話的聲音分辨出他，在一個嘈雜的環境中，我們也可以在分辨出那個人在說什麼，看來，人腦也是一個實時的濾波系統。
當然，人并不是能聽到所有的聲音，聲音歸根結底仍然是波在媒體中傳播，人隻能聽到20~20000HZ的聲波，是以，低于20HZ的波我們叫次聲波，高于20000HZ的波叫超音波，根據上一章節所屬的香農采樣定律，如果我們要記錄一段聲波例如演唱會，那麼我們的采樣頻率至少就應該是40000HZ
實際上大部分的MP3，WAV，OGG等音樂媒體檔案，使用的采樣率是44100HZ，如果一個樣本點用16bits也就是2位元組來計算的話，每分鐘大約就需要5M位元組的資料量。
在多媒體的音樂檔案中，我們一般需要以下幾個參數，來确定多媒體檔案的聲音是如何采樣并如何播放的：
- SampleRate 采樣速率
- Channel 聲道數量
- BitPerSample 每個樣本的位數
依靠這幾個參數，我們就可以播放音樂檔案了
PCM碼流與WAV檔案格式
PCM 脈沖編碼調制是Pulse Code Modulation的縮寫，簡單來說就是抽樣、量化和編碼，也就是上章節分别對應的采樣速率、每個樣本的位數（量化）的概括了，它相當于信号的RAW（原始資料格式）。
WAV檔案（波形檔案）應該是最接近這種原始格式的檔案了，它沒有對資料額外的壓縮，由一個檔案頭構成後，剩下的都是原始的PCM資料，對播放參數進行配置後将它寫入聲霸卡就可以直接播放出聲音來了。
WAV的檔案頭如下（C++ 引用PainterEngine Audio代碼）
struct WaveHeader
{
pt_uchar riff[4]; //data exchange flag
pt_dword size; //filesize-header
pt_uchar wave_flag[4]; //wave
pt_uchar fmt[4]; //fmt
pt_dword fmt_len; //
pt_word tag; //format
pt_word channels; //
pt_dword samp_freq; //sample rate
pt_dword byte_rate; //samplerate*bytes of per sample
pt_word block_align; //channles * bit_samp / 8
pt_word bit_samp; //bits per sample
};
其中，最關鍵的幾個參數是channels，samp_freq和bit_samp。讀取資訊頭之後，往後搜尋data字元串，之後的資料都是PCM碼流了。
WAV檔案與頻譜分析
WavFreq是由C++編寫，UI架構使用Qt 4.8.6，播放接口使用DirectSound的一款簡易的聲波頻譜分析器。
每隔0.1秒，它會對聲波取樣4096個點并做FFT後将頻譜圖顯示。
作為本文的附錄程式，WavFreq的源代碼你可以在附件中找到，同時源碼遵循GPLv3開源協定。
軟體界面如圖g.1
圖g.1
然後打開一個WAV波形檔案，如圖g.2
圖g.2
打開後将顯示這個波形檔案的一些基本資訊，如圖g.3
圖g.3
點選菜單上的Start Analyse，觀察波形頻譜，如圖g.4
圖g.4
琴鍵分析與聲學模組化
為了更好的講解頻域在聲學分析上舉足輕重的作用，筆者以一段鋼琴按鍵音作為示範，範例檔案你同樣可以在附錄的檔案當中找到。
鋼琴中分别标注了七個不同的音，使用1,2,3,4,5,6,7,8（dou rui mi fa so la xi do）進行标注，它們的頻譜分别如下。
從上面的頻譜圖中我們可以看出，鋼琴音的頻率有着明顯的差别，随着音調的升高，主要的頻率也随之右移動。
在這裡，我們使用一些非常簡單的檢波濾波手段來區分這個時候按下的音是哪個。
首先我們先過濾掉那些并不主要的頻域信号，例如在上圖中，我們主要頻域在150~350的區間内，同時，我們取2500000的路徑成本作為其門檻值 以避免一些噪聲的幹擾。
這樣，鋼琴按鍵的識别過程就簡單變為了：
目前頻幅度量大于2500000時，在150~350HZ區間内找到幅值能量最高的點，依據該點所在頻率，濾出一些關鍵的特征頻率，并依此來判斷按鍵的類别。
具體的代碼與實作你仍然可以在附件中找到，運作結果如下圖所示（圖h.1）
圖h.1
被洩漏的密碼
如果你有看過那些間諜大片也許你對下面這個鏡頭場景并不陌生,例如《諜影重重》中就有這樣一個情節，特工想要進入目标的一個保險庫，于是他錄制了目标人物的聲音然後僞造了他的語音密碼，成功入侵了目标的保險庫。
另一個比較經典的鏡頭是語音的識别，特工錄制了目标的語音，很快就有裝置将語音識别成了文字和數字，假如目标在談話中或者是電話的按鍵音中洩露了密碼，那麼他就要倒大黴了。
不管情節如何或者這個到底具不具備可行性，上面兩個場景都和聲波的頻域分析有關，那麼從聲音還原出密碼具備多少的可行性呢。
在很多的直播錄制節目中，很多的主播在衆多的觀衆面前輸入自己的賬号與密碼，當然因為密碼是*字遮蓋的，觀衆無法直接看到密碼，但是不少的主播的鍵盤敲起來是非常的響的（例如青軸的機械鍵盤），我們可以很清楚的聽到敲擊鍵盤噼噼啪啪的聲音，有沒有可能通過對鍵盤音的分析模組化來還原出密碼呢。
筆者認為這絕對是有可能的，首先，Enter鍵 上檔鍵 空格鍵之類的鍵因為其模具的不同，其發出的敲擊聲音的頻率肯定是不同的，另外鍵盤經過長時間的敲打，因為磨損與雜質的不同其擊鍵音也有可能改變，最後一點是一個人長期的使用鍵盤，其對各個鍵的擊鍵的力度也是不一樣的，并且錄音裝置與鍵盤各鍵的距離不一樣，很可能從聲音的能量再做進一步的篩選。
當然，也許完全識别出密碼也許有困難，但是對主播擊鍵音長期的采樣分析建立模型，極大可能的縮小篩選範圍是絕對可行的。
當然，更好的聲音采樣器（采樣精度，采樣頻率）對分析的成功率也有更多的幫助，以此看來劣質的麥克風反而還保護了主播們的密碼隐私了。
筆者還未做過相關的實驗，讀者有興趣的話，不妨試試.
FFT2
二維沖擊采樣
還記得之前提到一維的采樣函數麼，我們對原函數進行了采樣操作，那麼拓展到二維方向，其沖擊采樣應該寫成了這種格式
如果對離散的函數進行采樣，那麼公式就由積分變為了累加
如果對特定樣本點進行采樣例如
隻需要對二維沖擊函數進行移位就行了
二維傅裡葉變換對
那麼我們很容易也将傅裡葉變換推到到二維的方面上來
其中，
是一個大小為M*N的矩陣，當然，逆變換同樣有
是以，二維的傅裡葉變換過程的算法就能簡單的概括為：
首先依次對一個二維矩陣的每一行做傅裡葉變換，得出變換的結果後，再對行變換結果的每一列做傅裡葉變換。
當然，為了友善計算機處理，其正變換僞代碼如下
1.設二維數組包含待變換的資料
2.對二維數組的每一行做傅裡葉變換，并将變換結果替換原數組。
3.将二維數組的行與列元素互換。
4.再對該二維數組的每一行做傅裡葉變換，并将變換結果放回原數組。
- 再将二維數組的行與列元素互換，其結果即為二維傅裡葉變換結果。
用C語言實作二維傅裡葉變換對
二維傅裡葉變換的C語言代碼如下，當中引用了之前編寫的一維傅裡葉變換函數，當然有了之前的基礎，二維傅裡葉變換的代碼簡單的多，您可以在本文的附錄中找到其完整的代碼：
void FFT_2(_IN complex x[],_OUT complex X[],int N)
{
for (int i=0;i<N;i++)
{
FFT(&x[i*N],&X[i*N],N);
}
//Matrix transpose
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
for(int i=0;i<N;i++)
{
FFT(&X[i*N],&X[i*N],N);
}
//Matrix transpose again
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
}
void IFFT_2(_IN complex X[],_OUT complex x[],int N)
{
//Matrix transpose
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
for(int i=0;i<N;i++)
{
IFFT(&x[i*N],&X[i*N],N);
}
//Matrix transpose again
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
for (int i=0;i<N;i++)
{
IFFT(&X[i*N],&X[i*N],N);
}
}
數字圖像的相位與頻譜
頻譜能量與相位
相信讀者從高中實體中經常看的到的能量公式，常常帶有平方關系，例如，能量是品質與速度的平方關系，能量是電壓的平方與電阻的關系。
在圖像中，我們常常也使用平方關系來表示“能量”這一關系，但這并不是說這張圖确實帶有實體上的多少能量，它更像一種比喻，就像古時群朝大臣跪地大喊吾皇萬歲萬歲萬萬歲一樣，或者避諱稱皇帝叫“萬歲”，但至少到今天，也沒有哪個皇帝有那能耐活得到一萬歲的，這頂多隻是口頭上的一個溜屁拍馬，或者說叫起來友善。
圖像上的能量也是類似這種關系，它更像是某種量化，并不是指它具體帶有的能量，不然得話，如果你想炸掉某樣東西，你隻要往他們那寄一張白紙就可以了（畢竟白色所帶的能量最高），更恐怖的是，要是你敢撕掉你的作業本，它甚至會引起一場爆炸。
回到正題，在圖像的頻譜中，我們應該如何表示這種能量的關系呢，顯然的，傅裡葉變換的結果是一個複信号，如果求其能量，僅需要對其實部與虛部求平方相加，當然，在頻譜圖中我們更多的是取複數的模。也就是對其能量進行一次開平方。
當然，僅僅有頻譜是無法表示一個複信号的，是以我們還需要引入相位的概念
其值是虛部除以實部的arctan值，當然，這也就表明它的範圍是
。
最後仍然值得提到的一個細節是，二維傅裡葉變換的結果是呈現中心共轭對稱的。（你仍然可以套用之前一維傅裡葉變換的證明方法來證明二維的共轭對稱性），證明過程在本文就不再複述了，但由共轭對稱性我們就很容易推到出頻域圖的中心對稱性（直流分量部分除外部分）。
圖像頻域與相位有什麼實體意義麼？
實際上《數字圖像處理》中有章節專門讨論這個問題，總結為圖像頻域表現的是灰階資訊，而相位則是位置資訊，頻域的低頻表示圖像基本的灰階變化，例如一個圖像灰階變化平緩也就意味着其頻率較低，而頻域的高頻則展現其灰階劇烈變化的部分例如物體的邊緣，星空的繁星（或者叫噪聲？）。
因為本文讨論的并不是如何對圖像進行哪些模糊銳化…之類的操作，是以更多的資訊我們這裡不再進行更多的讨論，我們僅僅需要知道，圖像的能量主要分布在了低頻當中，對頻域不過分的操作并不會特别影響圖像的視覺效果就可以了。
FFT2 DEMO程式示範
有了上述的理論基礎，對圖像進行二維離散傅裡葉變換的過程也就水到渠成了，筆者編寫了DEMO程式，你可以在本文的附錄中找到這個DEMO程式和它的代碼，在該程式中，首先打開的界面如圖所示（i.1）
圖i.1
選擇File-Open來打開一個圖檔檔案，如圖i.2
圖i.2
在這裡，就以筆者參照動漫《龍與虎》逢坂大河所繪制的同人圖為例（非專業畫師，不喜勿噴），因為是彩色圖檔，是以，我們分别對其的RGB（Red Green blue）分量進行分離，分别對它們做二維傅裡葉變換如圖i.3
圖i.3
其中，左邊第一張圖是原圖，第二列是原圖的RBG分量，第三列是其各分量的頻譜圖，第四列是其相位譜。
FFTShift
當然，在我們觀察到的頻譜圖中，我們更希望将頻譜顯示的更易于觀察，在上述的頻譜中，其低頻分量位于左上角，且因為傅裡葉變換的共轭對稱性，我們同樣很容易推導出二維傅裡葉變換是中心對稱的，假設我們将低頻分量移動到原點上向外為高頻，那麼圖像就會更加的易于觀察。
FFTShift的算法并不複雜，它實際上僅僅是對矩陣進行分割後（分割為4部分），對對角線的兩部分進行兩兩交換。如圖i.4
圖i.4
例如如下矩陣
經過FFTShift後就變為了
經過FFTShift後，DEMO的頻譜與相位圖亦發生相應的改變
圖i.5
可以看到，原本分散于四個角的低頻信号經過移位到中心後，變得更加的易于觀察了。
對數變換
有句話叫真理往往掌握在少數人的手中，在圖像的二維傅裡葉變換的頻譜能量圖中，從上個章節我們很容易觀察到，頻譜的能量往往在中心才較為的明顯（低頻域），顯然的，圖像的能量往往集中在低頻當中。因為低頻與高頻的能量差距過大，為了便于頻域圖像的觀察，我們對頻域進行對數操作，并對操作的結果縮放到0-255的灰階級别中，假設
F(x,y)
表示圖像頻域複信号的模值（就是能量開根号），那麼對數變換就是
下圖是經過了對數變換的頻域圖像，可以看到，其頻域灰階變化變得更加的平緩了（圖i.6）
圖i.6
魯棒盲水印
盲水印與版權
圖像水印被廣泛運用在了防僞，簽名，辨別等版權保護方面，簡單來說就是防止盜圖狗竊取自己的勞動成果或者将他抓個現行。畢竟網絡大了什麼樣的人都有，被人竊取成果反而署名上其他人的名字是一件極其令人憤慨的事情。
目前的水印主要是明水印（可見水印）居多，例如下圖（j.1）使用的就是明水印
J.1
明水印加起來簡單友善并且快捷，使用Photoshop來處理應該是一件非常簡單的事情，不需要太多的技術含量就可以給這張圖“冠名”了
當然，簡單的事情往往攻擊起來也很簡單，稍微懂點的盜圖者很容易就能把明水印删除并還原原圖，實在嫌麻煩可以直接裁剪圖檔，很多時候裁去水印并不會對圖像整體的視覺上造成多大影響（如圖j.2）。
圖j.2
當然，最煩人的是，水印常常破壞對原圖的視覺享受，這在影像原畫作品中，常常是不能被觀衆容忍的，一個作品突然冒冒失失的加上一個水印，常常導緻觀衆看起來就像心理有個疙瘩。
那麼有沒有水印技術，直接看不見，經過處理後就變得可見了呢。
這實際就是本文另一個要讨論的技術細節，實際上這種盲水印技術在國小的動手實驗上就有，最知名的應該是用米湯寫字，寫完後紙上是看不見的，如果要顯示資訊，那麼就用點碘液一洗，字就變藍色顯示出來了。
使用流程圖來表示這一過程（如圖j.3 j.4）
j.3
j.4
在數字圖像中，目前非常流行的一種盲水印隐寫技術，是對圖像像素操作的，其具體流程是，因為一個顔色具有RGB分量，一般來說，RGB三個顔色通道每個各占8位，也就是三位元組，對8位的最低位做修改，對原圖的影響是微乎其微的，那麼，每個像素就能帶有3位的資訊，一張100*100的24位位圖，就能夠帶有30000位3750位元組的資訊可插入，如果将水印資料插入到這裡，就可以實作一個簡單的盲水印了。
這種盲水印技術具有一定的可行性，但是，它的缺點仍然也是非常緻命的，對原圖像的旋轉，平移，縮放，改變色相都能非常輕易地摧毀水印。
于是，具備抗攻擊的盲水印技術，也就孕育而生了。
基于傅裡葉變換的魯棒頻域水印
從之前的章節我們已經了解了如何利用二維傅裡葉變換将圖像變換為頻域，并且我們也了解了圖像的能量主要集中在低頻當中，那麼如果我們将水印疊加在頻域的高頻中，理論上對原圖的視覺上并不會有過多的影響，同時，在頻域中的水印散列分部在空域（也就是逆變換後的圖像）的各個部分，對頻域水印的破壞将會變得更加的困難，同時，頻域水印對圖像的裁切，旋轉，平移，加噪都有一定的抗攻擊能力（詳細的推到可翻閱《數字圖像處理》一書），是以，頻域水印作為一種盲水印手段是擁有其理論依據的，對于這種具有一定抗攻擊能力的水印技術，我們又稱之為魯棒水印。
最後，作為讨論的細節，如何将水印疊加到頻域也是應該讨論的範疇，在這裡，長話短說地總結以下幾點：
- 對彩色圖像加水印，首先對圖像的RGB顔色通道分别分離，對R，G，B三個通道的顔色分别計算頻域，就和前幾個章節處理的那樣。
- 疊加混合的方式分為兩種，稱之為縮放疊加，一種為放大一種為縮小，因為二維傅裡葉變換後的結果是一個複信号，是以，如果我們僅僅修改頻譜能量而不影響其相位的話，應該将複信号的實部與虛部同比例放大或縮小就可以了。是以，水印也就是安裝水印輪廓對對應複信号進行放大或縮小
- 峰值信噪比（PSNR），歸一化相關系數（NC值）可用于判定水印對原圖的影響及水印的抗幹擾能力，兩者都是越高越好，當然，這常常是一個互相沖突的度量，往往抗幹擾能力強也意味着對原圖的幹擾大。
水印程式DEMO
ImageSigner是一款由筆者開發的圖像水印程式，作為一個示範的DEMO程式，它僅僅支援256*256大小的圖檔，你可以在本文的附錄中找到它的完整源代碼，您可以通過查閱源代碼，來了解頻域水印的具體算法及過程，它使用C++與Qt Framework 4.8.6編寫完成，遵循GPLv3開源協定。下面介紹其水印簽名及各種攻擊效果。
1.打開程式，界面如圖（k.1）
圖k.1
2.點選Open Source Image和Open Sign Image分别加載源圖與簽名圖檔（圖k.2 k.3）。
圖k.2
K.3
載入後界面如圖所示（k.4）
圖k.4
在簽名控制台，設定簽名的參數（圖k.5）
圖k.5
其中，第一欄表示簽名的通道，第二欄為簽名水印的混合模式，一般來說，enlarge對簽名圖像的抗幹擾能力更強，Reduce模式抗幹擾模式較弱，但Enlarge模式對原圖的影響較大。
Power一欄表示對簽名的混合能量（值越大表示水印越明顯），一般來說，能量越大對原圖的幹擾也較大。
選擇完成後，點選Do Sign進行簽名（圖k.6）
圖k.6
在右側你可以看到簽名後的圖檔，點選Save to file可以将簽名後的檔案儲存。 最後，我們選擇簽名的圖像，并分别使用裁剪，平移，旋轉，噪聲，塗抹，改變色相，縮放攻擊并檢視攻擊後水印的保留效果。
下面圖k.7為原圖，k.8 k.9分别為縮小及放大混合簽名後的圖像，可以看到視覺差別不大。
k.7
k.8
K.9
同時，簽名後圖像的頻譜如k.10 K.11（放大水印隻加在藍色通道）所示
K.10
K.11（為友善觀察依情況引用）
使用PhotoShop對圖像進行攻擊
裁剪攻擊與其頻譜
平移攻擊及其頻譜
旋轉攻擊及其頻譜
高斯噪聲與頻譜
塗抹攻擊與頻譜
改變色相攻擊與其頻譜
縮放攻擊（裁剪後放大至原尺寸）與頻譜
後記
文章到了這裡也許你對當中的技術意猶未盡，也許你看的一臉雲裡霧裡不知所雲，但不管如何，《從三角函數到語音識别再到數字水印》到了這裡也應該告一段落了。
我一直有個願望，把一個看上去複雜的東西，從它最原始的最簡單的根基開始，抽絲剝繭，循序漸進，一步一步地深入直到它開花結果，這個過程，就猶如見證了一場破繭化蝶。因為我相信，那些深入宇宙奧妙與規律的知識，是足以讓人震撼的，它并不輸于任何不可估價的藝術品。
遺憾的是，如今國内不論是學術還是技術都有一股浮躁之風，有人将其當做宣揚炫耀的籌碼，有人甚至剽竊他人的成果與勞動，作為自己謀利的肮髒台階，在我看來，他們都不是真正的“為學”者，我希望有更多的人，能夠沉下耐心，将前輩的知識打碎發揚，将那些看起來遙不可及的智慧，分解成一步一階的台階，讓更多的人更快更簡單地得到智慧的福澤。而不是光顧着告訴大家我掌握了這個技術有多麼“了不起”。
我并不知道讀者們閱讀本文時，是否能收到我所傳達的意思與願望，本文行文時間将近三個月，當中每一個文字每一行标點每一句代碼包括每一張例示圖檔都出自本人之手都凝聚着我的心血，但我并非聖賢，文章并不能做到句句嚴謹，行行精辟，相信有不少的缺漏錯誤還待讀者們斧正，但知識理當能夠有所發揚，如果您在我的字裡行間最終有所感悟并了解了“她”是怎麼一回事，那麼，您的一句“寫的不錯”将是我最大的欣慰與鼓舞。
最後，我還想再說些什麼，但現在我卻有些詞窮了，我将ImageSigner做更進一步的改善，讓它能夠支援到2次幂高寬的圖像數字簽名。當然，你仍然可以在附件中找到它的完整源代碼與Release版本并将它運用在實際的需求中。如果您有什麼相關的疑問或讨論，可以發送郵件到[email protected]
行文至此，那麼，就這樣吧。-------------------------------DBinary 2017年6月6日
github代碼下載下傳位址
https://github.com/matrixcascade/ImageSignergithub.com