魯曉帆，劉志剛，吳 峰

（西南交通大學(xué) 電氣工程學(xué)院，四川 成都 610031）

0 引言

基于小波變換的小波包變換，不僅在時(shí)頻域都具有良好的局部化性質(zhì)，而且對信號的高頻部分提供了更精細(xì)的分解，從而得到大量高頻信息[1-3]。小波包變換在電力系統(tǒng)多方面的應(yīng)用都有明顯的優(yōu)勢，例如在數(shù)據(jù)壓縮領(lǐng)域，小波包變換可以針對電力系統(tǒng)故障錄波信號中蘊(yùn)含豐富暫態(tài)特征信號的特點(diǎn)，將故障信號分解成各個(gè)頻帶的信息，選擇有重要信息的頻帶采用各種壓縮算法進(jìn)行壓縮。近些年，在小波包算法基礎(chǔ)之上，又提出了與嵌入式零數(shù)編碼相結(jié)合[4]、與快速傅里葉變換相結(jié)合[5]、混合小波包[6]等各種提高壓縮精度的壓縮算法。這些算法雖然在一定程度上使壓縮性能得到提升，但卻增加了算法復(fù)雜度。隨著電力系統(tǒng)數(shù)據(jù)記錄裝置的飛速發(fā)展、電力系統(tǒng)高采樣的逐步演變、濾波裝置的不斷更新，電力系統(tǒng)電能質(zhì)量數(shù)據(jù)逐漸呈現(xiàn)海量化趨勢，傳輸、存儲等環(huán)節(jié)對實(shí)時(shí)性的要求又進(jìn)一步提高，因此，電力系統(tǒng)數(shù)據(jù)壓縮效率問題得到越來越多的關(guān)注。以上各種壓縮算法所基于的小波包變換，其核心部分的卷積運(yùn)算過程復(fù)雜、運(yùn)算量大、實(shí)時(shí)性差，嚴(yán)重阻礙了其在實(shí)際工程中的應(yīng)用。為了提高小波包變換的計(jì)算速度，解決目前電力系統(tǒng)海量數(shù)據(jù)的壓縮效率問題，本文提出將小波包變換并行化處理。

目前，針對這一問題，國內(nèi)外學(xué)者提出利用各種并行技術(shù)來加快小波包變換的運(yùn)算速度[7-13]，大部分文獻(xiàn)都是利用多機(jī)網(wǎng)絡(luò)環(huán)境進(jìn)行并行化處理，該方法不可避免地受到系統(tǒng)維護(hù)及網(wǎng)絡(luò)傳輸環(huán)節(jié)所帶來的不利影響。隨著計(jì)算機(jī)的發(fā)展，單機(jī)多核處理器應(yīng)用越來越廣泛，因此基于多核技術(shù)的并行小波包研究具有十分重要的意義。

本文利用 Pthreads[14]和 OpenMP[15-17]2 種并行編程環(huán)境，針對提高小波運(yùn)算速度問題，提出小波包Mallat分解重構(gòu)算法不同的并行策略。在Pthreads環(huán)境下，提出小波包分解過程數(shù)據(jù)級并行、重構(gòu)過程任務(wù)級并行；在OpenMP環(huán)境下，分別提出小波包分解重構(gòu)特征嵌套和非嵌套并行策略。并利用C語言在單機(jī)雙核處理器平臺上進(jìn)行了海量數(shù)據(jù)壓縮實(shí)驗(yàn)研究，取得了明顯的效果，為解決電力系統(tǒng)數(shù)據(jù)壓縮效率問題提供一種新的方法。

1 小波包Mallat算法并行性分析

