小子域滤波程序

clear all
close all

% load hulin.mat;
T=load（‘hulin.dat‘）;

px=T（:1）;
py=T（:2）;
pg=T（:3）;
XXX=reshape（px105115）;
XX=XXX‘;
YYY=reshape（py105115）;
YY=YYY‘;
gg=reshape（pg105115）;
g=gg‘;

[mn]=size（g）;
x=0:2:n-1;y=0:2:m-1; xstep=x（2）-x（1）;ystep=y（2）-y（1）;   NN=30;


mt=floor（log2（m））+2;%扩边处理；
M=2^mt;
nt=floor（log2（n））+2;
N=2^nt;
LL1=floor（（M-m）/2）;
LL2=M-m-LL1;
L1=floor（（N-n）/2）;
L2=N-n-L1;
gk=zeros（mN）;gk（:L1+1:1:L1+n）=g;
for k=1:1:L1
    for row=1:1:m
        gk（rowk）=gk（rowL1+1）*（1+cos（pi*（L1+1-k）/L1））/2;
    end
end
for k=1:1:L2
    for row=1:1:m
    gk（rown+L1+k）=gk（rown+L1）*（1+cos（（pi*k）/（L2）））/2;
    end
end
gk2=zeros（NM）;gk2（:LL1+1:1:LL1+m）=gk‘;
for k=1:1:LL1
    for row=1:1:N
    gk2（rowk）=gk2（rowLL1+1）*（1+cos（pi*（LL1+1-k）