數(shù)字信號相似度方法研究及GPU并行加速

朱靜靜，謝曉堯*，劉志杰，王 培

(1．貴州師范大學(xué) 信息與計算科學(xué)重點實驗室，貴州 貴陽 550001；2.中國科學(xué)院國家天文臺，北京 100101)

0 引言

近年來，隨著社會經(jīng)濟(jì)和科技的迅猛發(fā)展，我們正處于大數(shù)據(jù)時代，而這些大數(shù)據(jù)都是具有時間屬性的時間序列，相似性度量在其中得到了廣泛的應(yīng)用。時間序列是由多項與時間順序有關(guān)的數(shù)據(jù)記錄組成的元素的有序集合[1]，它存在于金融收益、網(wǎng)絡(luò)安全、氣象研究、天文研究等[2]，相似性度量是指將兩個數(shù)據(jù)對象的相似度進(jìn)行比較，是衡量數(shù)據(jù)對象間相似程度的指標(biāo)，常常應(yīng)用于數(shù)據(jù)分析和數(shù)據(jù)挖掘過程中，對于模式識別、信息檢索、聚類、分類、序列模式、神經(jīng)網(wǎng)絡(luò)等研究中起著非常重要的作用，具有很高的研究價值[3]。對數(shù)據(jù)進(jìn)行深入研究，從而更好地理解數(shù)據(jù)信息間的關(guān)系以及更方便有效地提取出潛在的隱藏價值。

計算相似度的方法較多，例如各種距離測度法、相對誤差距離法、余弦相似度度量法、相關(guān)系數(shù)法、最大相異系數(shù)法、特征值相似比法[4-5]。經(jīng)典的時間序列相似性度量方法總體分為兩種：鎖步度量(lock-step measures)和彈性度量(elastic measures)[6]。鎖步度量是兩個時間序列進(jìn)行點對點或一對一的比較[7]，最常見的是歐氏距離[8]，Keogh等人和FALOUTSOSC等人分別提出利用歐氏距離[9][10]進(jìn)行相似性度量的方法。彈性度量是兩個時間序列進(jìn)行一對多或一對零的比較，常見的有動態(tài)時間規(guī)整和基于編輯距離的度量[11]。

信號相似度評估結(jié)果可為噪聲信號來源判定、仿真信號置信度評價等應(yīng)用提供一種判斷依據(jù)[12]。數(shù)字信號處理中，計算信號相似度的方法有很多，如按定義計算；通過快速傅里葉變換；互相關(guān)的庫函數(shù)等?；ハ嚓P(guān)是對兩個序列相似性的一種度量，是一個信號相對于另一個信號位移的函數(shù)，稱為滑動點積?；ハ嚓P(guān)常用于在長信號中尋找較短信號的已知功能，這種相似性度量方式在模式識別、電子斷層掃描、密碼分析等領(lǐng)域中都有應(yīng)用，互相關(guān)在含有周期成分的信號中提取特定頻率成分信息；圖像匹配；線性定位和相關(guān)測速以及計算信號時延等都有著至關(guān)重要的應(yīng)用意義。

近年來，提高數(shù)據(jù)處理效率一直是計算機(jī)研究領(lǐng)域的熱點，為了有效提高大規(guī)模計算的處理速度與性能，中央處理器(Central Processing Unit，CPU)+圖形處理器(Graphics Processing Unit，GPU)的異構(gòu)并行計算架構(gòu)得到了進(jìn)一步應(yīng)用與發(fā)展。

本文以500 m口徑球面射電望遠(yuǎn)鏡 (Five-hundred-meter Aperture Spherical Telescope， FAST) 19波束接收機(jī)的數(shù)據(jù)進(jìn)行實驗，實現(xiàn)了CPU上的4種計算信號相似度的方法，并進(jìn)行對比和分析，針對大規(guī)模數(shù)據(jù)采用互相關(guān)函數(shù)定義的方式計算相似度耗時較長，提出了2種利用CPU與GPU異構(gòu)的并行架構(gòu)方法提高大規(guī)模數(shù)據(jù)計算效率。FAST望遠(yuǎn)鏡19波束接收機(jī)數(shù)據(jù)計算各波束間或者某些信號之間的相似性，可以檢測、識別甚至提取一些特定成分的信息，同時，在信號線性定位以及計算信號時延等方面都有重要的研究價值。

1 相似度