小波包是小波概念的推廣。它是在小波變換（V0=V1⊕W1=V2⊕W2⊕W1=…）只將 V（尺度）空間進(jìn)行分解的基礎(chǔ)上，進(jìn)一步對Wj（小波空間）進(jìn)行分解，使正交小波變換中隨j的增大而變寬的頻譜窗口進(jìn)一步變細(xì)，達(dá)到最適合于待分析信號的時(shí)頻窗口或最優(yōu)基[18]。其分解過程如圖1所示，圖中h0與h1分別表示小波包分解后的高頻系數(shù)與低頻系數(shù)，a表示尺度空間數(shù)據(jù)，d表示小波空間數(shù)據(jù)。

圖1 小波包Mallat分解原理Fig.1 Principle of Mallat decomposition of wavelet packet

小波包分解重構(gòu)過程的核心是濾波器與待處理數(shù)據(jù)的卷積過程。從卷積計(jì)算的公式可以了解到，卷積計(jì)算結(jié)果與待處理數(shù)據(jù)序列直接相關(guān)，不同位置的待處理數(shù)據(jù)與濾波器的卷積過程得到不同的卷積結(jié)果。因此不同位置的卷積結(jié)果之間并沒有相關(guān)性，可以將待處理數(shù)據(jù)進(jìn)行分配并行卷積。

小波包Mallat重構(gòu)原理如圖2所示，圖中g(shù)表示由小波空間數(shù)據(jù)與尺度空間數(shù)據(jù)重構(gòu)后得到的上一層的重構(gòu)數(shù)據(jù)?？梢钥闯觯〔ò儞Q中對尺度空間和小波空間的處理也不存在關(guān)聯(lián)性，該部分同樣可以作為并行任務(wù)執(zhí)行。

圖2 小波包Mallat重構(gòu)原理Fig.2 Principle of Mallat reconstruction of wavelet packet

2 基于Pthreads與OpenMP的并行小波包實(shí)現(xiàn)

2.1 Pthreads與OpenMP

Pthreads是Linux操作系統(tǒng)的多線程接口標(biāo)準(zhǔn)，在Windows操作系統(tǒng)上，可以利用Pthreads-win32開放源代碼版本實(shí)現(xiàn)并行編程。與Pthreads庫不同，OpenMP是通過采用指導(dǎo)語句、庫函數(shù)調(diào)用和環(huán)境變量來給用戶提供在共享存儲計(jì)算機(jī)上創(chuàng)建和管理可移植并行程序的方法。它是一個(gè)在多處理機(jī)上用以編寫可移植的多線程應(yīng)用程序的API庫，采用Fork-Join并行模式。

2.2 基于Pthreads的并行小波包實(shí)現(xiàn)

小波包對信號數(shù)據(jù)處理的流程，需要經(jīng)過數(shù)據(jù)加載、延拓、分解、二抽取、二插值、重構(gòu)這一過程。利用Pthreads實(shí)現(xiàn)分解過程數(shù)據(jù)級并行分解、重構(gòu)過程任務(wù)級重構(gòu)。其中，延拓、二抽取與二插值的過程與分解、重構(gòu)過程融合為同一個(gè)流程。

a.數(shù)據(jù)延拓。為了執(zhí)行數(shù)據(jù)延拓，將每組與后一組重復(fù)NFilter-1個(gè)數(shù)據(jù)，其中NFilter是濾波器系數(shù)個(gè)數(shù)，本文采取周期延拓，最后一組與第1組數(shù)據(jù)重復(fù)前NFilter-1個(gè)數(shù)據(jù)，假設(shè)有p個(gè)線程參與執(zhí)行并行計(jì)算，則每組有N／p+NFilter-1個(gè)數(shù)，N為原始數(shù)據(jù)個(gè)數(shù)，卷積時(shí)濾波器第1位只對應(yīng)每組前N／p個(gè)數(shù)進(jìn)行計(jì)算。

