Pardon the ten-year delay: To build 1D and 2D Haar-wavelet basis functions. Unadulterated but functional code (X is just a test signal): %macro o(a,b); proc iml; R=&&h&a; C=&&h&b; X=shape(1:R*C,R,C); Y=colvec(X); %macro n(n); S&n=0; L&n=log2(&n); D&n=2**(-L&n/2); H&n=1; A={1,1}; B={1,-1}; do j=1 to L&n; k=j-1; I=I(2**k); H&n=H&n@A||I@B; *Haar builder; D&n=D&n||shape(2**(-(L&n-k)/2),1,2**k); S&n=S&n||2**j-1; end; *Scale builder; D&n=diag(D&n); H&n=H&n*D&n; %mend; %n(R); %n(C); H=HR@HC; print H; *End Haar; *2D@2D; HXH=HR`*X*HC; IX=HR*HXH*HC`; *Haar transform and recover X; *2D@1D; HHY=H*Y; IY=H`*HHY; *Vectorize for use in penalized regression;