 Home |
|  Home |
566KB
|
THE WAVELET IMAGE DENOISING FAST ALGORITHM STUDY BASED ON DSP Zhicheng Hao Jiyin Zhao Xiaoling Wang Yu Liu College of Communication Engineering, Jilin University, Changchun, P. R. China Abstract: The speed of wavelet transform algorithm have restrained its application all the times. To this, the paper propose a fast wavelet denoising algorithm, which was identified in mathematics, and the algorithm step was given in this paper. Through analysis of Mallat tower algorithm principle, some step was improved, and fast FFT and long sequence signal volume fast algorithm was used in this paper, so that operation amount can be reduced greatly. DSP emulation experimental result indicated, compared with denoising algorithm using direct wavelet transform, this method can improve the speed of operation greatly. As to denoising effect, this method is superior to other several kinds of methods not only at the index of denoising performance, but also in visual effect. Key words: Long sequence volume, Mallat algorithm, PCT denoising algorithm, PSNR , DSP Introduction: Wavelet theory is a new application theory developing recently. Due to its good time-frequency local characteristics, wavelet analysis have been applied extensively in many fields, such as signal analysis, image processing, speech synthesis, trouble diagnosing, geological prospecting, etc. But because of the huge operation amount of direct discrete wavelet transform and its complexity, its application in actual system is hindered greatly, so the further research to discrete wavelet transform fast algorithm is necessary. Long sequence signal fast volume: For image signal, signal x(n) are long and filter h(n) are short generally. If computing correlation or volume directly, h(n) must be filled with many zeros, that is too wasteful. So x(n) can be divided into segments xi(n) which each length is similar to h(n). Supposing length of h(n) is L, so xi(n)~ i=0,1,...N/L-1, N is the length of x(n). Computing the volume yi(n)=xi(n)*h(n), according to the recurrence formula, yi(n) are combined together, then final output y(n) can be obtained. The formula is derived as follows. M Suppose, the lengths of h(n) and x(n) is M and N respectively, and N >>M. Volume 0 l , n=0,1,2,..., N+M-2. Divide x(n) into m segments, the length of each segment is L, so N=m×L, suppose p=m-1, In order to analyse the question, specially explain and derive the course for example. 1 ) ( l n h l x n y For example: the length of x(n) is N=10, the length of h(n) is M=4, suppose L=5, so m=2~p=1Ñ+L-1=13, xi(n), i=0,1 is the segments. In x0(n): x0(0)=x(0)~ x0(1)=x(1)~ x0(2)=x(2)~ x0(3)=x(3)~ x0(4)=x(4)~ In x1(n): x1(0)=x(5)~ x1(1)=x(6)~ x1(2)=x(7)~ x1(3)=x(8)~ x1(4)=x(9)~ y0(n), y1(n) (0 n L+M-2) is the volume of x ≤ ≤ 0(n), x1(n) and h(n) respectively, which can be indicated: y0(n)=x0(n)*h(n) , y1(n)=x1(n)*h(n). And y(n) is the final output which be indicated: y (n)=x (n)*h(n). Through the calculation, we conclude: y(0)=y0(0)ŷ(1)=y0(1)ŷ(2)=y0(2)+y1(5)ŷ(3)=y0(3)+y1(6)ŷ(4)=y0(4)+y1(7)ŷ(5)=y1(0)ŷ(6)=y1(1)ŷ(7)=y1(2 )ŷ(8)=y1(3)ŷ(9)=y1(4)~ y(10)=y0(5)ŷ(11)=y0(6)ŷ(12)=y0(7)~ Through summing up, the recurrence formula can be get: When k~0: y(0)=y0(0), and y(L-1)=y0(L-1); When 1 : y(kL)=y p k ≤ ≤ k(0)+yk-1(L), and, y(kL+M-2)=yk(M-2)+yk-1(L+M-2); = - ∑- = ) ( ( ) y(kL+M-1)= yk(M-1), and, y(kL+L-1)= yk(L-1); When k=p=N/L-1: y(kL+L)=yk(L), and, y(kL+L+M-2)=yk(L+M-2); While calculating the volume, a unavoidable thing is to use FFT, here a fast FFT algorithm can adopt in order to improve the transform speed.[1] principle of Mallat algorithm: The basic thought of Mallat tower algorithm is as follows: If we have already calculated out the discrete approaching coefficient and the discrete detail coefficient d of the signal under resolution 2 , the next level wavelet coefficients and can be calculated as follows: ) (n c j j ) (n ) (n f j - ) 1 n c j+ ( 1 n d j+ ) ( n c ) 2 ( ) ( ) ( j n m h m c + - 1 = ∑Z m j ∈ 1 (1) and n d ) 2 ( ) ( ) ( j n m g m c + - = ∑Z m j ∈ n c j ' ∑∈ + j - 1 n h n c m n h m c ) ( = ) ( )] ( [ = - - ) ( ∗ ) ( Z m j n d j ' ∑∈ + (2) h, g are both wavelet filters of real coefficient whose lengths are both L(If the lengths of them are not equal, the short one can be filled with zeros), h is low-pass filter and g is high-pass one. By formula (1) and (2) we can see that and d are obtained through dyaddown of c j - 1 n g n c m n g m c ) ( = ) ( )] ( [ = - - ) ( ∗ ) ( Z m j ' respectively. Formula (1) and (2) are the Mallat decomposition algorithm formula. Utilizing the course adverse to above, wavelet reconstruction algorithm formula is obtained as follows: ∑1 n c j+ ) ( j+ 1 n ) ( j+ ' and ) 1 n ) ( 1 n d j+ ( = ) ( j (3) The figure 1 show the process of Mallat algorithm wavelet decomposition and reconstruction.[2] g(-n) 2 d h(-n) 1(n) ∑+ + - +- 1 m d n m g m c n m h n c j j ) ( ) 2 ( ) ( 1 ) 2 ( ' (n) d1 f (n) c (n) 0 c (n) 1 (n) ' d (n) 2 2 d ' c 2 (n) (n) c2 1 g(-n) 2 (n) c' 2 2 h(-n) (a) 2h(-n) (n) d1 2 g(-n) 1(n) d ' d 2 (n) '(n) d (n) c1 2 h(-n) 2 2g(-n) + f (n) (n) 0 c + c (n) ' 1 2 c2 (n) ' c (n) (b)Fig.1 Mallat algorithm structrure Improved Mallat algorithm: If calculating c ) 1 n ' , d ) ' directly using Mallat algorithm, each need NL multiplication once , 2NL multiplication once altogether. In order to analyse expediently, suppose L =N, therefore 2N2 multiplication is need if calculate directly, the bigger N is, the bigger the calculating amount is, unfavorable to the real-time processing of the signal. Provide improved analysis and derive the course as follows.[3] First of all, The two real coefficients filters h(n), g(n) are formed into a complex sequence x(n)=h(n)-ig(n), and calculate its Fourier transform value X(k)=H(k)-iG(k) with fast FFT algorithm. cycle continuation to h(n),g(n), ( 1 n ( form a cycle sequence ) ( ~ n x = x ((n))N= ) n - h . and x (n)= ) ) ( ~ n g i ( ~ ( ~ n x n c ) ( RN(n), while n=0,1,...,N-1, RN(n)=1; else RN(n)=0. Carry on the same cycle continuation to c , namely ) n ( n c )) ) ( ~ ) ( 0 0 n R n c n N = c . Seeing that the length of the signal and length of the wavelet filters are all N, and that is N too after carrying on cycle continuation. So the formula (2) is turned into: 0 N ~0 = , ) (( ( 0 n m h m c n c 1 ) = ' ) ( N m ∑- = 1 0 ) ( ~ ( ~ 0 = - ) N m ∑- = 1 0 n m h m c ) ( ( ~ - 0 n m g m c n d (4) 1 ) = ' ) ( N m ∑- = 1 0 ) ( ~ ( ~ 0 = - ) N m ∑- = 1 0 n m g m c ) ( ( ~ - 0 here, n=0,1,...,N-1.and when , 0 - ≤ ≤ N m ) ( ) ( ~ 10 m c m c = .So = ' ) ( 0 1 k C .Suppose ∑N n W n c n m r - -11 ) ' ( nk N = , well them r m n - = , m r ≤ N m ≤ - - 1 ( m N ) . So - m - N k C 1 ) = ' ) ( ∑∑1 0 0 W r h m c ) ( ) ( ~ k r m N W m c ( - = ) ∑ ∑W r h ) ( ~ rk N 1 0 0 - ( ) ( ~ N m r m + - = = 1 N m r m nk N + - = = 1 - - - - 1 1 1 1 N N N N = ∑∑∑∑W m c W r h (5) rk N W m c - = 0 ) ) ( W r h ) ( rk N - ( r m nk N = = 0 0 r m nk N = = 0 0 0 Because , c is real, we can get: ) (r h ) ( 0 m * k C 1 k = ' ∑) ( 0 ∑N r - = 1 0 W r h ) ( rk N = N m - = 1 0 C k H W m c ) ( * 0 nk N ) ( ) ( (6a) According to this principle, we can get too: ' (6b) So we can get the following: C k G k D = 1 k ) ( 0 * ) ( ) ( = = C ~ (7) G k iC k H k C k iG k H k C k X k C + 0 k 0 0 0 ) ( * ) ( )[ ( * ) ( + * )] ( = ) ( * ) ( ) ( * ) ( ) ( ) ( 1 1 k + ' k ' D i Carry on adverse Fourier transform to formula (7), we can get ) ( n c ' ) ( 1 1 n d i + ' , and carry on dyaddown to the result, ) ( then extract its real and image part, the discrete approaching coefficient c and the discrete detail coefficient of the signal can be get. The step of improved Mallat deconstruction algorithm is as follows: 1 n 1 n d ) ) ( (n f (1) Select wavelet filters h(n), g (n), and they form complex sequence x(n)=h(n)-ig(n); Calculate out its Fourier transform X(k); And then compute its conjugation X*(k). (2) The signal is f(n), divide it into segments fi(n), then utilizing N/2 point FFT calculate Fi(k)=DFT{fi(n)}. (3) Calculate out the product of Fi(k) and X*(k), express it for Yi(k). (4) Calculate Y(k) using the long sequence volume fast algorithm. (5) Calculate IDFT{Y(k)}= =y(n) using N/2 point IFFT. n c ' ) ( 1 1 n d i + ' ) ( ) ( (6) Carry on dyaddown to y(n), then extract its real and image part, the discrete approaching coefficient and the discrete detail coefficient d of the signal can be get. 1 n c ) (n f The step of improved Mallat reconstruction algorithm under scalê+1 is as follows: (1) Carry on dyadup to approaching coefficient and detail coefficient , get c and . 1 n ) ( 1 n c j+ ) ( 1 n d j+ ) ( ' j+ 1 n ) ( 1 n d j+ ' ) ( (2) Form complex sequence y(n)= + 1 n c j+ ' ) ) ( 1 n d j+ ' , calculate out its Fourier value Y(k) using fast FFT. ( (3) Carry on the fast volume of y(n) and x(n), (namely calculate the product of Y(k) and X*(k), then calculate its IFFT). (4) Extract the real part of the volume value, it is the reconstructive approaching coefficient under scale j. ) (n c j Denoising algorithm: Threshold method and proportion-shrinking method are both commonly used method in wavelet denoising field. The former zeroize the wavelet coefficients smaller than threshold, save the bigger ones than threshold; the latter estimate the processed coefficients with wavelet coefficients multiplied by a proportional coefficients. Both of them have their own merits and demerits: Threshold method zeroize many smaller coefficients, it will bring up much facility for follow-up image processing, such as image compression. But meanwhile it will delete some useful signal wavelet coefficients so that reconstructive image produce distortion and twist. Proportion-shrinking method keep all of the image wavelet coefficients, make the reconstructive image to be undistorted, but because of it the phenomenon of the burr will come forth. Seeing that, this paper propose a new denoising algorithm - PCT algorithm (Proportion-shrinking Combined with Threshold). The essence of PCT algorithm is putting threshold policy into proportion-shrinking in order to obtain lubricous rebuilding image signal. Suppose the observed signal is: j i n j i s j i x ⋅ ] , [ = , ] , [ + ] , [ j i 2 , 1 , = ⋅ ⋅ ⋅ ⋅ ⋅ N is original image, (i, j) is pel location, is the Gauss white noise added to that point. Suppose , indicate the wavelet coefficients of , under scale k respectively, so the denoised wavelet coefficients can be get as follows: ) , ( j i s , ( j i Sk ) , ( j i n ) ) , ( j i X k ) , ( j i s ) , ( j i x σ ˄ ) , ( ˄ 2 ) , ( j i k (8) Formula (8) is the function of PCT algorithm. Here, indicates the estimation of noise variance, indicates the estimation of signal wavelet coefficient variance, and their calculation method are introduced in reference [4]. 2 j i S = k k j i X j i X ˄ ) , ( ) , ( ≥ λ σ , 2 , 2 k s k s ) , ( j i + σ ˄ 0 j i X ) , ( < λ n k σ ) ˄ n ik s σ ˄ , 2 j , ( λ is wavelet threshold, which is chosen locally adaptive threshold, and its calculation method is introduced in reference [5]. According to formula (8), if is bigger than threshold ) , ( j i X k λ , the function is similar to that of proportion-shrinking algorithm; if is smaller than ) , j i ( X k λ , the function is similar to that of soft-thresholding method. Results: We combined PCT denoising algorithm with improved Mallat algorithm, and experiment it on DSP emulation system. Choose the standard static Woman image with 256*256 size as experimental object. Daubechies wavelet filter was selected, whose length is 10. The operation time of total algorithm (the following time including image exports time) and two performance index (MSE and PSNR) were evaluated in the experiment (The noise added to original image is Gauss white noise which obeys N(0,202) distribution). The material data is showed in table 1~2 and figure 2~7.
Fig.2 Original image Fig.3 Image with noise Fig.4 Hard thresholding
Fig.5 Soft thresholding Fig.6 Proportion-shrinking Fig.7 PCT algorithm Discussion: We can see from table 1 that when length of filters is short the effect of improved Mallat algorithm is not so good (this mainly because fast FFT algorithm is unadaptable to short signal), but with the increasing of length, under equal conditions, the improved Mallat algorithm will save much more time than computing using Mallat algorithm directly. Especially when the wavelet filters are relatively large in length, the closer to the length of the signal, the better the effect is. From table 2 we conclude the PCT denoising algorithm can improve PSNR more 1.7dB than soft thresholding method and more 0.5dB than proportion-shrinking method. At the same time from Fig.2- Fig.7 we also find that visual effect of PCT algorithm is better than that of others. Conclusion: In sum, in multi-resolution analysis, the length difference between the low-pass filter and high-pass filter should be as little as possible, so that the method above can be used effectually. The algorithm proposed in this paper has good parallel capability, is very suitable for the parallel computation of the digital signal processor (DSP). In addition, if the signal sequence contains big prime factor arithmetic Fourier transform (AFT) can be used.[4] Seeing that the fast development of VLSI technology and maturity of DSP technology based on FFT, the hardware system of wavelet denoiding fast algorithm is easy to come true, so it will become more efficient in the future. Reference 1. Peiqing Chen. Digital signal processing [M]. Beijing: Tsinghua University Press, 1998: 215-252. 2. Fusheng Yang. The engineering analysis and application of wavelet transform. Beijing: Science Press, 1999: 1- 176. 3. Mallat S. A theory for multiresolution signal decomposition: the wavelet representation [J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 1989, 11(4): 674-693. 4. M. K. Mihçak, I. Kozintsev, K. Ramchandran, P. Moulin. "Low-complexity image denoising based on statistical modeling of wavelet coefficients". IEEE Trans. on Signal Processing, December 1999, 6(12): 300-303 5. Ching P C, So H C, Wu S Q, "On wavelet denoising and its application to time delay estimation", [J]. IEEE Trans. Signal Processing, 1993, 47(10): 2879-2882 |
| © NDT.net |