?2012年9月的時候，一個南京的大學生從電視臺播放的一段記者采訪360總裁周鴻祎的視頻中破解了周鴻祎的手機號碼，一時間被網絡熱炒。后來，又聽說某人買車的時候使用電話銀行付款，結果被人錄下聲音，破解了銀行卡號和密碼，導致存款被盜。最近居委會在小區里散發傳單，提醒一種新的詐騙方式：電話聽音破密碼詐騙。各種網絡和媒體對這些事情炒作的很熱鬧，但是趨于兩種極端，要么將其說的出神入化，要么將其貶的一錢不值，唯獨缺少理性的聲音解釋一下到底是什么原理。貓貓一直想說，“其實就是離散傅立葉變換（DFT）嘛，至于嗎？”，這次就借著這種所謂的新型詐騙方式，一并說道說道。

??????? 在動筆寫這篇文章之前，貓貓上網Google了一下“電話號碼破解”，居然發現某寶上有人叫賣“電話撥號音破解器”，價格從幾百到上千大洋的都有，貓貓不說了，只笑笑。本文就準備從理論上解釋一些破解電話撥號音的數學原理，看了本文的同學，你有福了，只需要依照本文給出的百十行演示代碼，稍加改造，再加上一個麥克風，你就可以擁有自己的“電話撥號音破解器”了，廢話少說，現在就開始吧。

??????? 首先要說的是，根據電話撥號音破解電話號碼，只適用于使用“雙音多頻技術（DTMF）”的電話設備，老式的撥號盤電話（脈沖式電話機）您就別試了（估計您也沒有這玩意了）。雙音多頻技術是貝爾實驗室的發明，就是將電話機的撥號鍵盤分成4X4的矩陣，每一行對應一個低頻信號，每一列對應一個高頻信號，如圖（1）所示：

圖（1）電話鍵盤雙音頻對照表

?

其中低頻信號和高頻信號的頻率都在人耳可以識別的頻率范圍之內。打電話撥號的時候，每按下一個鍵，就產生一個高頻信號和低頻信號的正弦信號組合，局端的電話交換機從這個組合信號中解出兩個頻率，就知道是那個按鍵被按下了。說到這里，您可能已經猜到了，只要從按鍵音中解析出這兩個頻率，就能根據上表查出是對應的哪個按鍵了。

??????? 現在該離散傅立葉變換（DFT）隆重登場了！關于離散傅立葉變換算法的來龍去脈，請上網Google之，這里只簡單介紹一下它的基本原理，以及它在電話號碼破解的過程中能做什么。

??????? 傅里葉是一位法國數學家和物理學家的名字，他于1807年在法國科學學會上發表了一篇論文，提到了一個觀點：任何連續周期信號可以由一組適當的正弦曲線組合而成，換言之，滿足一定條件的連續函數（周期函數）都可以表示成一系列三角函數（正弦和/或余弦函數）或者它們的積分的線性組合形式，這個轉換就被稱為“傅立葉轉換”。一般情況下，若“傅里葉變換”一詞的前面未加任何限定語，則指的是“連續傅里葉變換”，與之對應的自然就是“離散傅立葉變換”。事實上，變換原始信號可以是周期的，也可以不是周期的，根據原始信號類型的不同，傅立葉變換可以分成四種形式，分別是“非周期性連續信號傅立葉變換”、“周期性連續信號傅立葉變換”、“非周期性離散信號傅立葉變換”和“周期性離散信號傅立葉變換”，其中“周期性連續信號傅立葉變換”實際上就是傅里葉級數的推廣，因為積分其實是一種極限形式的求和算子而已。

??????? 一維傅立葉變換的數學意義就是將時域信號轉換為頻域信號，當然，也存在稱為傅立葉逆變換的變換，可以將頻域信號轉換成時域信號。連續時域信號以及對應的連續傅里葉變換都是連續函數（且是無限長度的時域信號），但是用計算機進行數字處理只能處理有限長度的離散信號，所以必須將原始信號離散化，同時建立對應的離散信號傅里葉變換。所謂的“離散傅里葉變換（Discrete Fourier Transform，縮寫為DFT）”，是指傅里葉變換在時域和頻域上都以離散的形式呈現。原始信號離散化的過程其實就是以一定的周期對原始信號進行采樣的過程，至于如何從連續傅立葉變換推導出離散傅立葉變換，不是本文的重點，您可以Google之。因為離散傅立葉變換計算量相當大，有很多提高效率的算法理論，其中應用最廣泛的就是快速傅立葉變換（FFT），關于快速傅立葉變換的理論依據，也不是本文的重點，您可以Google之。總之，您只需要知道，存在一種方法，可以將基于時間軸的時域信號，轉換成基于頻率的頻域信號即可，本文的重點是這種轉換的物理意義，以及如何應用這種轉換帶給我們的結果。