b.數(shù)據(jù)級并行分解。由于小波包分解時(shí)各數(shù)據(jù)序列之間沒有相關(guān)性，可以將各個(gè)數(shù)據(jù)分配給各個(gè)線程同時(shí)執(zhí)行。為保證各個(gè)執(zhí)行線程負(fù)載平衡，本文根據(jù)小波包分解層數(shù)將近似系數(shù)與細(xì)節(jié)系數(shù)的數(shù)據(jù)分別分為N／p組，每層N的大小為上一層的1／2，數(shù)據(jù)分組如圖3所示。

圖3 數(shù)據(jù)分組Fig.3 Data grouping

將分解的前q組數(shù)據(jù)分給q個(gè)線程，每個(gè)線程同時(shí)執(zhí)行分解過程，待線程掛起后，主線程繼續(xù)向這些子線程派發(fā)新的q組數(shù)據(jù)。分解一層后，先對近似系數(shù)進(jìn)行分組并行分解，再對細(xì)節(jié)系數(shù)分組執(zhí)行。此時(shí)需要注意的是，由于延拓方式的影響，對應(yīng)于每一層中近似與細(xì)節(jié)系數(shù)的延拓值將會不同，因此在分組完成后先要計(jì)算此次延拓后的余值項(xiàng)，再采用新的余值進(jìn)行延拓與卷積計(jì)算。

c.二抽取。為了減少小波包變換計(jì)算時(shí)間，可以將分解中二抽取過程與分解過程相結(jié)合[20]，即不再計(jì)算所有數(shù)據(jù)，而是選擇奇數(shù)位數(shù)據(jù)計(jì)算，偶數(shù)位數(shù)據(jù)直接跳過。圖4采用濾波器系數(shù)為4說明此過程。圖中當(dāng)數(shù)據(jù)延拓分組后，不同的線程針對本組數(shù)據(jù)進(jìn)行卷積計(jì)算時(shí)，每個(gè)線程執(zhí)行第1、3、5、…位為首的運(yùn)算。

圖4 數(shù)據(jù)二抽取分解原理Fig.4 Principle of second extraction decomposition

圖5為小波包Pthreads并行分解流程，圖中Ti為不同的線程，di為不同組的數(shù)據(jù)，主線程完成數(shù)據(jù)載入與延拓后，將數(shù)據(jù)分組，分出一組數(shù)據(jù)后創(chuàng)建一個(gè)線程執(zhí)行該組的分解運(yùn)算，完成近似系數(shù)的分解，接著繼續(xù)分解第i層細(xì)節(jié)系數(shù)，直到所有系數(shù)分解完成。

由于數(shù)據(jù)回存需用到上一步線程的結(jié)果，并且是在同一堆棧中進(jìn)行讀寫，為了保證數(shù)據(jù)回存的正確性，各線程完成卷積計(jì)算以后，需在分解數(shù)據(jù)回存之前與主線程合并，主線程完成分組后在此等待，創(chuàng)建的線程得到釋放。

d.二插值與任務(wù)級重構(gòu)。重構(gòu)的過程中可以發(fā)現(xiàn)，每一層對應(yīng)的近似系數(shù)和細(xì)節(jié)系數(shù)個(gè)數(shù)是相同的，并且不同的系數(shù)與不同的濾波器進(jìn)行運(yùn)算，不同類型系數(shù)之間沒有相關(guān)性，讀取數(shù)據(jù)實(shí)現(xiàn)卷積的過程就不會出現(xiàn)同步的問題，因此可以將近似系數(shù)和細(xì)節(jié)系數(shù)重構(gòu)的過程看成是不同的任務(wù)，實(shí)現(xiàn)任務(wù)分解方式及任務(wù)大小相同，這也滿足負(fù)載平衡要求。

圖5 小波包Pthreads并行分解Fig.5 Parallel wavelet packet decomposition on Pthreads

