%  本程序實現小波在矩陣方程求解中的應用
clear;clc;

%  構造小波矩陣
[h,g]= wfilters('db7','d');  %  分解低通和高通濾波器
N=512;               %  矩陣維數(大小為2的整數冪次)
L=length(h);         %  濾波器長度
rank_max=log2(N);    %  最大層數
rank_min=double(int8(log2(L)))+1;  %  最小層數
ww=1;   %  預處理矩陣

%  矩陣構造
for jj=rank_min:rank_max
    
    nn=2^jj;
    
    %  構造向量
    p1_0=sparse([h,zeros(1,nn-L)]);
    p2_0=sparse([g,zeros(1,nn-L)]);
    
    %  向量圓周移位
    for ii=1:nn/2
        p1(ii,:)=circshift(p1_0',2*(ii-1))';
        p2(ii,:)=circshift(p2_0',2*(ii-1))';
    end
    
    %  構造正交矩陣
    w1=[p1;p2];
    mm=2^rank_max-length(w1);
    w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);
    ww=ww*w;
    
    clear p1;clear p2;
end

norm(eye(N,N)-ww*ww.')
norm(eye(N,N)-ww.'*ww)
figure(1)
spy(ww)
title('小波預處理矩陣非零元素圖')

%  構造測試矩陣
for i=1:N;
    for j=1:N;
        z(i,j)=1/(abs(i-j)+1);
    end
end
% z=randn(N,N);

v=eye(N,1);
result=z\v;
 
%  小波變換(稀疏矩陣乘法)
z_t=ww*sparse(z)*ww.';
v_t=ww*sparse(v);
 
%  稀疏化矩陣
threshold=3*10^(-3);
z_t=z_t.*sparse(abs(z_t)>(threshold*norm(abs(z_t),1)));  %  大于閾值的元素保留,否則置零
 
%  求解小波域矩陣方程
j_t=z_t\v_t;
 
%  反小波變換
result_app=ww.'*j_t;
 
%  相對誤差
norm(result_app-result)/norm(result)
 
%  稀疏度圖
figure(2)
spy(z_t)
title('矩陣經小波變換稀疏化圖')