時間序列是將研究數(shù)據(jù)的變化過程按時間先后順序記錄下來而形成的序列[13]，時間序列相似性的定義由Agrawal等人在1993年提出[14]，數(shù)字信號的相似性度量是計算時間序列間的相似性。

在信號處理中經(jīng)常要研究兩個信號的相似性，或一個信號經(jīng)過一段延遲后自身的相似性，相關(guān)函數(shù)廣泛應(yīng)用于信號中隱含周期性的檢測, 確定未知參數(shù)的線性系統(tǒng)的頻域響應(yīng), 噪聲中信號中的檢測、提取等[15]。信號x(n)和y(n)的互相關(guān)函數(shù)的定義為:

(1)

該式表示Rxy(m)在時刻m時的值，等于將x(n)保持不動而y(n)左移m個抽樣周期后兩序列對應(yīng)相乘再相加的結(jié)果[16]，它是描述隨機(jī)信號x(n)和y(n)在任意兩個不同時刻t1，t2的取值之間的相似程度。

當(dāng)y(n)=x(n)時，上面定義的互相關(guān)函數(shù)變成自相關(guān)函數(shù)，即：

(2)

Rxx(m)反映了信號x(n)和其自身經(jīng)過一段延遲后的x(n+m)的相似程度[16]。自相關(guān)函數(shù)表達(dá)了同一過程不同時刻的相互依賴關(guān)系，而互相關(guān)函數(shù)表示不同過程的某一時刻的相互依賴關(guān)系。

信號x(n)和y(n)的線性卷積的定義為：

(3)

為了和(1)式的互相關(guān)函數(shù)比較，將(3)式中的m和n互換，則

(4)

而信號x(n)和y(n)互相關(guān)函數(shù)：

(5)

故相關(guān)與卷積的時域關(guān)系為：Rxy(m)=x(-m)*y(n)。

(6)

相關(guān)函數(shù)和線性卷積在計算形式上非常相似，它們的區(qū)別在于計算x(n)和y(n)的互相關(guān)時，兩個序列均不需要翻轉(zhuǎn)，只是將y(n)在時間軸上移動后對應(yīng)相乘再相加；而計算卷積時，需要先將一個序列翻轉(zhuǎn)后再移動。二者所表示的物理意義是截然不同的，線性卷積表示LSI系統(tǒng)的輸入、輸出和單位抽樣響應(yīng)之間的一個基本關(guān)系，而相關(guān)反映的只是兩個信號之間的相似程度，與系統(tǒng)無關(guān)[16]。

2 GPU、Pycuda和scikit-cuda

在高級游戲中，每個視頻幀需要執(zhí)行大量的浮點運(yùn)算，為了滿足這一需求且減小經(jīng)濟(jì)壓力，圖形處理器(Graphics Processing Unit，GPU)的設(shè)計理念油然而生[17]。GPU最初是一種圖形加速器，隨著科技的發(fā)展，GPU逐漸演化成一個強(qiáng)大的、多用途的、完全可編程的，以及任務(wù)和數(shù)據(jù)并行的處理器，廣泛應(yīng)用于大規(guī)模的并行計算；CPU適合處理控制密集型任務(wù)，GPU適合處理包含數(shù)據(jù)并行的計算密集型任務(wù)[18]，將GPU與CPU結(jié)合有效地提高大規(guī)模并行計算問題的處理速度與性能。

CUDA(Compute Unified Device Architecture)是由NVIDIA在2007年推出的一種編程模型，旨在支持GPU/CPU應(yīng)用程序的執(zhí)行[17]，它利用NVIDIA GPU中的并行計算引擎能更有效地解決復(fù)雜的并行計算問題。CUDA平臺可以通過CUDA加速庫、編譯器指令、應(yīng)用編程接口以及行業(yè)標(biāo)準(zhǔn)程序語言的擴(kuò)展(包括C、C++、Fortran、Python)來使用[18]。

PyCUDA是通過Python訪問NVIDIA的CUDA并行計算的API，PyCUDA可以充分利用CUDA驅(qū)動程序API的全部功能，具體是把CUDA C 程序嵌套在Python腳本中訪問CUDA并行計算的API，通過提高代碼抽象性來降低編程復(fù)雜性和改善編碼效率。

PyCUDA對NumPy的數(shù)組提供了很好的支持,相對純C/C++的CUDA程序而言,PyCUDA簡化了初始化工作,可以直接將NumPy的數(shù)組作為參數(shù)傳遞給GPU核函數(shù)使用,并且直接用該數(shù)組存儲GPU計算之后的結(jié)果[19]。