小波包重構(gòu)并行過程，先根據(jù)每一層循環(huán)與分解層數(shù)的關(guān)系，計(jì)算數(shù)據(jù)延拓后的系數(shù)值，二插值時(shí)不進(jìn)行補(bǔ)零操作，而是在非補(bǔ)零位置計(jì)算待處理數(shù)據(jù)與濾波器奇數(shù)位的系數(shù)卷積，補(bǔ)零位置計(jì)算待處理數(shù)據(jù)與濾波器偶數(shù)位的系數(shù)卷積。而整個(gè)過程中尺度空間與小波空間同時(shí)進(jìn)行重構(gòu)運(yùn)算，其任務(wù)級并行重構(gòu)流程如圖6所示。

圖6 小波包Pthreads并行重構(gòu)Fig.6 Parallel wavelet packet reconstruction on Pthreads

2.3 基于OpenMP的并行小波包實(shí)現(xiàn)

OpenMP標(biāo)準(zhǔn)可以直接在串行程序的基礎(chǔ)上調(diào)整程序的并行性，因此利用OpenMP實(shí)現(xiàn)小波包并行必須首先編寫小波包分解的串行程序，本文在Visual studio 2008平臺上利用C語言完成串行編程。OpenMP最顯著的用處在于，可以將小波包for循環(huán)分配到不同處理器上同時(shí)執(zhí)行，即根據(jù)處理器的個(gè)數(shù)p將循環(huán)次數(shù)分為p組，每組循環(huán)由一個(gè)線程執(zhí)行，這樣由p個(gè)線程同時(shí)完成不同次循環(huán)體中卷積、二抽取計(jì)算，如圖7所示。

圖7 OpenMP并行卷積過程Fig.7 Process of parallel convolution on OpenMP

小波包分解串行程序核心部分是一個(gè)嵌套循環(huán)過程。嵌套循環(huán)的含義是把小波包分解層次與計(jì)算低頻、高頻系數(shù)的次數(shù)相聯(lián)系。為處理邊界問題，對邊界進(jìn)行周期延拓，共享數(shù)據(jù)作為私有副本分配給各個(gè)線程。根據(jù)圖2，假設(shè)分解n層（n≥1），外層循環(huán)需要進(jìn)行2n-1次。假設(shè)外層分解至第j（j=0，1，2，…，2n-2）次，內(nèi)層循環(huán)則從 jN／（2n-1）開始至（j+1）N／（2n-1）次結(jié)束。

a.嵌套方案。在OpenMP內(nèi)，嵌套是一個(gè)默認(rèn)不使能的布爾屬性。當(dāng)并行區(qū)域嵌套出現(xiàn)在線程已經(jīng)運(yùn)行在并行區(qū)域又遇到另一個(gè)并行區(qū)域時(shí)，若嵌套使能，那么程序就會按照之前動(dòng)態(tài)線程的規(guī)則生產(chǎn)一個(gè)新的線程組。可以通過omp_set_nested（i）來設(shè)置嵌套使能狀態(tài)，i=0表示嵌套未使能，i=1表示嵌套使能。

小波包分解過程利用嵌套使能并行執(zhí)行嵌套循環(huán)。當(dāng)主線程遇到外層循環(huán)時(shí)，派生出一組線程，這些線程同時(shí)繼續(xù)向內(nèi)層循環(huán)執(zhí)行下一層并行，這時(shí)派生出的線程和原來的主線程成為針對內(nèi)存循環(huán)的主線程，繼續(xù)按照之前的方式派生新的線程從而并行執(zhí)行內(nèi)層循環(huán)，當(dāng)內(nèi)層循環(huán)執(zhí)行結(jié)束時(shí)，這些線程逐級掛起，最終回到最初的主線程上。嵌套循環(huán)并行執(zhí)行線程運(yùn)行過程如圖8所示，m為實(shí)際計(jì)算機(jī)核數(shù)。

圖8 嵌套循環(huán)并行線程運(yùn)行圖Fig.8 Threads of nested loop executing in parallel