圖（2）440Hz正弦波在時域和頻域的形態

?

??????? 現在看幾個轉換對比圖，理解一下這種轉換的物理意義。首先看圖（2），圖（2-a）是440Hz的正弦波在時域中的形態，采樣率是8000Hz。圖（2-b）就是通過離散傅立葉轉換后得到的頻域中的形態，理想狀態下，圖（2-b）應該在440Hz處一條直線，其他位置的值都是0，但是受原始信號雜波和轉換后的頻域分辨率影響，實際顯示的是一個呈金字塔形狀的圖形，不過還是可以明顯地看到在440Hz的時候功率（相對強度）值最大，其他位置的值都明顯小于440Hz位置的值。

圖（3）按鍵“1”對應的音頻在時域和頻域的形態

?

??????? 圖（3-a）是電話按鍵“1”對應的音頻波形，圖（3-b）是離散傅立葉轉換后得到的頻域形態，雖然受錄音雜波影響很大，但是還是可以很明顯看到在697Hz和1209Hz位置上，相對能量強度達到了最大值。

??????? 前面提到過，原始信號離散化的過程其實就是以一定的周期對原始信號進行采樣的過程，這里就提到了一個很重要的參數，就是采樣率T。因為這個采樣率T決定了轉換成頻域信號后的分辨率，另外，每次進行轉換的時域信號的個數N（也稱為離散傅立葉變換的點數）不能小于T，否則無法完整映射到頻域信號。實際上，是T和N共同決定了轉換后的頻域信號的分辨率。以本文要用到的聲音信號為例，轉換前的數據是每個采樣點上的聲波的振幅（各個頻率疊加后的振幅），橫坐標軸是T周期的時間點，縱坐標軸是振幅。轉換成頻域信號后，橫坐標軸就是頻率，縱坐標軸就是該頻率的產生強度，通過適當的轉換，可以理解為該頻率的功率。剛才提到，時域信號的采樣率T和參與轉換的信號點數N決定了頻域信號的分辨率，這是因為離散傅立葉轉換會將時域信號一對一的轉換為頻域信號，也就是說，N個時域信號轉換后會得到N個頻域信號。但是傅立葉轉換有個特點，就是時域信號的周期性和頻域信號的對稱性，所謂頻域信號的對稱性，是指轉換后的N個頻域信號在T/2位置呈現軸對稱特征，對信號分析有用的實際頻域信號點只有N/2 +1個。其中第一個點對應的是0Hz，也就是時域信號中直流分量強度，其他N/2個點分別對應(0,T/2]之間的頻率，頻率分辨率是T/N。由此可見，增加N可以提高頻域信號的分辨率，但是增加N會顯著增加每次傅立葉變換的計算量，因此對N的選擇需要做一個折衷。

??????? 使用離散傅立葉變換進行信號分析的時候，一般都會采用快速傅立葉變換（FFT）算法，網上有很多高效的快速傅立葉變換算法，包括基于實數的快速傅立葉算法，大家可以參考。本文給出的算法是貓貓上學時寫的一個快速傅立葉變換算法，輸入和輸出參數都是復數，算不上高效，但是很簡單，只有二十幾行代碼，適合研究原理使用。另外，這個算法采用原位變換的方式，可以減少傅立葉轉換算法使用過程中的數據復制操作。

