Compile the MATLAB program for wavelet signal decomposition and reconstruction

P.s. () thanks to the netizen 'Li Ming yangyan ', he pointed out a major bug in mydwt and myidwt, see the one-dimensional signal wavelet decomposition and reconstruction program published today. P.s.: (last November) a series of articles on Wavelet Transform and image processing were published in, sharing their experiences and programming in the process of learning wavelet. Over the past six months, I have been grateful for your attention and help. During discussions and exchanges, I have continuously expanded my knowledge from the questions raised by everyone, he has a deep understanding of wavelet theory and its application. Modify the program in the blog Based on the problems found during discussion and discussion. The procedures shared in the two articles on wavelet image decomposition and reconstruction have the following problems: (1) the wavelet functions used by the program are only non-standard Haar wavelet, its filter group is lo_d = [1/2 1/2], hi_d = [-1/2 1/2], which is fixed in the mydwt2.m program and cannot be other wavelet functions. (2) in non-standard Haar wavelet, the details of high-frequency coefficients (contour, edge, and other features) are not obvious in the decomposed coefficient matrix. (3) in the mydwt2 function, the matrix object of column transformation is the input matrix. This is incorrect. The matrix object should be the cache matrix after row transformation. (4) the output of mydwt2 function uses [ll, HL, LH, HH], which is not very standard, should be changed to [Ca, CV, CH, CD], that is, the first-level wavelet transform output coefficient matrix has four parts: average, vertical, horizontal, and diagonal. (5) the output y of the mywavedec2 function is a matrix of the same size as the input matrix X, and the average and detail coefficients after N-level decomposition have been combined into one. In fact, this definition is only valid for Haar wavelet. (6) In the original program, the modmat function should be called to trim the image matrix so that it can be divisible by the nth power of 2, mainly to generate a tower structure image, after the above problem is corrected, this modmat function is no longer needed. In response to the above problems, I have corrected the program and published it in three articles today. Please click to view it. The newly modified program is more concise and easy to understand, and its functions are also enhanced. It can use any wavelet function for wavelet decomposition, based on the wavelet decomposition coefficient matrix, the low-frequency coefficients and original images at the specified decomposition level can be reconstructed. 1. Problems and Solutions in the Process of wavelet image decomposition and reconstruction
MATLAB program-V2.0
MATLAB program-V2.0
Http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2514119.aspx () MATLAB wavelet analysis toolbox rich functions and powerful simulation function for us to learn wavelet, use the wavelet provides a convenient and fast way, however, if we need to have a deep understanding of the principles of wavelet analysis and learn well and make good use of wavelet, we should try to use our own programs to implement wavelet transformation and signal analysis, try to call less functions provided by MATLAB in your own program and write relevant wavelet functions with your own understanding. This process is a process of exploration and knowledge, it makes us feel the power of wavelet and the pleasure of learning. Next, I will share the one-dimensional and two-dimensional wavelet signal decomposition and reconstruction MATLAB programs. I also hope that some friends can share the self-developed programs and learn together to improve the program efficiency and simplicity. The first thing to note is that although you compile your own MATLAB program, it does not mean that you do not need the built-in functions of MATLAB at all. What we need to write is the main function to implement wavelet transformation, while the basic functions such as plotting still need to use MATLAB functions. In addition, according to the filter group principle of wavelet transform, the original signal must be processed by low-pass and high-pass filters. convolution is involved here. Convolution-the implementation of the FFT algorithm, I believe many friends can use MATLAB, C language, and so on, but compared with the FFT program compiled in the machine language that comes with MATLAB, the computing speed is generally several times slower and dozens of times slower. Therefore, the Conv () function of MATLAB is called directly when convolution is involved in my program. We know that the first-level Decomposition Process of wavelet transform is that the original signal performs low-pass and high-pass filtering respectively, and then performs binary subsample separately to obtain low-frequency and high-frequency (also known as average and detail) the two-part coefficient, while the multi-level decomposition is a recursive process that implements wavelet decomposition of the low-frequency coefficient obtained by the upper-level decomposition. The procedures for one-dimensional wavelet decomposition are as follows:

