信號調製方式識別在無線電管理、電子對抗等應用中佔據了十分重要的地位，從1969年Waver C S等人發表第一篇調製方式自動識別的論文以來，在該領域不斷有人提出新方法，例如Liang Hong、K.C.Ho採用小波變換識別FSK、PSK、16QAM三種數字調製信號[1]；Gardner將循環譜分析用於信號調製識別[2]；Assaleh等人把信號建模為一個兩階AR過程，並利用參數統計方法識別CW、PSK、FSK三類信號[3]。
信號調製識別一般包括兩個重要的部分，即類間（Inter－class）識別和類內（Intra－class）識別，本文著重研究FSK、PSK和多音FDM三類信號的類間識別問題，由於多音FDM是多載波信號，需採用時頻方法進行分析，但單純使用時頻不能很好地反應信號的特徵，為此，本文首先介紹了將傅氏變換應用於時頻分布各頻帶的聯合頻率分析方法，並通過DSP信號闡述了該方法的特性，然後根據譜相關循環頻率軸譜表徵信號循環平穩特性的優點，將其取代傅氏變換得到聯合信號頻率與循環頻率的兩維處理方法，並用於三種數位訊號的特徵分析。最後詳細介紹了基於聯合頻率分析的特徵提取及識別過程，給出了仿真測試結果。
1 聯合頻率分析
1.1 基於傅氏變換的聯合頻率分析
設信號為x（t）,首先對其瞬時自相關函數做傅氏變換，得到關於時間和頻率的兩維函數，即著名的Wigner－Ville時頻分布[4]：
其中μ為調製頻率，下面以雙邊帶幅度調製信號（DSB）為例，說明聯合頻率分析的特性，DSB信號的表達式為：
從（4）式可以看出，在信號頻率與調製頻率聯合平面上，存在多個非0值，且這些值出現的位置具有對稱性（本文只分析μ≥0，f≥0的情況）。顯然，在|μ|≤2fm範圍內的非0值，反應了調製信號的一些特徵。相比信號載頻，這些特徵一般集中於較低的調製頻率（μ）處，其他範圍內也有非0值，主要因為Wigner－Ville分布（SWD）、平滑偽Wigner－Ville分布（SPWD）[4]等。由於採用FFT計算傅氏變換會造成聯合分析平面包含大量高頻冗餘信息，降低分析效率，且通信信號一般具有循環平穩特徵，因此本文設計採用譜相關μ截面分析取代第二代傅氏變換，下面介紹其基礎知識。
1.2 譜相關理論
設信號x（t）微循環平穩且功率有限，則其在時間區間[－T，T]上的循環自相關函數為：
信號的譜相關是一種形式上的兩維傅立葉頻譜，兩個變量分別是F1和a，令f1＝0，得到循環頻率a軸上的譜相關Sxa（0），簡稱為a軸譜，a軸譜一般包含了與信號載頻、符號速率等有關的重要信息[5]。
1.3 基於軸譜的數字調製信號聯合頻率分析
按（1）式計算出信號的Wigner－Ville時頻分布後，分別固定每個頻率f，沿時間軸方向計算其a軸譜，得到變量分別為信號頻率f和循環頻率a的兩維頻譜Sx（a，f）。在本文中，聯合頻率分析將用於移頻鍵控（FSK）、移相鍵控（PSK）、多音頻分復用（M－tone FDM）三類數字調製信號的類間識別，下面結合時頻分布來分析採用軸譜後信號的聯合頻率特性。
M進位FSK信號採用M個不同的頻率來傳輸信息，圖1（a）為一個碼速為80bps的BFSK信號時頻分布圖，從圖中可看出，不同的信號頻率分布在不同的頻帶內，每一時刻只出現一個頻率。同時，在其聯合頻率分布圖1（b）中相應的頻帶（這裡指與a軸平行的頻率子帶，下同）內，出現一些譜峰。在a＝0時出現的譜峰（俯視圖中的黑點，下同），對應的信號頻率表徵了BFSK信號的頻率參數；a＝80Hz時也出現譜峰，其反應信號的碼速信息，顯然，MFSK（M＞2）信號存在M個頻帶具有此類似的特徵。
M進位PSK信號通過對載波的相位調製得到，圖2（a）一個載頻和碼速分別為1400Hz和200bps的BPSK信號時頻分布圖，信號相位的變化使得在載頻及其附近頻帶內，周期性出現一些頻率分量。顯然在載波頻帶內，這種周期性表現尤為強烈。如聯合頻率分布圖2（b）所示，a＝200Hz時，在載波頻帶內出現的譜峰表徵了碼速信息，MPSK（M＞2）信號的特徵與此類似。
多音FDM信號一般由多個相互獨立或正交的子信號和一個單載波導頻疊加而成。圖3（a）為一個12音FDM信號，每個子信號的調製方式為BPSK，符號速率為100bp，從其時頻圖中可看出，導頻在其頻帶內分布是均勻的，每路BPSK與前面分析的情況類似，聯合頻率分布圖如圖3（b）所示，由於導頻不包含調製信息，其在相應頻帶內分布也是均勻的，而各路BPSK在a＝0和a＝100Hz時出現譜峰。
從上述分析中可以得出，三類信號的聯合頻率分布有明顯的區別，如MFSK在M個頻帶內存在譜峰，MPSK信號僅在載波頻帶內存在譜峰，多音FDM信號雖在多個頻帶內存在譜峰，但導頻頻帶與子信號頻帶內譜峰分布不同。因此，聯合頻率分布很好地表徵了信號的特徵，可用於調製方式識別及信號的參數估計。
2 特徵提取及仿真實現
2.1 特徵提取
待分析信號的採樣率為16kHz，持續時間為1s，採用平滑偽Wigner－Ville分布計算信號時頻分布，平滑窗長度256點，窗移20定，對得到的時頻分布矩陣，利用幅度平方處理分別檢測每個頻帶的包絡，並使其通過低通濾波器，以濾除一些毛刺，同時，為了避免各頻帶大直流分量的幹擾，利用時頻矩陣的平均值，對整個矩陣進行歸一化，然後設定循環頻率a的分析區間[0，ax]，分別對每個頻帶計算器a軸譜，最終得到以信號頻率m和循環頻率n為變量的聯合頻率矩陣Sx（m，n）。
同其他的二維分析一樣，聯合頻率分析Sx（m，n）所需的處理時間也較長，與傳統的一維譜估計方法所提供的數據量相比一般要大得多，即使利用信號帶寬範圍來選擇Sx（m，n）中的分析區間，得到的特徵矩陣仍然太大而無法提供分類器使用，由於矩陣的奇異值是矩陣的固有特徵，比例不變性，因此，將奇異值分解（SVD）方法應用於聯合頻率矩陣特徵的提取，得到特徵矢量：C＝[σ1，…，σp]，這裡，p＝min（m，n），σ1，i＝1，…，p為奇異值，同時，Sx（0，n）中包含的譜峰數也將作為一個重要特徵用於調製識別。
2.2 識別性能測試
本文在Matlab平臺上對三類信號的類間識別進行了仿真實驗，在0－20db（步進為5dB）的信噪比範圍內，按隨機消息序列分別產生MFSK（M＝2，4，8）、MPSK（M＝2，4、8）及M音頻分復用（M＝8，12，16）三類信號，每一類信號在每個信噪比下的樣本數均為1000，然後按2.1節中的方法隨機選擇500個分別提取特徵組成訓練集，剩餘500個樣本的特徵組成測試集，分類器採用徑向基（RBF）神經網絡，在對分類器訓練之前，先根據信噪比的不同，將各組訓練集兩兩交叉組合，分別得到[0dB，5dB]、[5dB，10db]、[10dB，20dB]、[5dB、15db]和[0db、20dB]等多個訓練集組合，經不同組合訓練出的分類器，對三類信號的平均識別成功率有所不同，其他[5dB，15dB]組合對應的分類器性能最好，相應的識別結果如圖4。
本文提出的聯合信號頻率和循環頻率分析方法，將時頻分析與譜相關理論有機地結合起來，是描述信號特徵的一種有效工具，已成功應用於FSK、PSK和多音FDM三類信號的類間識別，由於聯合頻率平面包含了信號載頻、碼速等重要參數信息，也可用於後續的參數估計環節。