scikit-cuda為CUDA設(shè)備或運(yùn)行時、CUBLAS、CUFFT和CUSOLVER庫中的許多函數(shù)提供了Python接口，這些庫是作為NVIDIA的CUDA編程工具包的一部分發(fā)布的。

3 計算信號相似性的方法

計算信號相似性的方法很多，本文主要討論與分析了按互相關(guān)定義計算、傅里葉變換、Python中的Numpy和Scipy庫以及CPU+GPU異構(gòu)并行架構(gòu)的兩種方法計算信號相似度，得到一系列相關(guān)系數(shù)或相關(guān)函數(shù)值。Numpy是一種開源的數(shù)值計算擴(kuò)展，可用來存儲和處理大型矩陣。Scipy是基于Numpy的一個用于數(shù)學(xué)、科學(xué)、工程領(lǐng)域的常用軟件包，可以處理插值、積分、優(yōu)化、圖像處理、常微分方程數(shù)值解的求解、信號處理等問題。

3.1 互相關(guān)定義與快速傅里葉變換計算相似度

1)互相關(guān)函數(shù)定義方法

2)快速傅里葉變換方法

快速傅里葉變換(Fast Fourier Transformation , FFT)和逆變換(Inverse Fast Fourier Transformation,IFFT)計算互相關(guān)函數(shù)。對于N點有限長信號，通過FFT將信號的時域信息轉(zhuǎn)換為頻域信息，把X(n)信號頻域的共軛乘以Y(n)信號的頻域得到兩個信號的互相關(guān)函數(shù)的頻域信息，然后把頻域恢復(fù)到時域，此時得到的數(shù)據(jù)信息范圍是[0 ～N-1, -(N-1) ～ 0]，而直接時域計算得到的數(shù)據(jù)范圍是[-(N-1) ～N-1]，所以需要利用fftshift函數(shù)對fft處理后得到的結(jié)果的順序進(jìn)行調(diào)整。

3.2 Python中的Numpy與Scipy庫計算相似度

Python中的numpy.correlate()和scipy.signal.correlate()兩個函數(shù)計算互相關(guān),一維的數(shù)據(jù)可以直接計算二者的相似度，但是實際處理的數(shù)據(jù)往往都是多維數(shù)據(jù)，需要對這些數(shù)據(jù)降維，對于二維數(shù)據(jù)，為了計算出完整的相關(guān)系數(shù)序列，需要用到一些“不存在”的點，這就需要對這些數(shù)據(jù)進(jìn)行補(bǔ)零操作，即對超出原始數(shù)據(jù)存儲區(qū)的部分取值為零。前端和末端的數(shù)據(jù)并沒有完全“嵌入”到原始數(shù)據(jù)的全部信息，會受到補(bǔ)零操作的影響，這些數(shù)據(jù)往往是不可用的。

3.3 CPU+GPU異構(gòu)并行計算相似度

按定義計算信號相關(guān)的方法計算復(fù)雜度約O(N2)，當(dāng)數(shù)據(jù)量大時，耗時較多，由于沒有進(jìn)行補(bǔ)零等操作，定義法計算的結(jié)果精度較高，采用CPU+GPU的異構(gòu)并行計算方式提高數(shù)據(jù)處理效率。

1)利用PyCUDA計算

Python是一種廣泛使用的解釋型、高級編程、通用型的編程腳本語言，PyCUDA是把CUDA C 程序嵌套在Python程序中，利用GPU+CPU即把串行處理數(shù)據(jù)放在CPU上運(yùn)行，把可以并行處理的數(shù)據(jù)放在GPU上運(yùn)行，提高按定義計算信號相似度的處理效率。共享內(nèi)存(Shared memory)是片上內(nèi)存，帶寬僅次于寄存器的存儲器，速度快于一般的全局內(nèi)存(Global memory)和常量內(nèi)存(Constant memory)等顯存，實驗中使用共享內(nèi)存減少線程塊中的線程必須訪問的數(shù)據(jù)總量，根據(jù)數(shù)據(jù)規(guī)模設(shè)置相應(yīng)的最優(yōu)的網(wǎng)格和塊大小計算相似度，進(jìn)一步提高數(shù)據(jù)處理速度。

2)利用scikit-cuda

