非下采样剪切波变换_NSST函数库_matlab版本

/****************************************************************** * atrousc.c - Written by Arthur Cunha. This routine builds up on * zconv2D_OS.c written by Jason Laska * * Inputs: x - A 2D signal * h - 2D filter * m - separable upsampling matrix * * Outputs: y - 2D result of convolution with filter * upsampled by a m, only the 'valid' part is returned. * Similar to conv2(x,h,'valid'), where h is the upsampled * filter. * * * * Usage: y = zconv2D_O(x,h,m); * * Notes: This function does not actually upsample the filter, * it computes the convolution as if the filter had been * upsampled. This is the ultimate optimized version. * Further optimized for separable (diagonal) upsampling matrices. * * This is a MEX-FILE for matlab * /********************************************************/ #include "mex.h" #include <math.h> #define OUT plhs[0] #define SIGNAL prhs[0] #define FILTER prhs[1] #define MMATRIX prhs[2] #define LINPOS(row,col,collen) (row*collen)+col void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *FArray,*SArray,*outArray,*M; /* FArray - Filter coefficients SArray - Signal coefficients outArray - Output coefficients M - upsampling matrix */ int SColLength,SRowLength,FColLength,FRowLength,O_SColLength,O_SRowLength; int SFColLength,SFRowLength; int n1,n2,l1,l2,k1,k2,f1,f2, kk2, kk1; double sum; int M0,M3,sM0,sM3; SColLength = mxGetM(SIGNAL); SRowLength = mxGetN(SIGNAL); FColLength = mxGetM(FILTER); FRowLength = mxGetN(FILTER); SFColLength = FColLength-1; SFRowLength = FRowLength-1; FArray = mxGetPr(FILTER); SArray = mxGetPr(SIGNAL); M = mxGetPr(MMATRIX); M0 = (int)M[0]; M3 = (int)M[3]; sM0 = M0-1; sM3 = M3-1; O_SColLength = SColLength - M0*FColLength + 1; O_SRowLength = SRowLength - M3*FRowLength + 1; OUT = mxCreateDoubleMatrix(O_SColLength, O_SRowLength, mxREAL); outArray = mxGetPr(OUT); /* Convoluyion loop */ for (n1=0;n1<O_SRowLength;n1++){ for (n2=0;n2<O_SColLength;n2++){ sum=0; kk1 = n1 + sM0;; for (k1=0;k1<FRowLength;k1++){ kk2 = n2 + sM3; for (k2=0;k2<FColLength;k2++){ f1 = SFRowLength - k1; /* flipped index */ f2 = SFColLength - k2; sum+= FArray[LINPOS(f1,f2,FColLength)] * SArray[LINPOS(kk1,kk2,SColLength)]; kk2+=M3; } kk1+=M0; } outArray[LINPOS(n1,n2,O_SColLength)] = sum; } } return; }