?42?/* ?43?window = 1 -> hanning window ?44?*/ ?45?bool?InitFft(FFT_HANDLE?*hfft,?int?count,?int?window) ?46?{ ?47?????int?i; ?48? ?49???? hfft->count?=?count; ?50???? hfft->win?=?new?float[count]; ?51?????if(hfft->win?==?NULL) ?52?????{ ?53???? ????return?false; ?54?????} ?55???? hfft->wt?=?new?COMPLEX[count]; ?56?????if(hfft->wt?==?NULL) ?57?????{ ?58???? ????delete[]?hfft->win; ?59???? ????return?false; ?60?????} ?61?????for(i?=?0;?i?<?count;?i++) ?62?????{ ?63???? ??? hfft->win[i]?=?float(0.50?-?0.50?*?cos(2?*?M_PI?*?i?/?(count?-?1))); ?64?????} ?65?????for(i?=?0;?i?<?count;?i++) ?66?????{ ?67???? ????float?angle?=?-i?*?M_PI?*?2?/?count; ?68???? ??? hfft->wt[i].re?=?cos(angle); ?69???? ??? hfft->wt[i].im?=?sin(angle); ?70?????} ?71? ?72?????return?true; ?73?} ?74? ?75?void?FFT(FFT_HANDLE?*hfft,?COMPLEX?*TD2FD) ?76?{ ?77?????int?i,j,k,butterfly,p; ?78? ?79?????int?power?=?NumberOfBits(hfft->count); ?80? ?81?????/*蝶形運算*/ ?82?????for(k?=?0;?k?<?power;?k++) ?83?????{ ?84???? ????for(j?=?0;?j?<?1<<k;?j++) ?85???? ????{ ?86???? ??? ??? butterfly?=?1<<(power-k); ?87???? ??? ????for(i?=?0;?i?<?butterfly/2;?i++) ?88???? ??? ????{ ?89???? ??? ??? ??? p=j?*?butterfly; ?90???? ??? ??? ??? COMPLEX t?=?TD2FD[i?+?p]?+?TD2FD[i?+?p?+?butterfly/2]; ?91???? ??? ??? ??? TD2FD[i?+?p?+?butterfly/2]?=?(TD2FD[i?+?p]?-?TD2FD[i?+?p?+butterfly/2])?*?hfft->wt[i*(1<<k)]; ?92???? ??? ??? ??? TD2FD[i?+?p]?=?t; ?93???? ??? ????} ?94???? ????} ?95?????} ?96? ?97?????/*重新排序*/ ?98?????for(i?=?1,j?=?hfft->count/2;i?<=?hfft->count-2;?i++) ?99?????{ 100???? ????if(i?<?j) 101???? ????{ 102???? ??? ??? COMPLEX t?=?TD2FD[j]; 103???? ??? ??? TD2FD[j]?=?TD2FD[i]; 104???? ??? ??? TD2FD[i]?=?t; 105???? ????} 106???? ??? k?=?hfft->count/2; 107???? ????while(k?<=?j) 108???? ????{ 109???? ??? ??? j?=?j?-?k; 110???? ??? ??? k?=?k?/?2; 111???? ????} 112???? ??? j?=?j?+?k; 113?????} 114?}

?

??????? InitFft()函數初始化一個傅立葉轉換句柄，包括正/余弦系數表（wt）和窗口系數表（win），count參數是傅立葉轉換的點數，就是一次能轉換的最大信號個數，windows參數是決定使用什么窗函數，這段代碼只支持漢寧窗（hanning window）。有關離散傅立葉變換中窗函數的作用將在后面介紹。FFT()函數就是快速傅立葉變換的實現，TD2FD是個復數數組，每調用一次FFT()函數可以轉換count個數據。對于音頻數據這樣的實數，只需使用COMPLEX結構的實數分量，虛數分量可以置0。COMPLEX數據結構定義如下：

?8?struct?COMPLEX ?9?{ 10?????float?re; 11?????float?im; 12?}; 13?

COMPLEX數據結構重載了+、-和*運算符，所以FFT()函數的代碼看起來很簡練。

??????? 前文提到過，計算機不能處理無限長度的數據，離散傅立葉變換算法只能對數據一批一批地進行變換，每次只能對限時間長度的信號片段進行分析。具體的做法就是從信號中截取一段時間的片段（比如本文對按鍵音頻文件的分析，從音頻文件中讀取的數據就可以看做是從無限長度信號中截取的一秒鐘音頻信號），然后對這個片段的信號數據進行周期延拓處理，得到虛擬的無限長度的信號，再對這個虛擬的無限長度信號進行傅立葉變換。但是信號被按照時間片截取成片段后，其頻譜就會發生畸變，這種情況也被稱之為頻譜能量泄露。