scikit-cuda為CUDA設(shè)備或運(yùn)行時、CUBLAS、CUFFT和CUSOLVER庫中的許多函數(shù)提供了Python接口，這些庫是作為NVIDIA的CUDA編程工具包的一部分發(fā)布的，還提供了在CULA密集型工具包中選擇函數(shù)的接口，提供了與C類似的低級包裝函數(shù)和與NumPy和Scipy類似的高級函數(shù)。實驗中主要用到了skcuda.linala.dot()函數(shù)可以自適應(yīng)的計算兩個信號的外積，然后進(jìn)一步計算相關(guān)系數(shù)。

4 實驗與結(jié)果分析

4．1 實驗

本文實驗平臺為GeForce GTX 1080 Ti GPU，Intel(R) Core (TM) i7-6700 CPU@ 3.40 GHz，16 GB 內(nèi)存和Red Hat Enterprise Linux Server 7.1操作系統(tǒng)。使用 Python和 PyCUDA(Python+ CUDA C)編寫程序，集成開發(fā)環(huán)境為PyCharm，同時環(huán)境支持為CUDA Toolkit 9.2。

以FAST望遠(yuǎn)鏡19波束接收機(jī)的FITS格式數(shù)據(jù)進(jìn)行實驗，每個文件是單個波束在一次觀測中所記錄的數(shù)據(jù)，由于每個波束的FITS文件約2 G，19波束約38 G數(shù)據(jù)，且原始數(shù)據(jù)“SUBINT”數(shù)據(jù)中的DATA部分?jǐn)?shù)據(jù)是五維數(shù)組，由于計算速度和計算機(jī)內(nèi)存的需求和限制，需要對數(shù)據(jù)進(jìn)行降維、兩路偏振合并成一路偏振以及采樣等預(yù)處理操作，得到一個第一維表示時間，第二維表示頻率通道的19波束的二維數(shù)組數(shù)據(jù)。實驗中使用的是二維數(shù)組的某一頻率通道的19波束的部分時間序列數(shù)據(jù)，分別為5 120×1、10 240×1、15 360×1、20 480×1、25 600×1、30 720×1六個計算規(guī)模大小，并對按定義計算信號相似度的方法進(jìn)行加速。

4.2 結(jié)果分析

CPU上的4種方法包括FFT、Python中的2個庫函數(shù)以及互相關(guān)定義法，通過實驗發(fā)現(xiàn)Python中的2個庫函數(shù)Numpy.correlate()和Scipy.signal.correlate()兩種方法計算的結(jié)果基本一樣，如圖1所示，上面兩部分是FAST望遠(yuǎn)鏡19波束某頻率通道的部分時間序列經(jīng)過預(yù)處理后的二維數(shù)字信號數(shù)據(jù)，第三部分是這兩個數(shù)字信號的相似度的可視化表示，可以發(fā)現(xiàn)這兩個時間序列(數(shù)字信號)相似度很低，且兩種方法計算的結(jié)果基本相同。

圖1 Numpy和Scipy庫評估相似度結(jié)果Fig.1 The evaluate similarity results of Numpy and Scipy libraries

通過FFT方法和Python中的2個庫函數(shù)方法計算FAST望遠(yuǎn)鏡19波束接收機(jī)的一個時間段的某頻率通道的兩個信號的相似度，經(jīng)過預(yù)處理后的數(shù)據(jù)是二維數(shù)組，為了得到完整的計算結(jié)果，對于N個數(shù)據(jù)點，在計算過程中都需要對其補(bǔ)充N-1個零。對于FFT的方法是先分別對兩個信號進(jìn)行傅里葉變換，然后第一個信號傅里葉變換的復(fù)共軛乘以第二個信號的傅里葉變換值，最后再用fftshift函數(shù)調(diào)整數(shù)據(jù)的順序。

CPU上的4種方法計算兩個信號的相似度的結(jié)果基本一致，如圖2所示，按照定義的方式計算的結(jié)果只有時滯大于零的部分，而其他3種方法計算的相似度結(jié)果是有時滯大于零和小于零兩部分，因為計算過程中只有按照定義的方式?jīng)]有對數(shù)據(jù)進(jìn)行補(bǔ)零操作，而其他3種方法都進(jìn)行了補(bǔ)零操作，這些補(bǔ)零操作從某種程度上影響了相似度的精確度。

圖2 CPU上的4種方法評估相似度結(jié)果Fig.2 The evaluate similarity results of four ways on the CPU