b.非嵌套方案。另外一種方法即類似于小波分解并行執(zhí)行過程。對于嵌套循環(huán)的程序，可以選擇只對一層循環(huán)作并行處理，由于小波包分解中外層循環(huán)次數(shù)不多，而內(nèi)層循環(huán)復(fù)雜度較高，本文選擇只對內(nèi)層循環(huán)過程作并行處理，即利用多核處理器同時(shí)對不同尺度的系數(shù)做卷積和二抽取過程。非嵌套循環(huán)并行執(zhí)行線程運(yùn)行狀態(tài)如圖9所示。

圖9 非嵌套循環(huán)并行線程運(yùn)行圖Fig.9 Threads of non-nested loop executing in parallel

3 數(shù)據(jù)壓縮分析實(shí)驗(yàn)

3.1 實(shí)驗(yàn)基本過程

為了檢測并行小波包在數(shù)據(jù)壓縮中的性能表現(xiàn)，本文對數(shù)據(jù)進(jìn)行并行小波包壓縮解壓處理，由于小波包變換可以使信號經(jīng)過小波包變換后能量集中在少數(shù)的系數(shù)上，并且這些系數(shù)的能量能夠遠(yuǎn)大于其他多數(shù)小波包的系數(shù)，因此可采用閾值方法，現(xiàn)有的文獻(xiàn)閾值選取包括全局閾值與局部閾值，考慮到全局閾值能降低計(jì)算的復(fù)雜度，而本文分析重點(diǎn)在于研究算法的并行性，因此選取全局閾值壓縮。

數(shù)據(jù)壓縮方法采用閾值壓縮算法，閾值處理后需要記錄非零小波系數(shù)的位置以及系數(shù)的取值，此過程采用文獻(xiàn)[22]所提出的位圖壓縮法。

小波分解的最大級數(shù)J定義[23]如下：

其中，ff是信號的基礎(chǔ)頻率，fs是信號的采樣頻率。

實(shí)驗(yàn)主要步驟如下。

a.載入信號后，對信號作并行小波包分解；分解層數(shù)根據(jù)式（1）得出，本實(shí)驗(yàn)分解層數(shù)為3層，分解過程采用第2節(jié)提出的小波包并行分解方案。

b.將分解后的小波系數(shù)，經(jīng)過全局閾值處理，利用位圖壓縮算法將剩余后的小波系數(shù)依次排列，完成壓縮過程。解壓過程根據(jù)位圖壓縮法還原各系數(shù)原始位置，對其余位置進(jìn)行補(bǔ)零操作。

c.采用第2節(jié)提出的并行小波包重構(gòu)方案，對解壓后系數(shù)作小波包重構(gòu)，還原原始數(shù)據(jù)。

實(shí)驗(yàn)主要用于檢驗(yàn)并行小波包的計(jì)算效率，在保證分解重構(gòu)的結(jié)果與串行計(jì)算結(jié)果相一致的前提下進(jìn)行。

實(shí)驗(yàn)時(shí)間記錄從主線程執(zhí)行并行小波包分解開始，直到小波包重構(gòu)結(jié)束。由于實(shí)驗(yàn)針對小波包分解壓縮重構(gòu)的并行過程，因此不記錄載入信號與寫出信號的時(shí)間。

仿真信號源采用的是：

3.2 實(shí)驗(yàn)環(huán)境

本文采用的實(shí)驗(yàn)平臺為：處理器為Intel（R）Core（TM）2 Duo CPU T5750，內(nèi)存為 3 GB。 所有數(shù)據(jù)為5次計(jì)算時(shí)間平均值。

3.3 實(shí)驗(yàn)方案

為比較并行算法性能，采用本文所提出基于Pthreads及OpenMP 2種并行編程環(huán)境下的小波包并行算法與串行小波包算法分別進(jìn)行數(shù)據(jù)壓縮。實(shí)驗(yàn)在 210、215、216、220、221這 5 種數(shù)據(jù)量下，采用 db4 小波包，將信號分解為3層。實(shí)驗(yàn)結(jié)果見表1、2。