??????? 為了減少能量泄露，人們研究了很多截斷函數對信號進行截取操作，這些截斷函數就被稱為窗函數。窗函數w(t)被設計成頻帶無限的函數，所以即使原始信號是有限帶寬信號，被窗函數截取后得到的片段也會變成無限帶寬，也就是說，信號經過窗函數處理后，在頻域的能量與分布都被擴展了，有效地減少了頻譜能量泄露。

??????? 不同的窗函數對信號頻譜的影響是不一樣的，比如最簡單的矩形窗，實際上就是對信號不做任何處理，簡單地按照時間片段截取一定長度的信號進行處理。本文的演示算法選用了漢寧窗（hanning window），漢寧窗的數學定義如下：

其中R(t)是原始信號，t的范圍是 0 <= t < N-1，對于其他范圍的值，w(t) = 0。圖（4）顯示了漢寧窗對信號截取的示意圖，以及對頻域轉換結果的影響，藍色區域是窗口覆蓋的數據部分。漢寧窗的作用是分析帶寬加寬，但是降低了頻率分辨率，不過撥號音的低頻和高頻組中，任意兩個頻率之間的間隔都超過了70Hz，所以犧牲這點分辨率對按鍵聲音的識別沒有影響。

圖（4）漢寧窗信號截取示意圖

?

??????? 使用了窗口，就需要討論窗口的滑動問題，也就是窗口重疊的處理。用漢寧窗截取的信號片段，可以看出來窗口中部分信號被削弱了（造成衰減），為了抵消部分窗口對信號造成的衰減，各種窗函數都需要對信號進行相應的重疊處理。本文選擇的重疊處理方式是選取信號時每次滑動半個窗口位置，使得每個窗口的后半個窗口的衰減在下個窗口的前半個窗口中的到一定的補償，這個在計算一段時間信號總的功率譜時會有相應的體現。計算一段時域信號的總功率譜的函數代碼如下：

31?bool?PowerSpectrumT(FFT_HANDLE?*hfft,?short?*sampleData,?int?totalSamples,?intchannels,?float?*power) 32?{ 33?????int?i,j; 34? 35?????for(i?=?0;?i?<?hfft->count;?i++) 36???? ??? power[i]?=?(float)0.0; 37? 38???? COMPLEX?*inData?=?new?COMPLEX[hfft->count]; 39?????if(inData?==?NULL) 40???? ????return?false; 41? 42?????int?procSamples?=?0; 43?????short?*procData?=?sampleData; 44?????while((totalSamples?-?procSamples)?>=?hfft->count) 45?????{ 46???? ??? procData?=?sampleData?+?procSamples?*?channels; 47???? ????for(j?=?0;?j?<?hfft->count;?j++) 48???? ????{ 49???? ??? ??? SampleDataToComplex(procData,?channels,?&inData[j]); 50???? ??? ??? procData?+=?channels; 51???? ????} 52???? ??? procSamples?+=?(hfft->count?/?2);?/*每次向后移動半個窗口*/ 53? 54???? ??? FftWindowFunction(hfft,?inData); 55???? ??? FFT(hfft,?inData); 56? 57???? ????for(i?=?0;?i?<?hfft->count;?i++) 58???? ????{ 59???? ??? ??? power[i]?+=?float(20.0?*?log10(sqrt(inData[i].re?*?inData[i].re?+inData[i].im?*?inData[i].im)?/?(hfft->count?/?2))); 60???? ????} 61?????} 62? 63?????delete[]?inData; 64? 65?????return?true; 66?}

??????? PowerSpectrumT()函數要做的事情很簡單，就是每次從采樣數據中選擇count個，對其進行窗口處理，然后用傅立葉轉換得到頻域信號，最后計算并累加各個頻率的相對功率。處理完一批采樣數據后，向后偏移 count/2 個采樣數據（相當于滑動半個窗口），重復上述過程，直到剩余的采樣數據全部完成為止（剩余數據不足一個窗口時也退出計算）。sampleData參數是采樣數據，本文處理的是16位音頻數據，因此每個采樣的大小是兩個字節，totalSamples是總的采樣數據的個數，totalSamples必須大于count才能得到有意義的功率頻譜。channels是采樣數據的聲道數，對于有兩個聲道的立體聲采樣數據，則取兩個聲道的平均值做為采樣數據進行傅立葉轉換。SampleDataToComplex()函數將采樣數據賦值給復數對象的實數部分，虛數部分置0，對于雙聲道，則取左右聲道的平均值賦值給復數對象的實數部分，其實現如下：