在CPU上的4種方法對于不同數(shù)據(jù)規(guī)模所需時間的差異，經(jīng)過多次實驗，選取了一些相對理想的數(shù)據(jù)，不同方法所需時間如圖3所示，很明顯的看出計算速度最快的是Scipy.signal.correlate()方法，按照定義方法計算相對耗時，Numpy.correlate()和FFT方法計算效率介于二者之間。雖然按照定義方法計算相似度耗時，但是由于沒有進(jìn)行補(bǔ)零等操作，計算的精確度相對較高。

圖3 CPU上4種方法計算不同矩陣規(guī)模的時間Fig.3 The time of calculate different matrix scales of four methods on the CPU

由于按照相關(guān)函數(shù)定義的方法計算信號相似度相對耗時，所以采用CPU+GPU異構(gòu)并行的2種方式提高計算效率，計算結(jié)果如圖4所示，這兩種方法是對按定義計算信號相似度方法的改進(jìn)，計算的結(jié)果只有時滯大于零的數(shù)據(jù)，且這兩個時間序列相似度較高。PyCUDA方法需要根據(jù)數(shù)據(jù)規(guī)模調(diào)整網(wǎng)格和數(shù)據(jù)塊的大小，不能自適應(yīng)的計算，而scikit-cuda庫中的culinalg.dot()函數(shù)自適應(yīng)計算，降低編程復(fù)雜度，數(shù)據(jù)處理效率更高。

圖4 PyCUDA和scikit-cuda庫評估相似度結(jié)果Fig.4 The evaluated similarity results of PyCUDA and scikit-cuda libraries

多次實驗顯示，利用PyCUDA計算矩陣的外積。對于10 240×1規(guī)模，最高提速11.8倍，對于20 480×1規(guī)模矩陣，最高提速12倍，對于30 720×1規(guī)模矩陣，比只用CPU計算耗時，后面計算具體相關(guān)函數(shù)值通過PCI-E總線傳輸數(shù)據(jù)到CPU上進(jìn)行計算，數(shù)據(jù)規(guī)模不斷增大，顯存和內(nèi)存間大量的數(shù)據(jù)通信耗時較多，帶來了較大的時間損耗，這在一定程度上削弱了GPU 快速進(jìn)行高性能并行計算的能力，整體提速不明顯，需要進(jìn)一步研究。scikit-cuda庫函數(shù)方式計算信號相似度提速效果較好，如圖6所示，30 720×1信號矩陣最高提速7.3倍，并行處理的性能提升并不是線性增長，隨著數(shù)據(jù)規(guī)模的增加，受到PCI-E總線帶寬的影響和GPU核心數(shù)目的限制，加速比會逐漸趨于平緩，實驗中只使用了FAST望遠(yuǎn)鏡19波束某一頻率通道的部分時間序列數(shù)據(jù)，數(shù)據(jù)量逐漸增大時加速比會趨于平緩。

圖5 CPU+GPU的scikit-cuda庫與CPU定義方法加速比Fig.5 The acceleration ratio of CPU+GPU scikit-cuda library and CPU definition method

CPU與GPU的計算信號相似度方法結(jié)果對比如圖6所示，可以看出多種方法計算兩個數(shù)字信號相似度的結(jié)果基本一致，按照定義的方式計算的結(jié)果較精確，通過CPU+GPU的異構(gòu)并行的方式顯著提高數(shù)據(jù)的處理效率。

圖6 CPU與GPU的計算相似度方法結(jié)果對比Fig.6 The comparison of results of the CPU and GPU calculate the similarity methods

5 總結(jié)

本文討論并分析了CPU上的4種計算信號相似度的方法以及結(jié)合GPU的異構(gòu)并行方式計算信號相似度。實驗結(jié)果表明：CPU的4種方法中Scipy.signal.correlate()方法計算效率最高，按定義計算的精度較高，當(dāng)數(shù)據(jù)規(guī)模大時，按定義法計算相對耗時，所以采用CPU+GPU的異構(gòu)并行的2種方式提高數(shù)據(jù)的處理效率。實驗表明：基于定義的方式計算信號相似度，通過PyCUDA實現(xiàn)，整體提速效果不明顯，需要進(jìn)一步深入研究，Scikit-cuda庫的使用顯著提高了定義方式計算相似度的數(shù)據(jù)處理效率，在數(shù)據(jù)處理中具有良好的應(yīng)用前景。

計算FAST望遠(yuǎn)鏡19波束中各波束間數(shù)據(jù)的相似性，或者某些信號之間的相似性，可以檢測、識別甚至提取一些特定成分的信息，特殊信號線性定位以及計算信號時延等方面值得進(jìn)一步深入研究。