一种改进的双库MP信号稀疏分解方法

Two dictionaries Matching Pursuit for sparse decomposition of signalsPeng XuCenter of Neuroinformatics, School of Life Scienceand Technology, University of Electronic Science andTechnology of China, ChengDu, 610054, China leisure_xp@Dezhong Yao Center of Neuroinformatics, School of Life Science and Technology, University of Electronic Science and Technology of China, ChengDu, 610054, China dyao@AbstractSignal may be decomposed sparsely and power-focally in an over-complete dictionary with Matching Pursuit (MP). In this paper, proposed is a modified MP method named two dictionaries MP to decompose signal more sparsely. In the iteration procedure of the two dictionaries MP , the over-complete dictionary is classified into two separate dictionaries with the selected and unselected atoms, and in each iteration, the algorithm was designed to have more chances than the original MP to choose the atom in the selected atom dictionary as the optimal atom by a simulate annealing threshold function, thus the algorithm avails for a more sparse decomposition. The decomposition results for a cosine-modulated exponential signal and an actual speech signal showed that the proposed two dictionaries MP could decompose signal more sparsely.Keywords : Sparsity; Over-complete dictionary; Matching pursuit; Simulate annealing threshold function1. IntroductionHunting for a compact and effective signal expression is one of the research focuses in the signal processing (SP) because the more compact expression the signal has, the higher ratio compression the signal can attain [1-4]. Sparse component analysis (SCA), developed recently, is an effective signal processing method for signal compact representation. The principle of SCA could be briefly stated as the following. Once a certain signal only has several non-null entries in time or other transformation domains, it means that the power of signal is focally concentrated and this signal could be compactly represented in the corresponding domain with only a few non-null entries [2-5]. Sparsity is adopted to measure how compact a signal is. Signal sparsity could be measured with )1(≤p l p norm or the decaying-to-null speed of the sorted absolute signal entries [3-4, 6]. Most of the current SCA algorithms are to decompose signals sparsely in an over-complete dictionary because the decomposition in the over-complete dictionary is not unique and has great freedom for the basis selection, moreover when the over-complete dictionary is large enough, the signal must have the sparse representation in the over-complete dictionary if sparse constraints are emphasized [3-6]. The construction of over-complete dictionary is adapted with different applications and in this paper, the Symmlet wavelet was taken to construct the over-complete dictionary [3-4]. There are many approaches to the SCA decomposition, among which Matching Pursuit (MP) developed by Mallat and Zhang is a popularly used one [7-8]. MP decomposes signal into sparse representation iteratively in the over-complete dictionary and in thedecomposition procedure, the optimal basis (atom) for the signal representation could be adaptively selected according to the characteristic of signal.In this paper, to get a more sparse decomposition, the over-complete dictionary is classified into two separate dictionaries with the selected (old) and unselected (new) atoms, and similar to the selection of optimal atom in the original MP , two local optimal atoms in the old and new atom dictionaries will be firstly taken out. If the local optimal atom in the old-atom dictionary has larger coefficient (inner product) than that in the new-atom dictionary does, the former one will be considered as the global optimal atom for this iteration. Otherwise, the relative error (RE) between two residuals, which corresponds to these two local optimal atoms, is taken as a measure to decide the global optimal atom, if RE is smaller than a supposed threshold, the old one will be chosen, else the one in new atom dictionary will be considered, thus the algorithm will have more chances than the original MP to select the atom in the old-atom dictionary as the optimal atom, which avails for a more sparse decomposition.The two dictionaries MP (TDMP) was introduced in detail in Section 2 and a comparative study to original MP was conducted in Section 3. Discussions and conclusions concluded this paper in Section 4.2. Methods2.1 Sparsity of signalFor a signal denoted by X , if there are only a few non-null entries in its coefficients of representation and the number of non-null entries is relatively small compared with the length of signal, the signal can be supposed to be sparse. The sparseness of signal X could be measured with )1(≤p l p norm as follows,p n j pj p x X 11)||(∑== (1)where 10≤≤p , and n is the signal dimension of X . The sparser a signal is, the smaller )1(≤p l p norm is. At the same time, the sparseness can be expressed in another intuitionistic way. Let X and Y be two different signal arrays, if these two arrays are sorted in the descending order by their absolute values of coefficients, the one with faster decaying-to-null speed is sparser than the other one. As shown in Fig.1, X is sparser than Y . [3-4]Fig.1. An intuitionistic way to compare the sparseness of signals2.2 Construction of wavelet over-complete dictionaryIn signal processing, the usually adopted transform dictionary (matrix) is orthogonal and the decomposition is unique [4]. The current SCA methods are mainly based on the learning in a over-complete dictionary to decompose a signal sparsely. The over-complete dictionary can be derived from the complete dictionary by over-sampling. The difference between the over-complete dictionary and the complete dictionary is that the basis in the complete dictionary is orthogonal to each other, while the orthogonality may not be true in an over-complete dictionary. According to different problems, different over-complete dictionaries avail for various situations [2-4, 6]. In our work, we adopted the over-complete dictionary based on the wavelet. A wavelet basis function can be expressed as,)()(,sb t t W b s −Ψ= (2) where s,b are the scale and translation parameters, respectively. If over-sampling is implemented for scale and translation parameters in the wavelet function, a wavelet over-complete dictionary could be expanded from the traditional orthogonal dictionary [3-4], and in our work, the Symmlet wavelet was taken to construct the over-complete dictionary. Certainly, a large over-complete dictionary is good to sparse decomposition, but the decomposition time will increase a lot accordingly. To consider the decomposition effect and time, in our work, the number of atoms consisting of the over-complete dictionary was 6 times of signal length. In Fig.2, three S ymmlet wavelet atoms at scale s =64 with different translations b, 85, 95 and 105, are shown. It shows that these three atoms in the over-complete dictionary are not orthogonal.Fig.2. Three atoms in the Symmlet wavelet over-complete dictionary. The scales are all 64 and from leftto right, the translations are 85, 95 and 105, respectively. 2.3 Signal decomposed with Matching PursuitMP was developed by Mallat and Zhang to decompose signal into sparse representation iteratively in the dictionary Φ[7]. In each iteration, a certain atom with the maximum inner product between itself and the residual signal is selected as the optimal one. After several iterations, when the convergence condition is satisfied, a sparse decomposition of signal in the corresponding dictionary is achieved. The iterative MP algorithm can be briefly stated as follows.(1) Initialization: Set S R S k ===)0()0(,0,1, where S is the signal to be decomposed, R isthe residual signal during the iterations; the superscript is the current iteration number; initializing the stopping error ε with a small positive number; initializing the coefficients adders for all the atoms with 0.(2) Atom selection: Calculate inner product between current residual signal and each atom inthe over-complete dictionary, and find the atom with maximum inner product as follows: ><=−j k jk R φα,max )1(, where j φis the jth atom in the dictionary and <·> is the inner product operator.(3) Update: Accumulate the current maximum coefficient k α to the adder of current optimalatom; add the weight of the selected atom to the decomposed signal )(k S and remove the corresponding contribution of this atom from the residual signal )(k Ras,)()()1()(,k k k k k k S S R S S −=+=−φα(4) Convergence judgment: If ε≤||||)(k R , stop the iterative procedure; else 1+=k k ,jump to step 2 and go on.When algorithm converges, signal will be sparsely decomposed into a linear combination of some atoms in the over-complete dictionary.2.4 Two dictionaries Matching PursuitConsidering the sparsity of signal decomposition, the fewer the selected atoms are during the decomposition, the more sparsely the signal will be decomposed [3, 6]. During the MP decomposition procedure, the selected optimal atoms can be divided into two kinds: one has been considered as the optimal atom in the previous iterations, and the other belongs to such kind that has not been selected in the previous iterations. To be consistent with the above two kinds of atoms, the over-complete dictionary consists of two kinds of separate dictionaries, one is constructed by the “old atoms ” selected as the optimal atoms in iterations, the other includes those “new atoms ” never considered as the optimal atoms in iterations. For convenience to describe the algorithm, the former kind of dictionary is named the old-atom dictionary and the latter one is the new-atom dictionary. Accordingly, at the beginning, all atoms belong to the new-atom dictionary. In the previous iterations of MP decomposition procedure, most of the selected optimal atoms are in the new-atom dictionary and with the iterations on, the old-atom dictionary will be expanded. When the old-atom dictionary is large enough, the selected optimal atom may be either in the old-atom dictionary or in the new-atom dictionary. From the view of sparsity, the selection of the new atom does not avail for the sparsity of decomposition and to select the atoms in the old-atom dictionary as possible is expected for a sparse decomposition. Thus to get a more sparse decomposition, we propose that the criteria for optimal atom selection can be adjusted as follows.The over-complete dictionary is divided into the old-atom dictionary old Φand thenew-atom dictionary new Φ, then the inner product between the residual signal and each atom inthese two dictionaries is calculated respectively, thus two local maximum coefficients of inner product, old c and new c , can be achieved in the old-atom dictionary and new-atom dictionary respectively. Accordingly, the two local optimal atoms, old φ and new φcorresponding to old c and new c , can be determined in the old and new atom dictionaries too.If ||||new old c c ≥, the atom old φ is the global optimal atom for the current iteration and update the searching variables with old opt old opt c c φφ==,and opt opt c R R φ−=. Finally, accumulate current coefficient old c to the adder of the selected old atom old φ.If ||||new old c c <, the global optimal atom in the kth iteration will be decided with the following steps.(1) Calculate the two residual signals in the old and new atom dictionaries as:new new k new old old k old c R R c R R φφ−=−=−−11,, with 1−k R the residualsignal in the prior iteration.(2) Calculate the relative error r between the residual signals in the old and new atomdictionaries as:||R ||||R R ||r new new old −= (3) Decide the global optimal atom by the given threshold function T. If T r ≤, the old atomold φ is considered as the global optimal atom and update the variables with old opt old opt c c φφ==,; If T r >, the new atom new φis newly selected as the global optimal atom and update the variables with new opt new opt c c φφ==,, insert this new atom into the old-atom dictionary and remove the new atom from the new-atom dictionary. After the optimal atom is decided, the corresponding coefficients opt c is accumulated to the adder of the optimal atom.(4) Update the residual signal: opt opt 1k k c R R φ−=−.The main difference between the original MP and TDMP lies in that the selection of the global optimal atom is decided in the divided new and old atom dictionaries by a given threshold function T . If T=0, our method reduces to the original MP . To be convergent and stable, T is a descending function with iteration steps and in this paper, the annealing function used in the popular simulate annealing algorithm is adopted [9],N /1k 0T )k (T α×= (3) where 17.0<≤α, 0T is the initialization temperature and set to smaller than 1 in our work, k is the current iteration number, N is the annealing speed factor. With the increasing of iteration, T is close to 0, thus in the final iteration steps, the effect of threshold function on the atom selection is lowered and the decomposition coefficients are small, so the decompositions in the final iterations mainly correspond to the detailed information of signal and the convergence can beguaranteed too.3. Results3.1 Numeric comparative studyIn this section, we took a cosine-modulated exponential signal as the test signal to numerically compare TDMP with MP . The adopted test signal was,)64/2cos(985.0)(|128|x x f x π−= (4) where 2561≤≤x and the test signal was of 256 sampling points length. The over-complete dictionary with 1536 )62561536(×=atoms was constructed by the Symmlet wavelet. The Symmlet wavelet, original MP and TDMP were taken to decompose signal. In original MP and TDMP , the terminating error was both 1.0E-6. The simulate annealing function was defined as 4.0k 0.935*0.9=T(k), with 9.00=T , N =2.5 and 935.0=α. The decomposition coefficients of these three methods were descendingly sorted by their absolute values and the curves of sorted coefficients were shown in Fig.3. Among these three methods, the wavelet coefficients descended slowly and the MP coefficients decayed faster than those of wavelet, whereas, TDMP had the steepest descending coefficients curve.Fig.3. The sorted decomposition coefficients of wavelet, MP and modified MP .To compare the difference of atom selection in the decomposition procedure between the original MP and TDMP, we showed the index of selected atom in each iteration step in Fig.4.(a)(b)Fig.4. The indexes of selected atoms during the decomposition iterations. (a) original MP; (b) TDMP; the green circles are the atoms in the old-atom dictionary; the blue asterisk are the atoms in the new-atom dictionary.The signal was reconstructed with different number of the atoms (coefficients) and the relative error (RE) between the recovered signal and the original signal was calculated. In this paper, the reconstruction with n atoms (coefficients) was referred to that the former n large coefficients and the corresponding atoms were taken to recover the signal. The reconstruction results were shown in Fig.5.(a) 5 atoms(b) 10 atoms(c) 15 atomsFig.5. Signal recovered with different number of atoms (coefficients). The left column was the results recovered from wavelet coefficients; the middle column was the results recovered from the original MP coefficients; the right column was the results recovered from the modified MP coefficients. The green signal was the error signal, the blue one was the reconstructed signal and the red one was the original test signal.3.2 Decomposition of an actual speech signalThe speech signal usually needed to be highly compressed. In this section, a speech signal with 512 sampling points was decomposed and recovered with Symmlet wavelet, original MP and TDMP , respectively. The over-complete dictionary with 3072 )65123072(×= atoms was constructed with Symmlet wavelet and the simulate annealing threshold function was similar to that in Section 3.1. The signals recovered with 40 atoms (coefficients) for these 3 methods wereshown in Fig.6.Fig.6. The speech signal reconstructed with 40 atoms (coefficients). The green signal was the errorsignal and the blue one was the reconstructed signal.4. Discussion and ConclusionsAccording to the decaying speed of the sorted coefficients shown in Fig.3, the decomposition coefficients of wavelet were much scattered compared to other two methods and the coefficients of original MP and TDMP were relatively sparse and power concentrated. Compared to the original MP , TDMP could decompose signal more sparsely.As shown in Fig.4, the selection of optimal atom during the decomposition was largely different between the original MP and TDMP . During the decomposition procedure of the original MP , many of the selected optimal atoms were in the new-atom dictionary, which went against the decomposition sparsity. Whereas, in the decomposition procedure of TDMP , the atoms in the old-atom dictionary were selected as the optimal atoms in most of the iterations, except for the first and final several iterations. In the first several iterations, though the annealing temperatureinvolved in the threshold function was high, the new-atom dictionary was relatively large and the old-atom dictionary was very small, thus most of the selected optimal atoms in these iterations belonged to the new-atom dictionary. With the iterations on, the annealing temperature was descending and the old-atom dictionary was expanded fast. When the old-atom dictionary was large enough compared to the shrinking new-atom dictionary, though the annealing temperature was relatively low, most of the selected optimal atoms were in the old-atom dictionary in the middle stage of the decomposition procedure. In the final iterations before convergence, the simulate annealing threshold function was close to 0 and it has little effect on the atom selection, so some new atoms were selected. Furthermore, in these final iterations, the coefficients of optimal atom were small and they mainly represented the detail of signals, thus the convergence of algorithm could be guaranteed.In Fig.5, the reconstruction differences among these three methods were obvious. Based on the sparse decomposition in the over-complete dictionary, only 10 atoms were needed to recover the signal with small REs, 3.9% for the original MP and 0.5% for TDMP, while with 10 coefficients of wavelet to reconstruct, the RE was nearly 80%. Especially, when 15 atoms (coefficients) were taken for signal reconstruction, the RE (0.09%) of TDMP was much smaller than the REs of the MP (0.3%) and the wavelet (72.6%). So among these three decomposition methods, TDMP could decompose signals into a more sparse and power concentrated representation, and could recover the signal with high quality. In the decomposition test for an actual speech signal, a good recovery signal was achieved with TDMP too.The above results showed that a signal could be decomposed more sparsely with TDMP than with the original MP because during the decomposition procedure of TDMP, the selection of atom was affected by an introduced annealing strategy threshold function to select the old atoms availing for sparse decomposition more easily.MP is very time consuming when the greedy exhaustive search in the whole huge over-complete dictionary is adopted. How to improve the efficiency of sparse decomposition is still a challenging problem. With SCA decomposition in the over-complete dictionary, the power of signal will be concentrated on several atoms and by suitable choice of atoms for signal representation, the noise can be effectively removed [3-4, 6]. Furthermore, the visual processing system are proved to encode the visual information sparsely [10-11], so if an visual image or signal is decomposed and processed with SCA, a sparse representation of visual signal consistent with the visual physiology can be expected.References[1] A.N. Skodras, C.A. Christopoulos, T. Ebrahimi, JPEG2000: The upcoming still imagecompression standard, Pattern recognition letters, 22 (October 2001) 1337-1345.[2] R. Adelino, F.D. Silva, Atomic decomposition with evolutionary pursuit, Digital signalprocessing, 13 (April 2003) 317-337.[3] S. Chen, D.L. Donoho, M. Saunders, Atomic Decomposition by Basis Pursuit, SIAM J.SciComp, 20(1) (1998) 33-61.[4] X.M. Huo, Sparse Image Representation via Combined Transforms, PHD Paper of Stanforduniversity, 1999, available online at: /xiaoming.[5] P. Bofill, M. Zibulevsky, Underdetermined blind source separation using sparse representations,Signal Processing, 81 (November 2001) 2353-2362.[6] D.L. Donoho, Sparse Components of Images and Optimal Atomic Decompositions, 1998,available online at: http://www. /donoho/Reports.[7] Mallat S, Zhang Z. Matching Pursuit in a time-frequency dictionary. IEEE Transactions onSignal Processing, 41 (December 1993), 3397-3415.[8] M.R., McClure, L. Carin, Matching pursuits with a wave-based dictionary, IEEE Transactionson Signal Processing, 45(12) (December 1997) 2912-2927.[9] L. Ingber, V ery fast simulated reannealing, Math Comput Modeling, 12 (1989) 967-973.[10] B.A. Olshausen, D.J. Field, Sparse coding with an overcomplete basis set: a strategyemployed by V1, Vision Research, 37 (December 1997) 3311-3325.[11] D. Marr, Vision: A computational investigation into the human representation and processingof visual information, Freeman New York, 1982.。