表1 小波包Pthreads數(shù)據(jù)壓縮耗時(shí)Tab.1 Time consumption of data compression by wavelet packet on Pthreads

表2 小波包OpenMP數(shù)據(jù)壓縮耗時(shí)Tab.2 Time consumption of data compression by wavelet packet on OpenMP

為了顯示并行算法相對串行算法的性能改變，采用加速比描述：加速比=串行算法計(jì)算時(shí)間／并行算法計(jì)算時(shí)間。

表3中 S1、S2、S3分別表示 Pthreads環(huán)境下小波包并行壓縮解壓加速比、OpenMP環(huán)境下小波包嵌套并行壓縮解壓加速比、OpenMP環(huán)境下小波包非嵌套并行壓縮解壓加速比。

表3 加速比比較Tab.3 Comparison of speedup

3.4 實(shí)驗(yàn)結(jié)論

a.從實(shí)驗(yàn)可以看出基于Pthreads與OpenMP的并行小波包算法，在計(jì)算速度上有著較大的優(yōu)勢，這種性能的提升隨著數(shù)據(jù)量的增大而增大。相對于串行算法有著接近2倍的速度提升，甚至在用OpenMP進(jìn)行并行處理時(shí)出現(xiàn)超線性加速比的情況。

b.在數(shù)據(jù)量較少時(shí)，由于并行計(jì)算不可避免的開銷要比并行效果帶來的節(jié)約時(shí)間要大，并行計(jì)算的性能不如串行計(jì)算，在數(shù)據(jù)量為215時(shí)，各種方法的加速比表現(xiàn)較為一致，且都大于1，說明在此數(shù)據(jù)量附近，并行帶來的節(jié)約只略大于并行開銷，一旦計(jì)算機(jī)有其他程序的任務(wù)，就會嚴(yán)重影響此時(shí)加速比的表現(xiàn)。這一問題將隨著海量數(shù)據(jù)處理的發(fā)展而逐漸被忽略。

c.隨著數(shù)據(jù)量的增大，基于OpenMP的并行小波包加速比比基于Pthreads的并行小波包要好，這是由于Pthreads是直接控制線程調(diào)度，可以由程序員自主控制程序的底層運(yùn)行，一旦代碼優(yōu)化效果好，就可以達(dá)到較高的加速比，反之，由于優(yōu)化存在難度，在代碼復(fù)雜度較高的情況下，可能達(dá)不到OpenMP的并行效果。

d.嵌套循環(huán)的并行性能比非嵌套循環(huán)較弱，這是由于嵌套使能開啟后，內(nèi)部循環(huán)仍舊按照之前的方式派生線程，由于外層循環(huán)根據(jù)處理器個(gè)數(shù)已經(jīng)產(chǎn)生最大可并行執(zhí)行的線程，因此，在內(nèi)層循環(huán)中所派生出的新線程與其主線程并沒有同時(shí)執(zhí)行程序而是交替進(jìn)行。此外創(chuàng)建和退出線程的過程又帶來一定的時(shí)間消耗，使得性能下降。但是在處理器個(gè)數(shù)大于2時(shí)，通過動(dòng)態(tài)分配線程，改變外層和內(nèi)層線程分配方式，可以使嵌套并行循環(huán)的性能更優(yōu)。

4 結(jié)語

隨著單機(jī)多核處理器的廣泛使用，利用多核技術(shù)實(shí)現(xiàn)并行小波包計(jì)算可以大幅提高處理數(shù)據(jù)的速度，進(jìn)一步促進(jìn)小波包變換在電力系統(tǒng)中的實(shí)際應(yīng)用，使其可以滿足由于高速電力發(fā)展帶來的海量數(shù)據(jù)處理實(shí)時(shí)性的要求，保證大電力系統(tǒng)更安全、穩(wěn)定、經(jīng)濟(jì)運(yùn)行。