Function [Ca, CD] = mydwt (x, lpd, HPD, dim );
% Function [Ca, CD] = mydwt (x, lpd, HPD, dim) for one-dimensional discrete wavelet decomposition of input sequence X, output decomposition sequence [Ca, CD]
% Input parameter: X -- input sequence;
% LPD -- low-pass filter;
% HPD -- Qualcomm filter;
% Dim -- Wavelet Decomposition series.
% Output parameter: Ca-wavelet factorization coefficient of the average part;
% Cd -- Wavelet factorization coefficient of detail.

CA = x; % initialize ca, CD
Cd = [];
For I = 1: dim
Cvl = Conv (Ca, lpd); % low-pass filtering. To improve the running speed, call the Conv function Conv () provided by Matlab ()
DNL = downspl (cvl); % calculate the decomposition coefficient of the average part by subsampling
Cvh = Conv (Ca, HPD); % Qualcomm Filter
Dnh = downspl (cvh); % The detailed coefficient after decomposition of this layer is obtained through subsample
CA = DNL; % The average Partial Coefficient after sampling enters the next layer of Decomposition
Cd = [cd, dnh]; % store the detailed coefficients obtained from the layer into the sequence CD
End

Function Y = downspl (X );
The % Function Y = dowmspl (x) samples the input sequence and outputs the sequence y.
% Sampling is to take its even bits for the input sequence and discard odd bits. For example, if X = [x1, x2, X3, X4, X5], then Y = [X2, X4].

N = length (x); % read input sequence length
M = floor (n/2); % the length of the output sequence is half of the length of the input sequence (take the integer part when decimal)
I = 1: m;
Y (I) = x (2 * I );

Reconstruction is the inverse process of decomposition. The low-frequency and high-frequency coefficients are sampled on top, low-pass, and high-pass filters respectively. Note that the number of low frequency and high frequency coefficients at the same level must be equal during reconstruction.

Function Y = myidwt (Ca, CD, LPR, HPR );
% Function myidwt () performs inverse Discrete Wavelet Transformation on the input wavelet decomposition coefficient to reconstruct the signal sequence y
% Input parameter: ca -- Wavelet factorization coefficient of the mean part;
% Cd -- Wavelet factorization coefficient of details;
% LPR, HPR-rebuild the low-pass and high-pass filters used.

LCA = length (CA); % find the length of the decomposition coefficient of the average and detail parts
LCD = length (CD );

While (LCD)> = (LCA) % in each layer reconstruction, the length of Ca and CD must be equal, so after each layer reconstruction,
% If the LCD is smaller than the LCA value, the reconstruction is stopped. At this time, the Ca is the reconstruction signal sequence y.
Upl = upspl (CA); % samples the average Partial Coefficient
Cvl = Conv (UPL, LPR); % low-pass convolution

Cd_up = Cd (LCD-LCA + 1: LCD); % obtain the coefficient of detail required for the current layer reconstruction. The length is equal to the length of the average coefficient of the current layer.
UPH = upspl (cd_up); % samples the coefficient of detail
Cvh = Conv (UPH, HPR); % Qualcomm convolution

CA = cvl + cvh; % use the sequence of current layer reconstruction to update CA for next layer Reconstruction
Cd = Cd (1: LCD-LCA); % discard the detailed coefficient used for this layer reconstruction, and update the CD
LCA = length (CA); % obtain the average length of the coefficient and the coefficient of the detail part used for the next layer reconstruction.
LCD = length (CD );
End % LCD <LCA, reconstruction complete, end cycle
Y = Ca; % the output reconstruction sequence y is equal to the average Partial Coefficient sequence CA after reconstruction.

Function Y = upspl (X );
The % Function Y = upspl (X) is used to sample the input one-dimensional sequence X, that is, to sample each element in sequence X.
% Insert zero, for example, x = [x1, x2, X3, X4]. After sampling, y = [X1, 0, X2, 0, X3, 0, X4].

N = length (x); % read input sequence length
M = 2 * N-1; % the length of the output sequence is twice the length of the input sequence and then minus one
For I = 1: M % The even bits of the output sequence are 0, and odd bits are input sequence elements whose order is equal to the corresponding position.
If Mod (I, 2)
Y (I) = x (I + 1)/2 );
Else
Y (I) = 0;
End
End