17?static?void?SampleDataToComplex(short?*sampleData,?int?channels,?COMPLEX?*cd) 18?{ 19?????if(channels?==?1) 20?????{ 21???? ??? cd->re?=?float(*sampleData?/?32768.0); 22???? ??? cd->im?=?0.0; 23?????} 24?????else 25?????{ 26???? ??? cd->re?=?float(*sampleData?+?*(sampleData?+?1)?/?65536.0); 27???? ??? cd->im?=?0.0; 28?????} 29?}

??????? 現在來講一下對PowerSpectrumT()函數轉換結果的處理，PowerSpectrumT()函數的power參數返回最終計算好的功率譜，其分辨率是T/N。每個點對應的頻率可使用以下公式計算：

?

Freq = n * (T / N)??????????????????????????????????????????????? （公式2）

?

以音頻信號采樣率8000Hz，傅立葉轉換的count是2048為例，其頻譜分辨率是3.90625Hz，power[0]是頻率0的相對功率，也就是直流信號的功率，power[1]是頻率3.9Hz對應的相對功率，power[2]是頻率7.8Hz對應的相對功率，對應關系以此類推。如果在power[240]處出現頻率最大值，則其對應的頻率應該是（8000/2048）*240=937Hz。

??????? 雙音頻電話按鍵音中有兩個主頻率，因此轉換后的power數組中應該存在兩個相對極大值，需要從power數組中找出這兩個相對極大值。從M個數中選最大的N個數是一個很常見的算法，其主要思想就是維護一個大小是N的有序組，然后從剩下的M-N個數中逐個與有序組中的最小值比較，如果小于有序組中的最小值，則繼續比較下一個數，如果大于有序組中的最小值，則替換有序組中的最小值，并對有序組進行適當的變換，使其保持有序。使用堆來維護有序組被認為是最高效的方式，大家有興趣可以從網上找相關的代碼學習。本文設計的從power數組中選出最大的兩個值的算法也采用這種思想，但是有序表只有兩個元素，因此對有序表的維護不需要復雜的排序處理，如果有序性失效，只要交換兩個元素的位置即可，其具體算法如下：

290?void?ExchangeIndex(int?*index,?float?*power) 291?{ 292?????if(power[index[1]]?>?power[index[0]]) 293?????{ 294???? ????int?t?=?index[0]; 295???? ??? index[0]?=?index[1]; 296???? ??? index[1]?=?t; 297?????} 298?} 299? 300?void?SearchMax2FreqIndex(float?*power,?int?count,?int&?first,?int&?second) 301?{ 302?????int?max2Idx[2]?=?{?0?}; 303? 304???? ExchangeIndex(max2Idx,?power); 305?????for(int?i?=?2;?i?<?count;?i++) 306?????{ 307???? ????if(power[i]?>?power[max2Idx[1]]) 308???? ????{ 309???? ??? ??? max2Idx[1]?=?i; 310???? ??? ??? ExchangeIndex(max2Idx,?power); 311???? ????} 312?????} 313? 314???? first?=?min(max2Idx[0],?max2Idx[1]); 315???? second?=?max(max2Idx[0],?max2Idx[1]); 316?}

??????? SearchMax2FreqIndex()函數從power數組中選擇最大的兩個值，并返回它們在power數組中的位置（power數組的下標），并根據頻率排序，first總是較小的那個頻率的位置，second總是較大的那個頻率的位置。根據first和second指示的位置，利用“公式2”就可以計算出兩個頻率，但是因為頻譜分辨率的關系，這兩個計算出來的頻率不會剛好就是圖（1）中所示的標準頻率，但是沒關系，它們與標準頻率差的不多，通過查表就可以確定是哪個標準頻率。

??????? 至此，所有的核心代碼都完成了，創建一個例子試試吧，我的例子很簡單，如圖（5）所示，就是打開按鍵音頻文件，然后顯示是哪個按鍵。您可以做的更復雜一些，直接從麥克風獲得音頻數據，實時分析按鍵聲音，這樣就可以實時破解電話號碼了。

圖（5）按鍵音識別實例程序

總結

以上是生活随笔為你收集整理的听声音破解电话号码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。