離散小波變換（DiscreteWavelet Transform, DWT）和離散傅立葉變換（Discrete FourierTransform, DFT）不一樣，在Matlab中確實有dwt函式，但它與一般書中講的DWT不一樣，dwt是基於Mallat（法國學者，音譯為馬拉特）演算法實現的，針對的離散時間訊號，而DWT指的是將連續小波變換（Continuous WaveletTransform, CWT）中的尺度引數a和時移引數b離散化，即將下式中的a和b離散化，但t仍然是連續的。

Mallat演算法框圖如下：

        小波變換跟普通的正交變換不同，它是一個多層分解。如Mallat演算法框圖所示，一個長為N的訊號x

         不同的小波母函式對應的小波基不同，在特定小波基下最多可以分解多少層在Matlab中可以使用函式wmaxlev去求解。

          在Mallat演算法框圖中，多層分解的過程相當於濾波的過程，例如由x0分別得到d1和x1相當於讓x0分別通過濾波器h1和h0，然後再下采樣，每兩個點抽取一個點。濾波的過程實際上就是卷積了，但這裡就有一個問題：N點長的序列與M點長的序列線性卷積後長度為N+M-1，並不等於N，也就是說下采樣後每兩點抽取一點得到的結果並不是N/2，而且線性卷積的四個步驟（反轉、平移、相乘、相加）在計算卷積結果兩側的點時，反轉平移後兩個卷積序列只有部分點對應上，到底如何處理邊界的這些點？這裡就涉及到了訊號的擴充套件，Matlab支援多種擴充套件模式，可以使用函式dwtmode查詢、更改。值得注意的是，只有Periodization模式分解結果長度為N/2，其它模式長度均為floor((N+M-1)/2)，為了保證小波分解後總體長度不變，因此這裡必須採用Periodization模式；Matlab預設訊號擴充套件模式為Symmetrization(half-point)，可以通過dwtmode(‘per’)設定為Periodization模式。

        在Matlab中，dwt只能實現單層分解，更強大功能的函式是wavedec函式，可以設定分解層數。前面說了小波分解實際上一個濾波的過程，濾波器的係數的可通過函式wfilters得到，濾波可以通過卷積實現，卷積可以表示為矩陣運算，因此小波分解可以表示為輸出向量等於小波變換矩陣乘以輸入向量。在這篇部落格裡，我要得到一個小波變換矩陣，使其等價於Matlab自帶的函式dwt和wavedec功能。我的Matlab版本是Version7.9.0.529(R2009b)，不知道不同的版本是否會不一樣。

        在Mallat演算法框圖中，G₁(n)=x(n)*g(n)，H₁(n)= x(n

設x(n)的長度為N，h(n)的長度為M，則H₁(n)的長度為N+M-1(假設N>M)，以上卷積關係寫為矩陣形式：（H₁(n)以y(n)表示）

卷積之後還要抽樣，每兩個點抽一個點，這裡dwtmode設定為Periodization模式，因此抽樣只有N/2個點，實際上就是對上面由h構成的矩陣的N+M-1行當中抽取N/2行，但抽取哪N/2行呢？通過驗證，在MATLAB中是這樣子做的，將矩陣的前M/2-1行和後M/2-1行略掉，這樣還剩N+M-1-( M/2-1)*2=N+1行，然後抽取其中的偶數行即可得到N/2行，這裡假設了M為偶數，而在小波變換中濾波器的長度(h和g的長度)就是偶數。

       由於這裡dwtmode設定為Periodization模式，即對x進行週期延拓（注：這裡x的長度取為2的整數次冪，因此與Periodized Padding模式(可通過dwtmode(‘ppd’)設定)是相同的），因此矩陣每一行實際上應該是h各元素的一個迴圈移位，因為卷積要反轉移位，所以每一行是h各元素倒序的迴圈移位，當然，由於h的長度小於x，因為要補零的。

        從由h構成的矩陣抽取N/2行，同理，再從由g構成的矩陣抽取N/2行，上下放在一起可以組成一個N×N的矩陣，即我們要找的小波變換矩陣。

       另外，為了實現多層分解，只要用一個矩陣連乘就可以了。從x得到D₁和A₁的過程（即第一層分解）可以用矩陣表示為：

第二層分解是由A₁得到D₂和A₂的過程，可用矩陣表示如下：

若把第一層分解和第二層分解寫在一起，可用分塊矩陣的概念表示如下：

第三層分解以此規律類推即可。

         理論就闡述這麼多，最重要的是給出Matlab程式：

function [ ww ] = dwtmtx( N,wtype,wlev )
%DWTMTX Discrete wavelet transform matrix
%   This function generates the transform matrix ww according to input
%   parameters N,wtype,wlev .
%Detailed explanation goes here
%   N is the dimension of ww
%   wtype is the wavelet type
%   wlev is the number of decomposition level
%NOTE: The extension mode must be Periodization('per')
[h,g]= wfilters(wtype,'d');         %Decomposition low&high pass filter
L=length(h);                        %Filter length
h_1 = fliplr(h);                    %Flip matrix left to right
g_1 = fliplr(g);
loop_max = log2(N);
loop_min = double(int8(log2(L)))+1;
if wlev>loop_max-loop_min+1
    fprintf('\nWaring: wlev is too big\n');
    fprintf('The biggest wlev is %d\n',loop_max-loop_min+1);
    wlev = loop_max-loop_min+1;
end
ww=1;
for loop = loop_max-wlev+1:loop_max
    Nii = 2^loop;
    p1_0 = [h_1 zeros(1,Nii-L)];
    p2_0 = [g_1 zeros(1,Nii-L)];
    p1 = zeros(Nii/2,Nii);
    p2 = zeros(Nii/2,Nii);
    for ii=1:Nii/2
        p1(ii,:)=circshift(p1_0',2*(ii-1)+1-(L-1)+L/2-1)';
        p2(ii,:)=circshift(p2_0',2*(ii-1)+1-(L-1)+L/2-1)';
    end
    w1=[p1;p2];
    mm=2^loop_max-length(w1);
    w=[w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)];
    ww=ww*w;
    clear p1;clear p2;
end

%The end!!!
end

%驗證函式dwtmtx的正確性
clear all;close all;clc;
N = 2^7;
% 'db1' or 'haar', 'db2', ... ,'db10', ... , 'db45'
% 'coif1', ... , 'coif5'
% 'sym2', ... , 'sym8', ... ,'sym45'
% 'bior1.1', 'bior1.3', 'bior1.5'
% 'bior2.2', 'bior2.4', 'bior2.6', 'bior2.8'
% 'bior3.1', 'bior3.3', 'bior3.5', 'bior3.7'
% 'bior3.9', 'bior4.4', 'bior5.5', 'bior6.8'
% 'rbio1.1', 'rbio1.3', 'rbio1.5'
% 'rbio2.2', 'rbio2.4', 'rbio2.6', 'rbio2.8'
% 'rbio3.1', 'rbio3.3', 'rbio3.5', 'rbio3.7'
% 'rbio3.9', 'rbio4.4', 'rbio5.5', 'rbio6.8'
wtype = 'rbio6.8';
wlev_max = wmaxlev(N,wtype);
if wlev_max == 0
    fprintf('\nThe parameter N and wtype does not match!\n');
end
dwtmode('per');
for wlev = 1:wlev_max
    ww = dwtmtx(N,wtype,wlev);
    x = randn(1,N);
    y1 = (ww*x')';
    [y2,y2l] = wavedec(x,wlev,wtype);
    y_err = sum((y1-y2).*(y1-y2));
    fprintf('wlev = %d: y_err = %f\n',wlev,y_err);
end