We know that two-dimensional wavelet decomposition and reconstruction can be implemented using a series of one-dimensional wavelet decomposition and reconstruction. The following procedures describe and reconstruct two-dimensional wavelet based on Haar wavelet:

Function [ll, HL, LH, HH] = mydwt2 (X );
% Function mydwt2 () implements two-dimensional wavelet decomposition on input R * C dimension matrix X, and outputs four decomposition coefficient submatrices [ll, HL, LH, HH]
% Input parameter: X -- input matrix, which is the R * C dimension matrix.
% Output parameter: LL, HL, LH, hh -- is the four equal child matrices of the decomposition coefficient matrix. The size is R/2 * C/2 dimension.
% LL: Low-Frequency partial decomposition coefficient; hl: vertical decomposition coefficient;
% LH: horizontal decomposition coefficient; hh: diagonal decomposition coefficient.

LPD = [1/2 1/2]; HPD = [-1/2 1/2]; % default low-pass and high-pass filters
[Row, Col] = size (x); % read the size of the input matrix

For j = 1: row % first, one-dimensional discrete wavelet decomposition is performed on each row of the input matrix.
Tmp1 = x (J ,:);
[A1, CD1] = mydwt (tmp1, lpd, HPD, 1 );
X (J, :) = [a1, CD1]; % store the decomposition coefficient sequence in matrix X to obtain [L | H]
End
For k = 1: Col % then perform one-dimensional discrete wavelet decomposition on each column of the input matrix.
Tmp2 = x (:, k );
[Ca 2, CD2] = mydwt (tmp2, lpd, HPD, 1 );
X (:, K) = [Ca, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]; %
End

LL = x (1: Row/: COL/2); % ll is the upper left corner of matrix X
LH = x (Row/2 + 1: Row, 1: COL/2); % LH is the lower left corner of matrix X
Hl = x (1: Row/2, COL/2 + 1: COL); % HL is the upper right corner of matrix X
HH = x (Row/2 + 1: Row, COL/2 + 1: COL); % hh is the lower right corner of matrix X

Function Y = myidwt2 (LL, HL, LH, HH );
% Function myidwt2 () performs inverse wavelet transformation on the input sub-matrix sequence to reconstruct the matrix Y
% Input parameter: LL, HL, LH, hh -- four matrices whose sizes are R * C
% Output parameter: y -- is a matrix with a size of 2R * 2C.

LPR = [1]; HPR = [1-1]; % default low-pass and high-pass filters
Tmp_mat = [ll, HL; LH, HH]; % combines the input four Matrices Into a Matrix
[Row, Col] = size (tmp_mat); % to obtain the number of rows and columns in the composite matrix

For k = 1: Col % split each column of the tmp_mat composite matrix into upper and lower half
A1 = tmp_mat (1: Row/2, K); % the two parts separated are used as the average coefficient sequence cache and the detail coefficient sequence CD1 respectively.
CD1 = tmp_mat (Row/2 + 1: Row, k );
Tmp1 = myidwt (CA1, CD1, LPR, HPR); % reconstruction Sequence
Yt (:, K) = tmp1; % stores the reconstruction sequence in the corresponding column of the yt matrix to be output. In this case, y = [L | H]
End

For j = 1: row % splits each row of output matrix Y into two halves.
Ca 2 = yt (J, 1: COL/2); % the two separate parts are used as the average coefficient sequence.
CD2 = yt (J, COL/2 + 1: COL );
Tmp2 = myidwt (A2, CD2, LPR, HPR); % reconstruction Sequence
Yt (J, :) = tmp2; % stores the reconstruction sequence in the corresponding row of the yt matrix to be output, and obtains the final output matrix Y = yt.
End
Y = yt;

I try to implement modularization when writing programs to facilitate calls between different functional programs. The above program implements the functions of one-dimensional and two-dimensional wavelet decomposition and reconstruction. Next we can explore technologies such as wavelet signal compression and image fusion.