数字图像处理MATLAB学习笔记（六）

数字图像处理MATLAB学习笔记（六）

Wavelets

1 Background

这里我们简单回顾一下傅里叶变换，并对一些概念做一下定义：

设正向DFT为$T(u,v,\dots )$，根据一般关系DFT和IDFT关系式可以表示为：
$$\begin{gathered} T(u,v,\dots )=\sum _{x,y}f(x,y)g_{u,v,\dots }(x,y) \\ f(x,y)=\sum _{x,y}T(u,v,\dots )h_{u,v,\dots }(x,y) \end{gathered}$$
其中$g_{u,v,\dots }$和$h_{u,v,\dots }$分别表示为正变换核(Lowpass)和反变换核(Highpass)。我们用级数展开可以表达关系为：
$$h_{u,v}(x,y)=g_{u,v}^{\ast }(x,y)=\frac{1}{\sqrt{MN}}\exp (j2\pi (ux/M+vy/N))$$
其中，$j=\sqrt{-1}$，$\ast$表示为复共轭算子，$u=0,1,2,\dots ,M-1$且$v=0,1,2,\dots ,N-1$。u和v分别表示水平和垂直频率，变换核可分为
$$\begin{gathered} h_u(x,y)=\frac{1}{\sqrt{M}}\exp (j2\pi x/M) \\ {h}_v(x,y)=\frac{1}{\sqrt{N}}\exp (j2\pi y/N) \\ \text{且，二者是正交的，正交关系为：} \\ ⟨h_r,h_s⟩={\delta }_{rs}=\{\begin{array}{ll}1& r=s\\ 0& \text{otherwise}\end{array} \end{gathered}$$
DWT(Discrete Wavelets Transform)离散小波变换指的是基本小波的尺度和平移进行离散化。

针对小波核，具有的一般特性为：

• 可分离性、可伸缩性、可平移性：

变换核可用三个可分的二维小波来表示：
$$\begin{gathered} {\psi }^H(x,y)=\psi (x)\phi (y) \\ {\psi }^V(x,y)=\phi (x)\psi (y) \\ {\psi }^D(x,y)=\psi (x)\psi (y) \end{gathered}$$
其中，${\psi }^H(x,y)$、${\psi }^V(x,y)$和${\psi }^D(x,y)$分别为水平Horizontal、垂直Vertical、对角Diagonal小波。且是二维可分的尺度函数为$\phi (x,y)=\phi (x)\phi (y)$。

每个二维函数都是两个一维是函数的、平方可积的尺度和小波函数的乘积：${\phi }_{j,k}=2^{j/2}\phi (2^{j}x-k)$和${\psi }_{j,k}=2^{j/2}\psi (2^{j}x-k)$。平移参数$k$决定这些一维函数沿x轴的位置，尺度j决定他们的宽度并且按照$2^{j/2}$控制高度或振幅。

其中，与展开函数相关联的是母小波的二进制尺度和证书平移$\psi (x)={\psi }_{0,0}(x)$，并且尺度函数$\phi (x)={\phi }_{0,0}(x)$。

• Multiresolution Compatibility多分辨率的兼容性：一维尺度函数满足下要求即可用于多分辨率分析：
  • ${\phi }_{j,k}$对整数平移是正交的。
  • 在低尺度或低分辨率下，以一系列${\phi }_{j,k}$的展开函数来描述的一组函数包含在可以用高尺度描述的函数中。
  • 唯一可在每个尺度上描述的函数是$f(x)=0$。
  • 在$j\to \infty$，任何函数都可以以任意精度来描述。
• Orthogonality正交性

小波变换和Fourier Transform之间的关系：

原始信号中，时域的分辨率高，频域分辨率为零，意味着我们可以在时域中分辨出非常小的特征，而在频域中没有特征。

傅里叶变换在频域具有高分辨率，在时域具有零分辨率。

短时傅里叶变换相对于上面的两种方法，它在时域和频域具有中等的分辨率。

小波变换有如下的特点：

• 对于小频率值，频域分辨率高，时域分辨率低

• 对于大频率值，频域分辨率低，时域分辨率高

小波变换在频域分辨率和时域分辨率两者之间做了权衡：在与时间相关的特征上上，它在时域具有高分辨率；而在与频率相关的特征上，它在频域具有高分辨率。

2 The Fast Wavelets Transform(FWT)

根据小波核，我们有个重要推论：$\phi (x)$和$\psi (x)$可以用双分辨率副本的线性组合来表达。

通过级数展开，我们得到：
$$\begin{gathered} \phi (x)=\sum _n{h}_{\phi }(n)\sqrt{2}\phi (2x-n) \\ \psi (x)=\sum _n{h}_{\psi }(n)\sqrt{2}\phi (2x-n) \end{gathered}$$
其中，$h_{\phi }$和$h_{\psi }$分别是尺度sacling和小波矢量wavelet vectors。

关于FWT的迭代计算如图：

用公式表达如下：
$$\begin{gathered} W_{\psi }^D(j,m,n)=h_{\psi }(-m)\odot [h_{\psi }(-n)\odot W_{\phi }(j+1,m,n){|}_{n=2k,k\ge 0}]{|}_{m=2k,k\ge 0} \\ {W}_{\psi }^V(j,m,n)=h_{\phi }(-m)\odot [h_{\psi }(-n)\odot W_{\phi }(j+1,m,n){|}_{n=2k,k\ge 0}]{|}_{m=2k,k\ge 0} \\ {W}_{\psi }^H(j,m,n)=h_{\psi }(-m)\odot [h_{\phi }(-n)\odot W_{\phi }(j+1,m,n){|}_{n=2k,k\ge 0}]{|}_{m=2k,k\ge 0} \\ {W}_{\phi }(j,m,n)=h_{\phi }(-m)\odot [h_{\phi }(-n)\odot W_{\phi }(j+1,m,n){|}_{n=2k,k\ge 0}]{|}_{m=2k,k\ge 0} \end{gathered}$$
其中，$h_{\phi }$和$h_{\psi }$分别是lowpass 和 highpass decomposition filters。

从迭代图中可以看出

• $W_{\phi }$是基于低通滤波器产生的，所以被称为approximation coefficients近似系数。
• $W_{\psi }^i\text{for }i=H,V,D$分别是水平、垂直、对角的detail coefficients。

由于FWT是向创建低分辨率分量，$f(x,y)$表示高分辨率分量，所以在第一层叠戴循环时，我们可以令$W_{\phi }(j+1,m,n)=f(x,y)$

2.1 FWTs Using the Wavelet Toobox

在MATLAB中，我们若想生成FWT的filter，可以使用函数wfilters：

[Lo_D, Hi_D, Lo_R, Hi_R] = wfilters(wname)

其中，Lo_D、Hi_D、Lo_R、Hi_R都是行向量，分别为低通分解、高通分解、低通重建、高通重建滤波器；wname指的是小波滤波系数，具体如下表所示：

Wavelet	wfamily	wname
Haar	‘haar’	‘haar’
Daubechies	‘db’	‘db2’、‘db3’、…‘db5’
Coiflets	‘coif’	‘coif1’、‘coif2’、…‘coif5’
Symlets	‘sym’	‘sym1’、‘sym2’、…‘sym5’
Discrete Meyer	‘dmey’	‘dmey’
Biorthogonal	‘bior’	‘bior1.1’、‘bior1.3’、‘bior1.5’、‘bior2.2’…
Reverse Biorthogonal	‘rbio’	‘rbior1.1’、‘rbior1.3’、‘rbior1.5’、‘rbior2.2’…

想查询更多小波族相关信息，可以在MATLAB中输入命令：

waveinfo(wfamily)

Frequently coupled filter pairs频率耦合滤波器对用如下语句代替：

[F1, F2] = wfilters(wname, type)

其中，type设置为’d’、‘r’、‘l’、'h’分别得到分解、重建、低通、高通滤波器。分解或低通滤波器将返回到F1，结果放在F2中。

如果想获得归一化正交变换的尺度和小波函数的数字近似值，可以使用函数wavefun：

[phi, psi, xval] = wavefun(wname, iter)

其中，phi和psi为返回的近似向量，phi为尺度scale，psi为小波wavelet，xval为评估向量。iter比啊是迭代次数。对于双正交变换，命令则变为：

[phi1, psi1, phi2, psi2, xval] = wavefun(wname, iter)

其中，phi1和psi1是decomposition function分解函数，phi2和psi2是reconstruction function重建函数。

Example Haar filters, scaling, and wavelet functions

>> [Lo_D, Hi_D, Lo_R, Hi_R] = wfilters('haar')
Lo_D =
    0.7071    0.7071
Hi_D =
   -0.7071    0.7071
Lo_R =
    0.7071    0.7071
Hi_R =
    0.7071   -0.7071
>> [phi, psi, xval] = wavefun('haar', 10);
>> xaxis = zeros(size(xval));
>> subplot(121);plot(xval, phi, 'k', xval, xaxis, '--k');
>> axis([0 1 -1.5 1.5]); axis square
>> title('Haar Scaling Function');
>> subplot(122);plot(xval, psi, 'k', xval, xaxis, '--k');
>> axis([0 1 -1.5 1.5]); axis square
>> title('Haar Wavelet Function');

其中，Haar的尺度函数和小波函数都是不连续且紧密支撑的，支撑是1，支撑区间的优先区间之外是0。

有了wavefilter，计算wavelet transform小波变换一般使用函数wavedec2：

[C, S] = wavedec2(X, N, Lo_D, Hi_D)
[C, S] = wavedec2(X, N, wname)

其中，X是输入图像，N是被计算的尺度数，输出的C和S分别为行向量C和记录C中系数排列的矩阵S。

Example A simple FWT using Haar filters：

>> f = magic(4)
f =
    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1
>> [c1, s1] = wavedec2(f, 1, 'haar')
c1 =
  Columns 1 through 9
   17.0000   17.0000   17.0000   17.0000    1.0000   -1.0000   -1.0000    1.0000    4.0000
  Columns 10 through 16
   -4.0000   -4.0000    4.0000   10.0000    6.0000   -6.0000  -10.0000
s1 =
     2     2
     2     2
     4     4

其中，f是一个4x4的图像，被变换为1x16的小波分解向量c1和3x2的记录矩阵s1。

记录矩阵s1记录了矩阵的大小：每个向量的第一个元素是参考细节或近似矩阵的行数，第二个元素是列数。

• s1(end, :)记录了原始图像f的大小；
• s1(2, :)记录了3个细节系数矩阵的大小；
• s1(1, :)记录了计算近似矩阵的大小；

当我们的计算尺度为两个尺度时：

>> [c2, s2] = wavedec2(f, 2, 'haar') % 分解2次
c2 =
  Columns 1 through 9
   34.0000         0         0    0.0000    1.0000   -1.0000   -1.0000    1.0000    4.0000
  Columns 10 through 16
   -4.0000   -4.0000    4.0000   10.0000    6.0000   -6.0000  -10.0000
s2 =
     1     1
     1     1
     2     2
     4     4

这里，s2记录矩阵记录的数据为：

• s1(end, :)记录了原始图像f的大小；
• s1(3, :)记录了尺度/分解级别为2时，3个细节系数矩阵的大小；
• s1(2, :)记录了尺度/分解级别为1(第一次)时，3个细节系数矩阵的大小；
• s1(1, :)记录了最后近似矩阵的大小；

2.2 FWTs without the Wavelet Toolbox

function wavefilter：创建小波滤波器

function [varargout] = wavefilter(wname, type)
%WAVEFILTER Create wavelet decomposition and reconstruction filters .
%   [VARARGOUT] = WAVEFILTER (WNAME, TYPE) returns the decomposition
%   and/or reconstruction filters used in the computation of the
%   forward and inverse FWT(fast wavelet trans form).
% 
%   EXAMPLES : 
%       [ld, hd, lr, hr] = wavefilter('haar')   Get the low and highpass
%                                               decomposition (ld, hd) 
%                                               and reconstruction
%                                               (lr, hr) filters for
%                                               wavelet 'haar'.
%       [ld, hd] = wavefilter('haar', 'd')      Get decomposition filters
%                                               ld and hd.
%       [lr, hr] = wavefilter('haar', 'r')      Get recons truction
%                                               filters lr and hr.
% 
%   INPUTS:
%       WNAME                   Wavelet Name
%       -----------------------------------------------------------------
%       'haar' or 'db1'         Haar
%       'db4'                   4th order Daubechies
%       'sym4'                  4th order Symlets
%       'bior6.8'               Cohen-Daubechies- Feauveau biorthogonal
%       'jpeg9.7'               Antonini-Barlaud-Mathieu-Daubechies
% 
%       TYPE                    Filter Type
%       -----------------------------------------------------------------
%       'd'                     Decomposition filters
%       'r'                     Reconstruction filters

% Check the input and output arguments.
error(nargchk(1, 2, nargin));

if (nargin == 1 && nargout ~= 4) || (nargin == 2 && nargout ~= 2)
    error('Invalid number of output arguments.');
end

if nargin == 1 && ~ischar(wname)
    error('WNAME must be a string.');
end

if nargin == 2 && ~ischar(type)
    error('TYPE must be a string.');
end

% Create filters for the requested wavelet.
switch lower(wname)
    case {'haar', 'db1'}
        ld = [1 1]/sqrt(2);
        hd = [-1 1]/sqrt(2);
        lr = ld;
        hr = -hd;
    case 'db4'
        ld = [-1.059740178499728e-002 3.288301166698295e-002 ...
            3.084138183598697e-002 -1.870348117188811e-001 ...
            -2.798376941698385e-002 6.308807679295904e-001 ...
            7.148465705525415e-001 2.303778133088552e-001];
        t = (0:7);
        hd = ld;
        hd(end:-1:1) = cos(pi*t).*ld;
        lr = ld;
        lr(end:-1:1) = ld;
        hr = cos(pi*t).*ld;
    case 'sym4'
        ld = [-7.576571478927333e-002 -2.963552764599851e-002 ...
            4.976186676320155e-001 8.037387518059161e-001 ...
            2.978577956052774e-001 -9.921954357684722e-002 ...
            -1.260396726203783e-002 3.222310060404270e-002];
        t = (0:7);
        hd = ld;
        hd(end:-1:1) = cos(pi*t).*ld;
        lr = ld;
        lr(end:-1:1) = ld;
        hr = cos(pi*t).*ld;
    case 'bior6.8'
        ld = [0 1.908831736481291e-003 -1.914286129088767e-003 ...
            -1.699063986760234-002 1.193456527972926e-002 ...
            4.973290349094079e-002 -7.726317316720414e-002 ...
            -9.405920349573646-002 4.207962846098268e-001 ...
            8.259229974584023e-001 4.207962846098268e-001 ...
            -9.405920349573646e-002 -7.726317316720414e-002 ...
            4.973290349094079e-002 1.193456527972926e-002 ...
            -1.699063986760234e-002 -1.914286129088767e-003 ...
            1.908831736481291e-003];
        hd = [0 0 0 1.442628250562444e-002 -1.446750489679015e-002 ...
            -7.872200106262882e-002 4.036797903033992e-002 ...
            4.178491091502746e-001 -7.589077294536542e-001 ...
            4.178491091502746e-001 4.036797903033992e-002 ...
            -7.872200106262882e-002 -1.446750489679015e-002 ...
            1.442628250562444e-002 0 0 0 0];
        t = (0:17);
        lr = cos(pi*(t+1)).*hd;
        hr = cos(pi*t).*ld;
    case 'jpeg9.7'
        ld = [0 0.02674875741080976 -0.01686411844287495 ...
            -0.07822326652898785 0.2668641184428723 ...
            0.6029490182363579 0.2668641184428723 ...
            -0.07822326652898785 -0.01686411844287495 ...
            0.02674875741080976];
        hd = [0 -0.09127176311424948 0.05754352622849957 ...
            0.5912717631142470 -1.115087052456994 ...
            0.5912717631142470 0.05754352622849957 ...
            -0.09127176311424948 0 0] ;
        t = (0:9);
        lr = cos(pi*(t+1)).*hd;
        hr = cos(pi*t).*ld;
    otherwise
        error('Unrecognizable wavelet name (WNAME).');
end

% Output the requested filters.
if (nargin == 1)
    varargout(1:4) = {ld, hd, lr, hr};
else
    switch lower(type(1))
        case 'd'
            varargout = {ld, hd};
        case 'r'
            varargout = {lr, hr};
        otherwise
            error('Unrecognizable filter TYPE');
    end
end

function wavefast：实现FWT

function [c, s] = wavefast(x, n, varargin)
%WAVEFAST Computes the FWT of a '3-D extended' 2-D array.
%   [C, L] = WAVEFAST(X, N, LP, HP) computes ' PAGES' 2D N-level
%   FWTs of a 'ROWS x COLUMNS x PAGES' matrix X with respect to
%   decomposition filters LP and HP .
% 
%   [C, L] = WAVEFAST(X, N, WNAME) performs the same operation but
%   fetches filters LP and HP for wavelet WNAME using WAVEFILTER.
%   Scale parameter N must be less than or equal to log2 of the
%   maximum image dimension. Filters LP and HP must be even. To
%   reduce border distortion, X is symmetrically extended. That is,
%   if X = [c1 c2 c3 ... cn] (in 1D), then its symmetric extension
%   would be [... c3 c2 c1 c1 c2 c3 ... cn cn cn-1 cn-2 ...].
% 
%   OUTPUTS:
%       Vector C is a coef f icient decomposition vector:
% 
%           C = [ a1(n) ... ak(n) h1(n) ... hk(n) v1(n) ... vk(n)
%                 d1 (n) ... dk(n) h1 (n-1) ... d1(1) ... dk(1) ]
% 
%           where ai, hi, vi, and di for i = 0,1, ... k are columnwise
%           vectors containing approximation, horizontal, vertical, and 
%           diagonal coefficient matrices, respectively, and k is the
%           number of pages in the 3-D extended array X. C has 3n + 1
%           sections where n is the number of wavelet decomposi tions .
% 
%           Matrix S is an [(n+2) x 2] bookkeeping matrix if k = 1;
%           else it is [(n+2) x 3]:
% 
%               S =[ sa(n, :); sd(n, :); sd(n-1, :); ... ; sd(1, :); sx ]
% 
%           where sa and sd are approximation and detail size entries.

% Check the input arguments for reasonableness
error(nargchk(3, 4, nargin));
if nargin == 3
    if ischar(varargin{1})
        [lp, hp] = wavefilter(varargin{1}, 'd');
    else
        error('Missing wavelet name');
    end
else
    lp = varargin{1};
    hp = varargin{2};
end

% Get the filter length, 'lp', input array size, 'sx', and number of pages,
% 'pages', in extended 2-D array x
fl = length(lp);
sx = size(x);
pages = size(x, 3);

if ((ndims(x) ~= 2) && (ndims(x) ~= 3)) || (min(sx) < 2) || ~isreal(x) ...
        || ~isnumeric(x)
    error('X must be a real, numeric 2-D or 3-D matrix.');
end

if (ndims(lp) ~= 2) || ~isreal(lp) || ~isnumeric(lp) ...
        || (ndims(hp) ~= 2) || ~isreal(hp) || ~isnumeric(hp) ...
        || (fl ~= length(hp)) || rem(fl, 2) ~= 0
    error(['LP and HP must be even and equal length real, ' ...
    'numeric filter vectors.']);
end

if ~isreal(n) || ~isnumeric(n) || (n < 1) || (n > log2(max(sx)))
    error('N must be a real scaler between 1 and log2(max(size(X)))');
end

% Init the starting output data structures and initial approximation
c = [];
s = sx(1:2);
app = cell(pages, 1);
for i = 1:pages
    app{i} = double(x(:, :, i));
end

% For each decomposition ...
for i = 1:n
    % Extend the appriximation symmetrically.
    [app, keep] = symextend(app, fl, pages);
    
    % Convolve rows with HP and downsample. Then convolve columns with HP
    % and LP to get the diagonal and vertical coefficients.
    rows = symconv(app, hp, 'row', fl, keep, pages);
    coefs = symconv(rows, hp, 'col', fl, keep, pages);
    c = addcoefs(c, coefs, pages);
    s = [size(coefs{1}); s];
    coefs = symconv(rows, lp, 'col', fl, keep, pages);
    c = addcoefs(c, coefs, pages);
    % Convolve rows with LP and downsample. Then convolve columns with HP
    % and LP to get the horizontal and next approximation coeffcients.
    rows = symconv(app, lp, 'row', fl, keep, pages);
    coefs = symconv(rows, hp, 'col', fl, keep, pages);
    c = addcoefs(c, coefs, pages);
    app = symconv(rows, lp, 'col', fl, keep, pages);
end

% Append the final appriximation structures
c = addcoefs(c, app, pages);
s = [size(app{1}); s];
if ndims(x) > 2
    s(:, 3) = size(x, 3);
end

%----------------------------------------------------------------------
function nc = addcoefs(c, x, pages)
% Add 'pages' array coefficients to the wavelet decomposition vector.

nc = c;
for i = pages:-1:1
    nc = [x{i}(:)' nc];
end

%----------------------------------------------------------------------
function [y, keep] = symextend(x, fl, pages)
% Compute the number of coefficients to keep after convolution and
% downsampling. Then extend the 'pages' arrays of x in both dimensions.

y = cell(pages, 1);
for i = 1:pages
    keep = floor((fl + size(x{i}) - 1) / 2);
    y{i} = padarray(x{i}, [(fl - 1) (fl - 1)], 'symmetric', 'both');
end

%----------------------------------------------------------------------
function y = symconv(x, h, type, fl, keep, pages)
% For the 'pages' 2-D arrays in x convolve the rows or columns with h,
% downsample, and extract the center section since symmetrically extened.

y = cell(pages, 1);
for i = 1:pages
    if strcmp(type, 'row')
        y{i} = conv2(x{i}, h);
        y{i} = y{i}(:, 1:2:end);
        y{i} = y{i}(:, fl / 2 + 1:fl / 2 + keep(2));
    else
        y{i} = conv2(x{i}, h');
        y{i} = y{i}(1:2:end, :);
        y{i} = y{i}(fl / 2 + 1:fl / 2 + keep(1), :);
    end
end

function fwtcompare：自定义的FWT于wavedec2性能对比

function [ratio, maxdiff] = fwtcompare(f, n, wname)
%FWTCOMPARE Compare wavedec2 and wavefast.
%   [RATIO, MAXDIFF] = FWTCOMPARE(E, N, WNAME) compares the operation of
%   Wavelet Toolbox function WAVEDEC2 and custom function WAVEFAST.
% 
%   INPUTS:
%               F       Image to be trans formed.
%               N       Number of scales to compute.
%               WNAME   Wavelet to use.
% 
%   OUTPUTS:
%               RATIO Execution time ratio (custom/toolbox)
%               MAXDIFF Maximum coefficient di f ference .

% Get transform and computation time for wavedec2.
w1 = @() wavedec2(f, n, wname);
reftime = timeit(w1);

% Get transform and computation time for wavefast.
w2 = @() wavefast(f, n, wname);
t2 = timeit(w2);

% Compare the results.
ratio = t2 / reftime;
maxdiff = abs(max(w1() - w2()));

对比结果如下：

>> [ratio, maxdifference] = fwtcompare(f, 5, 'db4')
>> [ratio, maxdifference] = fwtcompare(f, 5, 'db4')
ratio =
    0.8640
maxdifference =
   1.6467e-12

3 Working with Wavelet Decomposition Sructures

我们知道自定义和Toolbox中的小波计算函数都会返回**[c,S]形式的结构变量。其中c**是transform coefficient vector变换系数向量，S是bookkeeping matrix记录矩阵。

这里的c我们是将多尺度二维小波变换系数集成为单一的一维向量：
$$c=[A_N(:{)}^{\prime }{H}_N(:{)}^{\prime }{V}_N(:{)}^{\prime }{D}_N(:{)}^{\prime }\dots H_i(:{)}^{\prime }{V}_i(:{)}^{\prime }{D}_i(:{)}^{\prime }\dots H_1(:{)}^{\prime }{V}_1(:{)}^{\prime }{D}_1(:{)}^{\prime }]$$
其中，$A_N$是第N个分解等级的近似系数矩阵，$H_i$、$V_i$、$D_i$分别是级别i的水平、垂直、对角的变换系数矩阵，i=1,2,…,N。

设定一共进行N级分解，也就是说，迭代循环N次，c中应包含3N+1个元素。当i=1时，计算最高尺度系数；当i=N时，计算最低尺度系数。所以，c的系数是由低尺度到高尺度排列的。

矩阵S时(N+2)x2的记录数组，形式为：
$$S=[s{a}_N;s{d}_N;s{d}_{N-1};\dots s{d}_i;\dots s{d}_1;sf]$$
$s{a}_N$表示第N级的近似矩阵$A_N$，$s{d}_i$表示第i级的细节矩阵大小（这里三个方向的细节矩阵都被转换为列向量进行存储），$sf$表示原始图像的大小。

对于RGB彩色图像，即三位数组作为输入时，我们将其作为extended 2-D array扩展的二维数组进行处理。也就是说，矩阵由矩阵构成，(实际上很简单)，表达如下：
$$\begin{gathered} A_N(:{)}^{\prime }=[A_N^1(:{)}^{\prime }{A}_N^2(:{)}^{\prime }\dots A_N^k(:{)}^{\prime }] \\ {H}_i(:{)}^{\prime }=[H_i^1(:{)}^{\prime }{H}_i^2(:{)}^{\prime }\dots H_i^k(:{)}^{\prime }] \\ {V}_i(:{)}^{\prime }=[V_i^1(:{)}^{\prime }{V}_i^2(:{)}^{\prime }\dots V_i^k(:{)}^{\prime }] \\ {D}_i(:{)}^{\prime }=[D_i^1(:{)}^{\prime }{D}_i^2(:{)}^{\prime }\dots D_i^k(:{)}^{\prime }] \end{gathered}$$
其中，k时扩展数组的页数，i时分解的级数。这样我们得到的**[c,S]**，c还是有3N+1个元素构成，S则变成了(N+2)x3个数组。在数组中，第三列指定了c中二位数组的数目。

Example Wavelet Toolbox functions for manipulating transform decomposition vector c：

>> f = magic(8);
>> [c1, s1] = wavedec2(f, 3, 'haar');
>> size(c1)
ans =
     1    64
>> s1
s1 =
     1     1
     1     1
     2     2
     4     4
     8     8
>> approx = appcoef2(c1, s1, 'haar')
approx =
  260.0000
>> horizdet2 = detcoef2('h', c1, s1, 2)
horizdet2 =
   1.0e-13 *
         0   -0.2842
         0         0
>> newc1 = wthcoef2('h', c1, s1, 2);
>> newhorizdet2 = detcoef2('h', newc1, s1, 2)
newhorizdet2 =
     0     0
     0     0

这里，function appcoef2是根据小波分解结构[C，S]，计算N级重构系数矩阵，即获得近似矩阵：

a = appcoef2(c, s, wname)

function detcoef2是获取detail matrix，语法如下：

d = detcoef2(o, c, s, n)

其中，o可以设置为’h’、’ v’、'd’分别表示水平、垂直、对角的细节矩阵。

function wthcoef2是用于对二维信号的小波系数阈值进行处理的函数：

nc = wthcoef2(type, c, s, n, t, sorh)

type针对阈值近似系数设置为’a’，细节阈值可以传参’h’、’ v’、'd’分别表示水平、垂直、对角。n时被阈值处理后的分解级别的向量，阈值以向量t中相应的阈值为基础，如果忽略t，所有满足type和n规范的系数都将被赋值为0。

3.1 Editing Wavelet Decomposition Coefficients without the Wavelet Toolbox

function wavework：

function [varargout] = wavework(opcode, type, c, s, n, x)
%WAVEWORK is used to edit wavelet decomposition structures .
%   [VARARGOUT] = WAVEWORK (OPCODE, TYPE, c, s, N, X) gets the
%   coefficients specified by TYPE and N for access or modi fication
%   based on OPCODE.
% 
%   INPUTS:
%       OPCODE      Operation to perform
%       -----------------------------------------------------------------
%       'copy'      [varargout] = Y = requested (via TYPE and N)
%                   coefficient matrix
%       'cut'       [varargout] = [NC, Y] = New decompos ition vector
%                   (with requested coefficient matrix zeroed) AND
%                   requested coefficient matrix
%       'paste'     [varargout] = [NC] = new decomposition vector with
%                   coefficient matrix replaced by X
% 
%       TYPE        Coefficient category
%       -----------------------------------------------------------------
%       'a'         Approximation coef ficients
%       'h'         Horizontal details
%       'v'         Vertical details
%       'd'         Diagonal details
% 
%       [c, S] is a wavelet toolbox decomposition structure.
%       N is a decomposition level (Ignored if TYPE = 'a').
%       X is a 2- or 3-D coefficient matrix for pasting.

error(nargchk(4, 6, nargin));

if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
    error('C must be a row vector.');
end

if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
        ((size(s, 2) ~= 2) && (size(s, 3) ~= 3))
    error('S must be a real, numeric two- or three-column array.');
end

elements = prod(s, 2);

if (length(c) < elements(end)) || ...
        ~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
    error('[C S] must form a standard wavelet decomposition structre.');
end

if strcmpi(opcode(1:3), 'pas') && nargin < 6
    error('Not enough input arguments');
end

if nargin < 5
    n = 1;
end
nmax = size(s, 1) - 2;

aflag = (lower(type(1)) == 'a');
if ~aflag && (n > nmax)
    error('N exceeds the decompositions in [C S].');
end

switch lower(type(1))   % Make pointers into C
    case 'a'
        nindex = 1;
        start = 1;
        stop = elements(1);
        ntst = nmax;
    case {'h', 'v' ,'d'}
        switch type
            case 'h'
                offset = 0;     % Offset to details
            case 'v'
                offset = 1;
            case 'd'
                offset = 2;
        end
        nindex = size(s, 1) - n;    % Index to detail infow
        start = elements(1) + 3 * sum(elements(2:nmax - n + 1)) + ...
            offset * elements(nindex) + 1;
        stop = start + elements(nindex) - 1;
        ntst = n;
    otherwise
        error('TYPE must begin with "a", "h", "v" or "d".');
end

switch lower(opcode)    % Do requested action
    case {'copy', 'cut'}
        y = c(start:stop);
        nc = c;
        y = reshape(y, s(nindex, :));
        if strcmpi(opcode(1:3), 'cut')
            nc(start:stop) = 0;
            varargout = {nc, y};
        else
            varargout = {y};
        end
    case 'paste'
        if numel(x) ~= elements(end - ntst)
            error('X is not sized for the requested paste.');
        else
            nc = c;
            nc(start:stop) = x(:);
            varargout = {nc};
        end
    otherwise
        error('Unrecognized OPCODE.');
end

function wavecut

function [nc, y] = wavecut(type, c, s, n)
%WAVECUT Zeroes coefficients in a wavelet decomposition s tructure .
%   [NC, Y] = WAVECUT (TYPE, C, s, N) returns a new decomposition
%   vector whose detail or approximation coefficients  (based on TYPE
%   and N) have been zeroed. The coefficients that were zeroed are
%   re turned in Y.
% 
%   INPUTS:
%       TYPE        Coefficient category
%       -----------------------------------
%       'a'         Approximation coefficients
%       'h'         Horizontal details
%       'v'         Vertical details
%       'd'         Diagonal details 
% 
%       [c, S] is a wavelet data structure.
%       N specifies a decomposition level (ignored if TYPE = 'a') .

error(nargchk(3, 4, nargin));
if nargin == 4
    [nc, y] = wavework('cut', type, c, s, n);
else
    [nc, y] = wavework('cut', type, c, s);
end

function wavepaste

function nc = wavepaste(type, c, s, n, x)
%WAVEPASTE Puts coefficients in a wavelet decomposition structure .
%   NC = WAVEPASTE(TYPE, c, S, N, X) returns the new decomposition
%   structure after pasting X into it based on TYPE and N.
% 
%   INPUTS:
%       TYPE        Coefficient category
%       ----------------------------------------
%       'a'         Approximation coefficients
%       'h'         Horizontal details
%       'v'         Vertical details
%       'd'         Diagonal details
% 
%       [c, S] is a wavelet data structure.
%       N specifies a decomposition level (Ignored if TYPE = 'a').
%       X is a 2- or 3-D approximation or detail coefficient
%           matrix whose dimensions are appropriate for decomposition
%           level N.

error(nargchk(5, 5, nargin));
nc = wavework('paste', type, c, s, n, x);

function wavecopy

function y = wavecopy(type, c, s, n)
%WAVECOPY Fetches coefficients of a wavelet decomposition structure.
%   Y = WAVECOPY (TYPE, c, S, N) returns a coefficient array based on
%   TYPE and N.
% 
%   INPUTS:
%       TYPE        Coefficient category
%       ----------------------------------------
%       'a'         Approximation coefficients
%       'h'         Horizontal details
%       'v'         Vertical details
%       'd'         Diagonal details
% 
%       [c, S] is a wavelet data structure .
%       N specifies a decomposition level (ignored if TYPE = 'a') .

error(nargchk(3, 4, nargin));
if nargin == 4
    y = wavework('copy', type, c, s, n);
else
    y = wavework('copy', type, c, s);
end

3.2 Displaying Wavelet Decomposition Coefficients

function wavedisplay

function w = wavedisplay(c, s, scale, border)
%WAVEDISPLAY Display wavelet decomposition coefficients.
%   W = WAVEDISPLAY (c, S, SCALE, BORDER) displays and returns a
%   wavelet coef ficient image.
% 
%   EXAMPLES:
%       wavedisplay(c, s);                      Display w/defaults.
%       foo = wavedisplay(c, s);                Display and return.
%       foo = wavedisp1ay(c, s, 4);             Magnify the details.
%       foo = wavedisplay(c, s, -4);            Magnify absolute values.
%       foo = wavedisplay(c, s, 1, ' append');  Keep border values.
% 
%   INPUTS/ OUTPUTS:
%       [c, S] is a wavelet decomposition vector and bookkeeping
%       matrix.
% 
%       SCALE           Detail coefficient scaling
%       -------------------------------------------------------
%       0 or 1          Maximum range (default)
%       2,3...          Magnify default by the scale factor
%       -1, -2...       Magnify absolute values by abs (scale)
% 
%       BORDER          Border between wavelet decompositions 
%       -------------------------------------------------------
%       'absorb'        Border replaces image (default)
%       'append'        Border increases width of image
% 
%       Image W:        ------- ------- --------------- ------------------
%                       |      |       |              |
%                       | a(n) | h(n)  |              |
%                       |      |       |              |
%                       ------- -------     h(n-1)    |
%                       |      |       |              |
%                       | v(n) | d(n)  |              |       h(n-2)
%                       |      |       |              |
%                       ------- ------- ---------------
%                       |              |              |
%                       |   v(n-1)     |    d(n-1)    |
%                       |              |              |
%                       ------- ------- --------------- ------------------
%                       |                             | 
%                       |           v(n-2)            |       d(n-2)
%                       |                             |     
% 
%       Here, n denotes the decomposition step scale and a, h, v, d are
%       approximation, horizontal, vertical, and diagonal detail
%       coefficients, respectively.

% Check input arguments for reasonableness .
error(nargchk(2, 4, nargin));

if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
    error('C must be a row vector.');
end

if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
        ((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
    error('S must be a real, numeric, two- or three-column array.');
end

elements = prod(s, 2);

if (length(c) < elements(end)) || ...
        ~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
    error('[C S] must be a standard wavelet decomposition structure');
end

if (nargin > 2) && (~isreal(scale) || ~isnumeric(scale))
    error('SCALE must be a real, numeric scalar');
end

if (nargin > 3) && (~ischar(border))
    error('BORDER must be character string.');
end

if nargin == 2
    scale = 1;  % Default scale
end

if nargin < 4
    border = 'absorb';
end

% Scale coefficients and determine pad fill
absflag = scale < 0;
scale = abs(scale);
if scale == 0
    scale = 1;
end

[cd, w] = wavecut('a', c, s);
w = mat2gray(w);
cdx = max(abs(cd(:))) / scale;
if absflag
    cd = mat2gray(abs(cd), [0, cdx]);
    fill = 0;
else
    cd = mat2gray(cd, [-cdx, cdx]);
    fill = 0.5;
end

% Build gray image one decomposition at a time
for i = size(s, 1) - 2:-1:1
    ws = size(w);
    
    h = wavecopy('h', cd, s, i);
    pad = ws - size(h);
    frontporch = round(pad / 2);
    h = padarray(h, frontporch, fill, 'pre');
    h = padarray(h, pad - frontporch, fill, 'post');
    
    v = wavecopy('v', cd, s, i);
    pad = ws - size(v);
    frontporch = round(pad / 2);
    v = padarray(v, frontporch, fill, 'pre');
    v = padarray(v, pad - frontporch, fill, 'post');
    
    d = wavecopy('d', cd, s, i);
    pad = ws - size(d);
    frontporch = round(pad / 2);
    d = padarray(d, frontporch, fill, 'pre');
    d = padarray(d, pad - frontporch, fill, 'post');
    
    % Add 1 pixel white border and concatenate coefficients.
    switch lower(border)
        case 'append'
            w = padarray(w, [1 1], 1, 'post');
            h = padarray(h, [1 0], 1, 'post');
            v = padarray(v, [0 1], 1, 'post');
        case 'absorb'
            w(:, end, :) = 1;
            w(end, :, :) = 1;
            h(end, :, :) = 1;
            v(:, end, :) = 1;
        otherwise
            error('Unrecongnized BORDER parameter.');
    end
    w = [w h; v d];
end

% Display result. If the reconstruction is an extended 2-D array
% with 2 or more pages, display as a time sequence.
if nargout == 0
    if size(s, 2) == 2
        imshow(w);
    else
        implay(w);
    end
end

Example Transform coefficient display using wavedisplay：

>> f = imread('Fig0704.tif');
>> [c, s] = wavefast(f, 2, 'db4');
>> wavedisplay(c, s);
>> subplot(131);wavedisplay(c, s);
>> subplot(132);wavedisplay(c, s, 8);
>> subplot(133);wavedisplay(c, s, -8);

在不加比例进行缩放时，细节系数在第一张图中很难看出；

在第二张图中，图像被放大了8倍得以强调，这里级别为1的系数子图像边界的中等灰度填充，插入是为了调节变换系数子图的维数变换；

第三张图显示了细节绝对值的结果。

4 The Inverse Fast Wavelet Transform

类似于FWT正向变换，只不过顺序是颠倒了一下，而且，FWT正向变换中用到的downsampling下采样要转换为相应的upsampling上采样。

用公式表达如下：
$$\begin{gathered} I_D=W_{\psi }^D(j,m,n){\uparrow }^{m=k/2,k\ge 0}\odot h_{\psi }(m) \\ {I}_V=W_{\psi }^V(j,m,n){\uparrow }^{m=k/2,k\ge 0}\odot h_{\phi }(m) \\ {I}_H=W_{\psi }^H(j,m,n){\uparrow }^{m=k/2,k\ge 0}\odot h_{\psi }(m) \\ {I}_A=W_{\phi }(j,m,n){\uparrow }^{m=k/2,k\ge 0}\odot h_{\phi }(m) \\ {W}_{\phi }(j+1,m,n)=[[I_D+I_V]{\uparrow }^{n=k/2,k\ge 0}]\odot h_{\psi }(n)+[[I_H+I_A]{\uparrow }^{n=k/2,k\ge 0}]\odot h_{\phi }(n) \end{gathered}$$
公式也可简化为：
$$\begin{gathered} D\to hp、hp \\ V\to lp、hp \\ H\to hp、lp \\ A\to lp、lp \end{gathered}$$
一般计算IFWT时，需要使用function waverec2：

g = waverec2(C, S, wname)

其中g是得到的重构后的二维图像。要求的重构滤波器可以使用如下语句：

g = waverec2(C, S, LO_R, HI_R)

我们也可以根据IFWT的迭代图，自定义一个逆快速小波变换的函数function waveback

function [varargout] = waveback(c, s, varargin)
%WAVEBACK Computes inverse FWTs for multi-level decomposition [c, S] .
%   [VARARGOUT] = WAVEBACK(C, s, VARARGIN) performs a 2D N-level
%   partial or complete wavelet reconstruction of decomposition
%   structure [C, S].
% 
%   SYNTAX:
%   Y = WAVEBACK(C, S, ' WNAME');   Output inverse FWT matrix Y
%   Y = WAVEBACK(C, S, LR, HR);     using lowpass and highpass
%                                   reconstruction filters (LR and
%                                   HR) or filters obtained by
%                                   calling WAVEFILTER with 'WNAME'.
% 
%   [NC, NS] = WAVEBACK(C, S, 'WNAME', N);  Output new wavelet
%   [NC, NS] = WAVEBACK(C, s, LR, HR, N);   decomposition structure
%                                           [NC, NS] after N step
%                                           reconstruction.

% Check the input and output arguments for reasonableness.
error(nargchk(3, 5, nargin));
error(nargchk(1, 2, nargout));

if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
    error('C must be a row vector.');
end

if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
        ((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
    error('S must be a real, numeric two- or three column array.');
end

elements = prod(s, 2);

if (length(c) < elements(end)) || ...
        ~(elements(1) + 3 * sum(elements(2:end- 1 )) >= elements(end))
    error('[C, S] must be a standard wavelet decomposition structure.');
end

% Maximum levels in [C, S].
nmax = size(s, 1) - 2;

% Get third input parameter and init check flags.
wname = varargin{1};
filterchk = 0;
nchk = 0;

switch nargin
    case 3
        if ischar(wname)
            [lp, hp] = wavefilter(wname, 'r');
            n = nmax;
        else
            error('Undefined filter.');
        end
        if nargout ~= 1
            error('Wrong number of output arguments.');
        end
    case 4
        if ischar(wname)
            [lp, hp] = wavefilter(wname, 'r');
            n = varargin{2};
            nchk = 1;
        else
            lp = varargin{1};
            hp = varargin{2};
            filterchk = 1;
            n = nmax;
            if nargout ~= 1
                error('Wrong number of output arguments.');
            end
        end
    case 5
        lp = varargin{1};
        hp = varargin{2};
        filterchk = 1;
        n = varargin{3};
        nchk = 1;
    otherwise
        error('Improper number of input arguments');
end

fl = length(lp);
if filterchk
    if (ndims(lp) ~= 2) || ~isreal(lp) || ~isnumeric(lp) || ...
            (ndims(hp) ~= 2) || ~isreal(hp) || ~isnumeric(hp) || ...
            (fl ~= length(hp)) || (rem(fl, 2) ~= 0)
        error('LP and HP must be even and equal length real numeric filter vectors');
    end
end

if nchk && (~isnumeric(n) || ~isreal(n))
    error('N must be a real numeric');
end

if (n > nmax) || (n < 1)
    error('Invalid number(N) of reconstructions requested.');
end

if (n ~= nmax) && (nargout ~= 2)
    error('Not enough output arguments.');
end

% Init decomposition
nc = c;
ns = s;
nnmax = nmax;

for i = 1:n
    % Compute a new approximation
    a = symconvup(wavecopy('a', nc, ns), lp, lp, fl, ns(3, :)) + ...
        symconvup(wavecopy('h', nc, ns, nnmax), ...
                  hp, lp, fl, ns(3, :)) + ...
        symconvup(wavecopy('v', nc, ns, nnmax), ...
                  lp, hp, fl, ns(3, :)) + ...
        symconvup(wavecopy('d', nc, ns, nnmax), ...
                  hp, hp, fl, ns(3, :));
    % Update decomposition
    nc = nc(4 * prod(ns(1, :)) + 1:end);
    nc = [a(:)' nc];
    ns = ns(3:end, :);
    ns = [ns(1, :); ns];
    nnmax = size(ns, 1) - 2;
end

% For complete reconstructions, reformat output as 2-D.
if nargout == 1
    a = nc;
    nc = repmat(0, ns(1, :));
    nc(:) = a;
end

varargout{1} = nc;
if nargout == 2
    varargout{2} = ns;
end

% --------------------------------------------------------------------
function w = symconvup(x, f1, f2, fln, keep)
% Upsample rows and convolve columns with f1; upsample columns and convolve
% rows with f2; then extract center assuming symmetrical extension

% Process each "page" (i.e., 3rd index) of an extended 2-D array
% separately; if 'x' is 2-D, size(x, 3) = 1.

% Preallocate w.
zi = fln - 1:fln + keep(1) - 2;
zj = fln - 1:fln + keep(2) - 2;
w = zeros(numel(zi), numel(zj), size(x, 3));
for i = 1:size(x, 3)
    y = zeros([2 1] .* size(x(:, :, i)));
    y(1:2:end, :) = x(:,:,i);
    y = conv2(y, f1');
    z = zeros([1 2] .* size(y));
    z(:, 1:2:end) = y;
    z = conv2(z, f2);
    z = z(zi, zj);
    w(:, :, i) = z;
end

在自定义一个与waverec2做比较的函数function ifwtcompare

function [ratio, maxdiff] = ifwtcompare(f, n, wname)
%IFWTCOMPARE Compare waverec2 and waveback.
%   [RATIO, MAXDIFF] = IFWTCOMPARE(F, N, WNAME) compares the
%   operation of Wavelet Too1box function WAVEREC2 and cus tom
%   function WAVEBACK.
% 
%   INPUTS:
%       F           Image to transform and inverse trans form.
%       N           Number of   scales to compute.
%       WNAME       Wavelet to use.
% 
%   OUTPUTS:
%       RATIO       Execution time ratio (custom/ toolbox).
%       MAXDIFF     Maximum generated image difference.

% Compute the transform and get output and computation time for waverec2.
[c1, s1] = wavedec2(f, n, wname);
w1 = @() waverec2(c1, s1, wname);
reftime = timeit(w1);

% Compute the transform and get output and computation time for waveback.
[c2, s2] = wavefast(f, n, wname);
w2 = @() waveback(c2, s2, wname);
t2 = timeit(w2);

% Compare the results.
ratio = t2 / reftime;
diff = double(w1()) - w2();
maxdiff = abs(max(diff(:)));

比较结果如下：

>> [ratio, maxdifference] = ifwtcompare(f, 5, 'db4')
ratio =
    1.1269
maxdifference =
   3.4106e-13

可以说差距比较小的。

5 Wavelet in Image Processing

基本过程如下：

1. 计算一幅image的2-D Wavelet Transform二维小波变换
2. 修改transform coefficients变换系数
3. 计算Inverse transform反变换

Example Wavelet directionality and edge detection：

>> f = imread('Fig0707(a).tif');
>> subplot(221);imshow(f);
>> [c, s] = wavefast(f, 1, 'sym4');
>> subplot(222);wavedisplay(c, s, -6);
>> [nc, y] = wavecut('a', c, s);
>> subplot(223);wavedisplay(nc, s, -6);
>> edges = abs(waveback(nc, s, 'sym4'));
>> subplot(224);imshow(mat2gray(edges));

Example Wavelet-based image smoothing or blurring：

与Fourier Transform对应的Wavelet Transform是有效的图像平滑和图像模糊的处理方法。

>> subplot(321);imshow(f);
>> [c, s] = wavefast(f, 4, 'sym4');
>> subplot(322);wavedisplay(c, s, 20);
>> [c, g8] = wavezero(c, s, 1, 'sym4');
>> [c, s] = wavefast(f, 4, 'sym4');
>> [c, g8] = wavezero(c, s, 1, 'sym4');
>> subplot(323);imshow(g8);
>> [c, g8] = wavezero(c, s, 2, 'sym4');
>> subplot(324);imshow(g8);
>> [c, g8] = wavezero(c, s, 3, 'sym4');
>> subplot(325);imshow(g8);
>> [c, g8] = wavezero(c, s, 4, 'sym4');
>> subplot(326);imshow(g8);

Example Progressive reconstruction：

>> [c, s] = wavefast(f, 4, 'jpeg9.7');
>> wavedisplay(c, s, 8);
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 1
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 2
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 3
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));
>> [c, s] = waveback(c, s, 'jpeg9.7', 1); % Approximation 4
>> f = wavecopy('a', c, s);
>> figure; imshow(mat2gray(f));