德国开元华人社区开元周游

标题: 做实习遇到的问题 [打印本页]

作者: loveless 时间: 10.6.2006 16:16
这次实习做的是关于用Wavelet变换来分析EMG的信号。观察Aktivitaet和它的关系。 现在碰到个问题，对于已有的Wavelet的C的源代码如何变成Simulink里的S-Function 呢？请懂的人说一下。十分感谢！

作者: loveless 时间: 10.6.2006 16:38
那你能不能给我看看C的源代码？很难理解，看不懂。

作者: loveless 时间: 11.6.2006 22:39
#define N0 128 #include "stdio.h" #include "stdlib.h" #include "math.h" #include "string.h" void db4(double *h,double *g,double *hh,double *gg); void wd(int N,double *h,double *g,double *c0,double *c,double *d); void wr(int N,double *h,double *g,double *c, double *d,double *cd); void main() { double fk[N0],c0[N0],c[N0],d[N0]; double h[8],g[8],hh[8],gg[8]; float fk0[N0]; FILE *fp; int i,k,j,n,l,N; fp=fopen("wdata.dat","rt"); fscanf(fp,"%d",&N); for(k=0;kfclose(fp); db4(h,g,hh,gg); for(k=0;kc0[k]=fk0[k]; c[k]=0; d[k]=0; } wd(N,hh,gg,c0,c,d); wr(N,hh,gg,c,d,c0); for(k=0;kreturn; } void wd(int N,double *h,double *g,double *c0,double *c,double *d) /* wavelet decomposition */ { int k,n,k2,l; double ck,dk; for(k=0;kck=0.0; dk=0.0; for(l=0;l<8;l++) { n=k+l; ck+=c0[n%N]*h[l]; dk+=c0[n%N]*g[l]; } c[k]=ck; d[k]=dk; } for(k=0;kk2=2*k; c0[k]=c[k2]; c0[N/2+k]=d[k2]; } return; } void wr(int N,double *h,double *g,double *c,double *d,double *c0) /* wavelet reconstruction */ { int k,n,l,k2; double ck,cn,dn; for(k=0;kk2=2*k; c[k2]=c0[k]; c[k2+1]=0; d[k2]=c0[N/2+k]; d[k2+1]=0; } for(k=0;kfor(k=0;kck=0.0; for(l=0;l<8;l++) { n=k-l; cn=c[(N+n)%N]; dn=d[(N+n)%N]; ck+=cn*h[l]+dn*g[l]; } c0[k]=ck; } return; } void db4(double *h,double *g,double *hh,double *gg) /* Daubechies 4 wavelet */ { int k,isgn; h[7]=-0.0105974017850890; h[6]= 0.0328830116668852; h[5]= 0.0308413818355607; h[4]=-0.1870348117190931; h[3]=-0.0279837694168599; h[2]= 0.6308807679398597; h[1]= 0.7148465705529154; h[0]= 0.2303778133088964; isgn=1; for(k=0;k<8;k++) { gg[k]=isgn*h[7-k]; isgn=-isgn; } for(k=0;k<8;k++) { g[k]=gg[7-k]; hh[k]=h[7-k]; } return; } float fun(float x) { float pi=3.1415926; float yx=30*exp(-x/40)*sin(2*pi*x/40); return(yx); } 这个用的是 DB4小波，周期延拓，可以实现精确重构的。

作者: loveless 时间: 11.6.2006 22:40
谢谢你给我看看！

作者: loveless 时间: 12.6.2006 22:02
谢谢阿，我再仔细看看。

作者: loveless 时间: 13.6.2006 16:08
问下斑竹，这个就是所谓的短波变换吗？ 我现在做实习里面的前提就是先要用Simulink进行仿真。 输入是时间和频率，或者只是频率。因为是离散的时间，时间间隔是常数。 那么如果用以上的C源代码。怎么变成Simulink里的S函数呢？ 数据是从From Workspace输入，经过变换和Scope之后，还要回到Workspace。就是to Workspace。而我现在就是要通过这个变化得到那个Spektrogramm。

作者: loveless 时间: 16.6.2006 13:46
的确是这样，现在就是这个S函数该怎么调整？能不能再说的详细一点啊？谢谢！

作者: 行走城市的猫 时间: 17.6.2006 08:43
eisenstange 好人加能人。 我文科生，一窍不懂得说。但是，非常欣赏这样相互帮助的态度。

作者: 行走城市的猫 时间: 17.6.2006 10:59
<!--emo&

--><img src='style_emoticons/<#EMO_DIR#>/smile.gif' border='0' style='vertical-align:middle' alt='smile.gif' /> <div class='quotetop'>QUOTE(eisenstange @ 17.06.2006, 09:07 )</div><div class='quotemain'>万事开头难，每个人都有一开始的时候，作为过来人，因该尽可能的扶一下，谁叫现在企业对毕业生的要求都那么高呢，大家学的好，对个人和国家社会都有利。牺牲个把小时，值得。有时候自己也能从中学到一些东西。利己利人何乐而不为。 [right][snapback]1008611[/snapback][/right] </div> 如果大家都能这么团结一致，共同研究，那么，我们中国学生勤奋好学，紧密联系会让周围的外国人敬佩的。

作者: loveless 时间: 17.6.2006 16:15
今天早上刚和一起做的人讨论好回来。就看到斑竹如此详细的解释。真是很感谢。马上试下。

作者: loveless 时间: 17.6.2006 16:19
Mex-Setup是在命令行里输入吗？

作者: loveless 时间: 17.6.2006 16:19
>> Mex-Setup Warning: Could not find an exact (case-sensitive) match for 'Mex'. C:\

rogram Files\MATLAB704\toolbox\matlab\general\mex.m is a case-insensitive match and will be used instead. You can improve the performance of your code by using exact name matches and we therefore recommend that you update your usage accordingly. Alternatively, you can disable this warning using warning('off','MATLAB:dispatcher:InexactMatch'). Select a compiler: [1] Lcc C version 2.4 in C:\

ROGRAM FILES\MATLAB704\sys\lcc [0] None Compiler: 1 Try to update options file: C:\Documents and Settings\JuliusWoo\Application Data\MathWorks\MATLAB\R14\mexopts.bat From template: C:\

ROGRAM FILES\MATLAB704\BIN\WIN32\mexopts\lccopts.bat Done . . . Usage: MEX [option1 ... optionN] sourcefile1 [... sourcefileN] [objectfile1 ... objectfileN] [libraryfile1 ... libraryfileN] or (to build an Ada S-function): MEX [-v] [-g] -ada <sfcn>.ads Use the -help option for more information, or consult the MATLAB API Guide. C:\

ROGRAM FILES\MATLAB704\BIN\WIN32\MEX.PL: Error: No file names given. ??? Undefined function or variable 'Setup'.

作者: loveless 时间: 17.6.2006 16:23
<div class='quotetop'>QUOTE(eisenstange @ 17.06.2006, 16:21 )</div><div class='quotemain'>可能是你的C编译器没有装好，你的机器上用的是什么编译器？ VC还是GCC？ [right][snapback]1008958[/snapback][/right] </div> 我也不清楚，这个是怎么看？

作者: loveless 时间: 17.6.2006 16:25
是不是直接把C代码，粘贴在下面。就是运行了Mex-Setup后的命令行下面？

作者: loveless 时间: 17.6.2006 16:31
// Onur G. Guleryuz 1995, 1996, 1997, // University of Illinois at Urbana-Champaign, // Princeton University, // Polytechnic University. #include <stdio.h> #include <stdlib.h> #include "wav_filters_extern.h" #include "wav_gen.h" #include "alloc.h" // Extend data by ofs pixels on both ends using periodic // symmetry. data should be "pointing in" to allow up to data[-ofs] // and data[N-1+ofs]. // ofs>=N is handled. void extend_periodic(float *data,int N,int ofs) { int i,k; k=N-1; for(i=1;i<=ofs;i++) { data[-i]=data[N-i]; data[k+i]=data[i-1]; } } // Extend data by ofs pixels on both ends using mirror // symmetry. data should be "pointing in" to allow up to data[-ofs] // and data[N-1+ofs]. // // There are two possible "mirrors", i.e., data[-1]=data[1] (mirror on 0) // or data[-1]=data[0] (mirror "between" -1 and 0). phase is either 0, 1, or 2 // to select among these cases. See if(phase==... below. void extend_mirror(float *data,int N,int ofs,int phase) { int i,k; int phase1,phase2; if((phase<0)||(phase>2)) { printf("extend_mirror: illegal phase\n"); exit(1); } if(phase==2){ // At left boundary the mirror is on 0. phase1=0; // At right boundary the mirror is on N-1. phase2=1; } else { // phase==0, at left mirror on 0, and at right it's between. // phase==1, at left mirror between, and at right it's on N-1. phase1=phase2=phase; } k=N-1; if(N==1) { for(i=1;i<=ofs;i++) { data[-i]=data[0]; data[k+i]=data[0]; } } else for(i=1;i<=ofs;i++) { data[-i]=data[i-phase1]; data[k+i]=data[N-phase2-i]; } } // This routine is used in evaluating forward wavelet transforms. // // data of dimension N contains the input. // float *filter contains the Nf filter taps. // ofs is the decimation factor. // // filtering starts at beg and continues until N+beg, yielding N/ofs samples. // data needs to be padded suitably at both ends. // This hoopla about the starting point is in place to ensure // correct boundary processing when evaluating wavelet transforms. // // coef returns the result of the filtering and decimation. // // Viewed as a transform the basis function that generates the // coefficient of index 0 is "<- fliplr(filter) -> (beg)", // i.e., the scalar product obtained by positioning // the flipped filter (formed into a row vector) // with the right side positioned on top of sample beg in data. int filter_n_decimate(float *data,float *coef,int N,float *filter, int Nf,int beg,int ofs) { int i,j,k; float temp; k=0; for(i=beg;i<N+beg;i+=ofs) { temp=0; for(j=i;j>i-Nf;j--) temp+=data[j]*filter[i-j]; coef[k]=temp; k++; } return(k); } // This routine is used in evaluating inverse wavelet transforms. // // Upsample and filter to implement a synthesis filter bank. // Basically the inverse of filter_n_decimate() above. // No actual upsampling etc., to avoid zero multiplies. // // coef is the input and contains coefficients resulting from // filtering and decimation by ofs. // // Because this routine is typically called twice for inverse wavelet // transforms, once for the low band and once for high, data is incremented via +=. // Hence data should initially be set to 0. // // Viewed as a transform, filter_n_decimate above starts by putting the // "flipped" forward filter at begf, i.e., <- fliplr(ffilter) ->(begf) // The resulting coefficient is at index 0. // Here we are upsampling by ofs and filtering with the inv. filter // time advanced by fbeg so the coefficient at index 0, multiplies the basis function // (-fbeg)<-(ifilter)->(Nf-1-fbeg=beg). // // How is that for confusing? All of this makes sense, // when you write it down carefully. But who has the time ...<!--emo&

--><img src='style_emoticons/<#EMO_DIR#>/smile.gif' border='0' style='vertical-align:middle' alt='smile.gif' /> void upsample_n_filter(float *coef,float *data,int N,float *filter, int Nf,int beg,int ofs) { int i,j,l,n,p,fbeg; float temp; fbeg=Nf-beg-1; for(i=0;i<N;i++) { l=0; n=(i+fbeg)%ofs; p=(i+fbeg)/ofs; temp=0; for(j=n;j<Nf;j+=ofs) { temp+=coef[p-l]*filter[j]; l++; } // += for successive calls for low and high bands. data+=temp; } } // Used in forward wavelet trf. // // All rows of the image are run through the dual filterbank // given by lp and hp. Results are returned in image "in place". // low freq. coefficients start at 0 and end at N/2, where high freq. // info starts. void filt_n_dec_all_rows(float **image,int Ni,int Nj,float *lp,int Nl, float *hp,int Nh) { int i,j,ext,ofs; float *data; ext=max(Nl,Nh); data=allocate_1d_float((Nj+2*ext),0); // Point in for required extensions. data+=ext; // offset for the location where high band coefficients // should be copied. ofs=Nj>>1; for(i=0;i<Ni;i++) { // Copy row. for(j=0;j<Nj;j++) data[j]=image[j]; // Extend. if(PS) extend_periodic(data,Nj,ext); else { // Mirrors at end points. extend_mirror(data,Nj,ext,2); } filter_n_decimate(data,image,Nj,lp,Nl,begflp,2); filter_n_decimate(data,image+ofs,Nj,hp,Nh,begfhp,2); } data-=ext; free((void *)data); } // Used in forward wavelet trf. // // Index swapped version of filt_n_dec_all_rows. void filt_n_dec_all_cols(float **image,int Ni,int Nj,float *lp,int Nl, float *hp,int Nh) { int i,j,ext,ofs; float *data,*data2; ext=max(Nl,Nh); data=allocate_1d_float((Ni+2*ext),0); // Point in for required extensions. data+=ext; // Need second array since cannot access columns directly via pointers. data2=allocate_1d_float(Ni,0); // offset for the location where high band coefficients // should be copied. ofs=Ni>>1; for(j=0;j<Nj;j++) { // Copy column. for(i=0;i<Ni;i++) data=image[j]; // Extend. if(PS) extend_periodic(data,Ni,ext); else { // Mirrors at end points. extend_mirror(data,Ni,ext,2); } filter_n_decimate(data,data2,Ni,lp,Nl,begflp,2); filter_n_decimate(data,data2+ofs,Ni,hp,Nh,begfhp,2); for(i=0;i<Ni;i++) image[j]=data2; } data-=ext; free((void *)data); free((void *)data2); } // Used in inverse wavelet trf. // // All rows of the image are run through the dual synthesis filterbank // given by lp and hp having number of taps Nl and Nh respectively. // Results are returned in image "in place". // // Here we actaully need to know the applied shift to the transform, // hence the passed parameters lev and shift_arr. void ups_n_filt_all_rows(float **image,int Ni,int Nj,float *lp,int Nl, float *hp,int Nh,int lev,int *shift_arr) { int i,j,k,ext,ofs; float *data1,*data2; ext=max(Nl,Nh); ofs=Nj>>1; data1=allocate_1d_float((ofs+2*ext),0); data2=allocate_1d_float((ofs+2*ext),0); data1+=ext;data2+=ext; for(i=0;i<Ni;i++) { for(j=0;j<ofs;j++) { k=j+ofs; // low pass and high pass coefficients. data1[j]=image[j];image[j]=0; data2[j]=image[k];image[k]=0; } // Take care of the extension at the boundaries. if(PS) { extend_periodic(data1,ofs,ext); extend_periodic(data2,ofs,ext); } else { // shift_arr[lev]=0 OR 1. // The symmetric banks used are positioned so that on the left side // the forward lowpass filter has its point of symmetry exactly on top of 0 // and forward high pass filter has its point of symmetry exactly on top of 1. // This is coordinated by calling filter_n_decimate() with the proper value of beg. // // The above is reversed at the right side assuming decimation by 2. // Hence shift_arr[lev]==0 => at left boundary the mirror is at 0, // at right boundary it's between, see extend_mirror(). // If shift_arr[lev]==1 then the positions for the low and high are reversed // since the filters are shifted so if // shift_arr[lev]==1 => at left boundary the mirror is between, // at right boundary it's at 0. extend_mirror(data1,ofs,ext,shift_arr[lev]); extend_mirror(data2,ofs,ext,1-shift_arr[lev]); } // invert low pass. upsample_n_filter(data1,image,Nj,lp,Nl,begilp,2); // invert high pass. upsample_n_filter(data2,image,Nj,hp,Nh,begihp,2); } data1-=ext;data2-=ext; free((void *)data1); free((void *)data2); } // Used in inverse wavelet trf. // // Index swapped version of ups_n_filt_all_rows. void ups_n_filt_all_cols(float **image,int Ni,int Nj,float *lp,int Nl, float *hp,int Nh,int lev,int *shift_arr) { int i,j,k,ext,ofs; float *data1,*data2,*data3; ext=max(Nl,Nh); ofs=Ni>>1; data1=allocate_1d_float((ofs+2*ext),0); data2=allocate_1d_float((ofs+2*ext),0); data1+=ext;data2+=ext; // Need third array since cannot access columns directly via pointers. data3=allocate_1d_float(Ni,0); for(j=0;j<Nj;j++) { for(i=0;i<ofs;i++) { k=i+ofs; // low pass and high pass coefficients. data1=image[j]; data2=image[k][j]; data3=data3[k]=0; } if(PS) { extend_periodic(data1,ofs,ext); extend_periodic(data2,ofs,ext); } else { extend_mirror(data1,ofs,ext,shift_arr[lev]); extend_mirror(data2,ofs,ext,1-shift_arr[lev]); } upsample_n_filter(data1,data3,Ni,lp,Nl,begilp,2); upsample_n_filter(data2,data3,Ni,hp,Nh,begihp,2); for(i=0;i<Ni;i++) image[j]=data3; } data1-=ext;data2-=ext; free((void *)data1); free((void *)data2); free((void *)data3); } 
作者: loveless 时间: 17.6.2006 16:33
以上的是Wavelet Basic的C代码。是不是把#Include <mex.h> 放在// Polytechnic University. #include <stdio.h> 之间？ 这个源代码本身已经包含了很多的Include,该怎么处理呢？ 
作者: loveless 时间: 17.6.2006 16:37
mex wav_basic.c wav_basic.obj .text: undefined reference to '_allocate_1d_float' wav_basic.obj .text: undefined reference to '_PS' wav_basic.obj .text: undefined reference to '_begflp' wav_basic.obj .text: undefined reference to '_begfhp' wav_basic.obj .text: undefined reference to '_begilp' wav_basic.obj .text: undefined reference to '_begihp' Specified export _mexFunction is not defined Missing exports. Aborting C:\ROGRAM FILES\MATLAB704\BIN\WIN32\MEX.PL: Error: Link of 'wav_basic.dll' failed. ??? Error using ==> mex Unable to complete successfully 
作者: loveless 时间: 17.6.2006 16:38
<div class='quotetop'>QUOTE(loveless @ 17.06.2006, 16:31 )</div><div class='quotemain'>// Onur G. Guleryuz 1995, 1996, 1997, // University of Illinois at Urbana-Champaign, // Princeton University, // Polytechnic University. #include <stdio.h> #include <stdlib.h> #include "wav_filters_extern.h" #include "wav_gen.h" #include "alloc.h" // Extend data by ofs pixels on both ends using periodic // symmetry. data should be "pointing in" to allow up to data[-ofs] // and data[N-1+ofs]. // ofs>=N is handled. void extend_periodic(float *data,int N,int ofs) { int i,k; k=N-1; for(i=1;i<=ofs;i++) {   data[-i]=data[N-i];   data[k+i]=data[i-1]; } } // Extend data by ofs pixels on both ends using mirror // symmetry. data should be "pointing in" to allow up to data[-ofs] // and data[N-1+ofs]. // // There are two possible "mirrors", i.e., data[-1]=data[1] (mirror on 0) // or data[-1]=data[0] (mirror "between" -1 and 0). phase is either 0, 1, or 2 // to select among these cases. See if(phase==... below. void extend_mirror(float *data,int N,int ofs,int phase) { int i,k; int phase1,phase2; if((phase<0)||(phase>2)) {   printf("extend_mirror: illegal phase\n");   exit(1); } if(phase==2){   // At left boundary the mirror is on 0.   phase1=0;   // At right boundary the mirror is on N-1.   phase2=1; } else {   // phase==0, at left mirror on 0, and at right it's between.   // phase==1, at left mirror between, and at right it's on N-1.   phase1=phase2=phase; } k=N-1;   if(N==1) {   for(i=1;i<=ofs;i++) {   data[-i]=data[0];   data[k+i]=data[0];   } } else   for(i=1;i<=ofs;i++) {   data[-i]=data[i-phase1];   data[k+i]=data[N-phase2-i];   } } // This routine is used in evaluating forward wavelet transforms. // // data of dimension N contains the input. // float *filter contains the Nf filter taps. // ofs is the decimation factor. // // filtering starts at beg and continues until N+beg, yielding N/ofs samples. // data needs to be padded suitably at both ends. // This hoopla about the starting point is in place to ensure // correct boundary processing when evaluating wavelet transforms. // // coef returns the result of the filtering and decimation. // // Viewed as a transform the basis function that generates the // coefficient of index 0 is "<- fliplr(filter) -> (beg)", // i.e., the scalar product obtained by positioning // the flipped filter (formed into a row vector) // with the right side positioned on top of sample beg in data. int filter_n_decimate(float *data,float *coef,int N,float *filter,       int Nf,int beg,int ofs) { int i,j,k; float temp; k=0; for(i=beg;i<N+beg;i+=ofs) {   temp=0;   for(j=i;j>i-Nf;j--)   temp+=data[j]*filter[i-j];   coef[k]=temp;   k++; } return(k); } // This routine is used in evaluating inverse wavelet transforms. // // Upsample and filter to implement a synthesis filter bank. // Basically the inverse of filter_n_decimate() above. // No actual upsampling etc., to avoid zero multiplies. // // coef is the input and contains coefficients resulting from // filtering and decimation by ofs. // // Because this routine is typically called twice for inverse wavelet // transforms, once for the low band and once for high, data is incremented via +=. // Hence data should initially be set to 0. // // Viewed as a transform, filter_n_decimate above starts by putting the // "flipped" forward filter at begf, i.e.,  <- fliplr(ffilter) ->(begf) // The resulting coefficient is at index 0. // Here we are upsampling by ofs and filtering with the inv. filter // time advanced by fbeg so the coefficient at index 0, multiplies the basis function // (-fbeg)<-(ifilter)->(Nf-1-fbeg=beg). // // How is that for confusing? All of this makes sense, // when you write it down carefully. But who has the time ...<img src='style_emoticons/<#EMO_DIR#>/smile.gif' border='0' style='vertical-align:middle' alt='smile.gif' /> void upsample_n_filter(float *coef,float *data,int N,float *filter,       int Nf,int beg,int ofs) { int i,j,l,n,p,fbeg; float temp; fbeg=Nf-beg-1; for(i=0;i<N;i++) {   l=0;   n=(i+fbeg)%ofs;   p=(i+fbeg)/ofs;   temp=0;   for(j=n;j<Nf;j+=ofs) {   temp+=coef[p-l]*filter[j];   l++;   }   // += for successive calls for low and high bands.   data+=temp; } } // Used in forward wavelet trf. // // All rows of the image are run through the dual filterbank // given by lp and hp. Results are returned in image "in place". // low freq. coefficients start at 0 and end at N/2, where high freq. // info starts. void filt_n_dec_all_rows(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh) { int i,j,ext,ofs; float *data; ext=max(Nl,Nh); data=allocate_1d_float((Nj+2*ext),0); // Point in for required extensions. data+=ext; // offset for the location where high band coefficients // should be copied. ofs=Nj>>1; for(i=0;i<Ni;i++) {     // Copy row.   for(j=0;j<Nj;j++)   data[j]=image[j];   // Extend.   if(PS)   extend_periodic(data,Nj,ext);   else {   // Mirrors at end points.   extend_mirror(data,Nj,ext,2);   }   filter_n_decimate(data,image,Nj,lp,Nl,begflp,2);   filter_n_decimate(data,image+ofs,Nj,hp,Nh,begfhp,2);   } data-=ext; free((void *)data); } // Used in forward wavelet trf. // // Index swapped version of filt_n_dec_all_rows. void filt_n_dec_all_cols(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh) { int i,j,ext,ofs; float *data,*data2; ext=max(Nl,Nh); data=allocate_1d_float((Ni+2*ext),0); // Point in for required extensions. data+=ext; // Need second array since cannot access columns directly via pointers. data2=allocate_1d_float(Ni,0); // offset for the location where high band coefficients // should be copied. ofs=Ni>>1; for(j=0;j<Nj;j++) {     // Copy column.   for(i=0;i<Ni;i++)   data=image[j];   // Extend.   if(PS)   extend_periodic(data,Ni,ext);   else {   // Mirrors at end points.   extend_mirror(data,Ni,ext,2);   }   filter_n_decimate(data,data2,Ni,lp,Nl,begflp,2);   filter_n_decimate(data,data2+ofs,Ni,hp,Nh,begfhp,2);   for(i=0;i<Ni;i++)   image[j]=data2;   } data-=ext; free((void *)data); free((void *)data2); } // Used in inverse wavelet trf. // // All rows of the image are run through the dual synthesis filterbank // given by lp and hp having number of taps Nl and Nh respectively. // Results are returned in image "in place". // // Here we actaully need to know the applied shift to the transform, // hence the passed parameters lev and shift_arr. void ups_n_filt_all_rows(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh,int lev,int *shift_arr) { int i,j,k,ext,ofs; float *data1,*data2; ext=max(Nl,Nh); ofs=Nj>>1; data1=allocate_1d_float((ofs+2*ext),0); data2=allocate_1d_float((ofs+2*ext),0); data1+=ext;data2+=ext; for(i=0;i<Ni;i++) {   for(j=0;j<ofs;j++) {   k=j+ofs;   // low pass and high pass coefficients.   data1[j]=image[j];image[j]=0;   data2[j]=image[k];image[k]=0;   }   // Take care of the extension at the boundaries.   if(PS) {   extend_periodic(data1,ofs,ext);   extend_periodic(data2,ofs,ext);   }   else {   // shift_arr[lev]=0 OR 1.   // The symmetric banks used are positioned so that on the left side   // the forward lowpass filter has its point of symmetry exactly on top of 0   // and forward high pass filter has its point of symmetry exactly on top of 1.   // This is coordinated by calling filter_n_decimate() with the proper value of beg.   //   // The above is reversed at the right side assuming decimation by 2.   // Hence shift_arr[lev]==0 => at left boundary the mirror is at 0,   // at right boundary it's between, see extend_mirror().   // If shift_arr[lev]==1 then the positions for the low and high are reversed   // since the filters are shifted so if   // shift_arr[lev]==1 => at left boundary the mirror is between,   // at right boundary it's at 0.   extend_mirror(data1,ofs,ext,shift_arr[lev]);   extend_mirror(data2,ofs,ext,1-shift_arr[lev]);   }   // invert low pass.   upsample_n_filter(data1,image,Nj,lp,Nl,begilp,2);   // invert high pass.   upsample_n_filter(data2,image,Nj,hp,Nh,begihp,2); } data1-=ext;data2-=ext; free((void *)data1); free((void *)data2); } // Used in inverse wavelet trf. // // Index swapped version of ups_n_filt_all_rows. void ups_n_filt_all_cols(float **image,int Ni,int Nj,float *lp,int Nl,       float *hp,int Nh,int lev,int *shift_arr) { int i,j,k,ext,ofs; float *data1,*data2,*data3; ext=max(Nl,Nh); ofs=Ni>>1; data1=allocate_1d_float((ofs+2*ext),0); data2=allocate_1d_float((ofs+2*ext),0); data1+=ext;data2+=ext; // Need third array since cannot access columns directly via pointers. data3=allocate_1d_float(Ni,0); for(j=0;j<Nj;j++) {   for(i=0;i<ofs;i++) {   k=i+ofs;   // low pass and high pass coefficients.   data1=image[j];   data2=image[k][j];   data3=data3[k]=0;   }   if(PS) {   extend_periodic(data1,ofs,ext);   extend_periodic(data2,ofs,ext);   }   else {   extend_mirror(data1,ofs,ext,shift_arr[lev]);   extend_mirror(data2,ofs,ext,1-shift_arr[lev]);   }   upsample_n_filter(data1,data3,Ni,lp,Nl,begilp,2);   upsample_n_filter(data2,data3,Ni,hp,Nh,begihp,2);   for(i=0;i<Ni;i++)   image[j]=data3; } data1-=ext;data2-=ext; free((void *)data1); free((void *)data2); free((void *)data3); } [right][snapback]1008971[/snapback][/right] </div> 我其实只需要Wavelet变换的代码。以上的是不是太多了？ 
作者: loveless 时间: 17.6.2006 17:09
那我该怎么办啊？我把东西都给你，你给我看看行吗？
作者: loveless 时间: 17.6.2006 17:31
好的，我再看下。
作者: loveless 时间: 17.6.2006 22:22
// Onur G. Guleryuz 1995, 1996, 1997, // University of Illinois at Urbana-Champaign, // Princeton University, // Polytechnic University. #include "wav_basic.h" #include "wav_filters_extern.h" #include "alloc.h" #include "wav_gen.h" // Macro that sets filter parameters. // Called via choose_filter(). #define set_filter_param(fname) \ { \ MFLP= fname ## _FLP; \ MFHP= fname ## _FHP; \ MILP= fname ## _ILP; \ MIHP= fname ## _IHP; \ \ /* Filter size. */ \ Nflp= fname ## _Nflp; \ Nfhp= fname ## _Nfhp; \ Nilp= fname ## _Nilp; \ Nihp= fname ## _Nihp; \ \ /* Filter symmetry */ \ PS= fname ## _PS; \ \ /* Filter starting points. */ \ begflp=mfl= fname ## _begflp; \ begfhp=mfh= fname ## _begfhp; \ begilp=mil= fname ## _begilp; \ begihp=mih= fname ## _begihp; \ } // Set wavelet filter parameters depending on name // and tap. See wav_filters.h and wav_filters_extern.h // This routine must be called before doing any transforms // to select which wavelet bank to use. void choose_filter(char name,int tap) { switch (name) { case 'A': switch (tap) { default: set_filter_param(Antonini); } break; case 'H': switch (tap) { default: set_filter_param(Haar); } break; case 'N': switch (tap) { case -1: set_filter_param(Nick_flip); break; default: set_filter_param(Nick); } break; default: switch (tap) { case 9: set_filter_param(D79); break; default: set_filter_param(D8); } } } // 2D "in place" wavelet transform. Input data is in // float **image and output is returned inside image. // // levs>=1 is the number of wavelet levels. // // Ni and Nj are the dimensions of the image. // Ni and Nj should be divisible by 2^levs. // // shift_arr_row and shift_arr_col (each of dimension levs) // contain the amount of shift (0 or 1) // that should be applied to filters in each level for transforms over rows and columns. // The use of shift_arr_row and shift_arr_col is mainly intended for // overcomplete filtering/transforms. // // forw != 0 for a forward transform (inverse is calculated otherwise). // // float *lp of dimension Nl contains the low pass wavelet filter. // float *hp of dimension Nh contains the high pass wavelet filter. // void wav2d_inpl(float **image,int Ni,int Nj,int levs,float *lp,int Nl, float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col) { int i,dimi,dimj; dimi=(forw)?NiNi>>(levs-1)); dimj=(forw)?NjNj>>(levs-1)); for(i=0;i<levs;i++) { if(forw) { // Do the rows. // Adjust filter parameters for overcomplete filtering. begflp=mfl+shift_arr_row; begfhp=mfh-shift_arr_row; filt_n_dec_all_rows(image,dimi,dimj,lp,Nl,hp,Nh); // Now the columns. begflp=mfl+shift_arr_col; begfhp=mfh-shift_arr_col; filt_n_dec_all_cols(image,dimi,dimj,lp,Nl,hp,Nh); } else { // Rows. begilp=mil+shift_arr_row[levs-1-i]; begihp=mih-shift_arr_row[levs-1-i]; ups_n_filt_all_rows(image,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row); // Columns. begilp=mil+shift_arr_col[levs-1-i]; begihp=mih-shift_arr_col[levs-1-i]; ups_n_filt_all_cols(image,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col); } dimi=(forw)?(dimi>>1)dimi<<1); dimj=(forw)?(dimj>>1)dimj<<1); } } // 2-D "in place" wavelet packet transform. // See wav2d_inpl() for an explanation of parameters. // // This routine will grow a full packet tree for levels upto LL_ONLY_AFTER_LEV // and just a wavelet tree afterward, i.e., if levs>LL_ONLY_AFTER_LEV // we grow wavelets only over LL band after LL_ONLY_AFTER_LEV. // LL_ONLY_AFTER_LEV is set in wav_gen.h Just set it to some large number // if you want packets all the way. // // Ni and Nj should be divisible by 2^levs. // shift_arr_row and shift_arr_col contain the amount of shift (0 or 1) // that should be applied to filters in each level. // This is mainly intended for overcomplete filtering. // See wav2d_inpl() comments. // void wavpack2d_inpl(float **image,int Ni,int Nj,int levs,float *lp,int Nl, float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col) { int i,k,l,dimi,dimj; int p,q,rcnt,ccnt; float **buffer; buffer=allocate_2d_float(Ni,Nj,0); dimi=(forw)?NiNi>>(levs-1)); dimj=(forw)?Nj:(Nj>>(levs-1)); for(i=0;i<levs;i++) { // Number of row bands. rcnt=Ni/dimi; // Number of column bands. ccnt=Nj/dimj; if(forw&&(i>=LL_ONLY_AFTER_LEV)) { // Grow wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } else if((!forw)&&(levs-1-i>=LL_ONLY_AFTER_LEV)) { // Inverse wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } // Do the next level of the wavelet tree over all of // the bands we are supposed to grow on. for(p=0;p<rcnt;p++) for(q=0;q<ccnt;q++) { for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) buffer[k][l]=image[p*dimi+k][q*dimj+l]; if(forw) { // Do the rows. // Adjust filter parameters for overcomplete filtering. begflp=mfl+shift_arr_row; begfhp=mfh-shift_arr_row; filt_n_dec_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh); // Now the columns. begflp=mfl+shift_arr_col; begfhp=mfh-shift_arr_col; filt_n_dec_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh); } else { // Rows. begilp=mil+shift_arr_row[levs-1-i]; begihp=mih-shift_arr_row[levs-1-i]; ups_n_filt_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row); // Columns. begilp=mil+shift_arr_col[levs-1-i]; begihp=mih-shift_arr_col[levs-1-i]; ups_n_filt_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col); } for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) image[p*dimi+k][q*dimj+l]=buffer[k][l]; } dimi=(forw)?(dimi>>1):(dimi<<1); dimj=(forw)?(dimj>>1):(dimj<<1); } free_2d_float(buffer,Ni); } // Forward complex wavelet transform using Kingsbury's paper. // Please check with Nick Kingsbury's paper to confirm that // I am doing things correctly here. // // Input is in float **im, and the results are returned in 4 // **float's, trf[0],...,trf[3]. // // Ni, Nj multiples of 2^levs. // levs>=1. // See wav2d_inpl() for an explanation of parameters. // void complex_wav_forw(float **im,float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp; // Choose Antonini. choose_filter('A',0); lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // First level, fully overcomplete resulting 4x redundancy. // Results are in trf. for(i=0;i<4;i++) { for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) trf[k][l]=im[k][l]; // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c); } // Now grow complex wavelets over the LL band // using Kingsbury's orthogonal bank. if(levs>1) { for(i=0;i<4;i++) { dimi=Ni/2; dimj=Nj/2; for(j=1;j<levs;j++) { if(i/2==0) { // Regular orth. bank. choose_filter('N',0); } else { // Flipped orth. bank (see Kingsbury). // Flipped is same as inverse. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over rows. filt_n_dec_all_rows(trf,dimi,dimj,lp,Nl,hp,Nh); if(i%2==0) { // Regular. choose_filter('N',0); } else { // Flipped. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over columns. filt_n_dec_all_cols(trf,dimi,dimj,lp,Nl,hp,Nh); dimi>>=1; dimj>>=1; } } } // Generate the complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // After this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=trf[0][k][l]+trf[3][k][l]; trf[3][k][l]=trf[0][k][l]-trf[3][k][l]; trf[0][k][l]=tmp; tmp=trf[1][k][l]+trf[2][k][l]; // It is possible that the below is 2 - 1. trf[2][k][l]=trf[1][k][l]-trf[2][k][l]; trf[1][k][l]=tmp; } } // Inverse complex wavelet transform. // trf is modified in intermediate computations. // Ni, Nj multiples of 2^levs // levs>=1. // // See complex_wav_forw() comments. // float **complex_wav_inv(float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp; float **im; // Clean the filter shift arrays. for(i=0;i<1<<levs;i++) shift_arr_r=shift_arr_c=0; // Get back the intermediate wavelet coefficients from the // complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // Prior to this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=(trf[0][k][l]+trf[3][k][l])/2; trf[3][k][l]=(trf[0][k][l]-trf[3][k][l])/2; trf[0][k][l]=tmp; tmp=(trf[1][k][l]+trf[2][k][l])/2; trf[2][k][l]=(trf[1][k][l]-trf[2][k][l])/2; trf[1][k][l]=tmp; } if(levs>1) { for(i=0;i<4;i++) { dimi=Ni>>(levs-1); dimj=Nj>>(levs-1); for(j=1;j<levs;j++) { if(i/2==0) { // Regular inverse, it just happens that orth. inverses are flipped. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over rows. ups_n_filt_all_rows(trf,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_r); if(i%2==0) { // Regular inverse. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over columns. ups_n_filt_all_cols(trf,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_c); dimi<<=1; dimj<<=1; } } } // Initialize final result with zeros. im=allocate_2d_float(Ni,Nj,1); // Choose Antonini. choose_filter('A',0); lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Now the first level. for(i=0;i<4;i++) { // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Inverse transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c); // Accumulate the results of 4x redundancy. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]+=trf[k][l]; } // Normalize results by 4 to account for 4x redundancy accumulation. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]/=4; return(im); } // Forward complex wavelet packet transform. // Ni, Nj multiples of 2^levs. // levs>=1. // // See complex_wav_forw() comments. // See wavpack2d_inpl() comments for LL_ONLY_AFTER_LEV. // // Not fully tested. Please check with literature. // void complex_wav_pack_forw(float **im,float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp,**buffer; int p,q,rcnt,ccnt; // Choose Antonini. choose_filter('A',0); lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // First level, fully overcomplete resulting 4x redundancy. // Results are in trf. for(i=0;i<4;i++) { for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) trf[k][l]=im[k][l]; // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c); } // Now grow using Kingsbury's orthogonal bank. if(levs>1) { buffer=allocate_2d_float(Ni/2,Nj/2,0); for(i=0;i<4;i++) { dimi=Ni/2; dimj=Nj/2; for(j=1;j<levs;j++) { // Number of row/column bands. rcnt=Ni/dimi; ccnt=Nj/dimj; if(j>=LL_ONLY_AFTER_LEV) { // Grow wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } for(p=0;p<rcnt;p++) for(q=0;q<ccnt;q++) { // Copy data to buffer. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) buffer[k][l]=trf[p*dimi+k][q*dimj+l]; if(i/2==0) { // Regular orth. bank. choose_filter('N',0); } else { // Flipped orth. bank (see Kingsbury). // Flipped is same as inverse. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over rows. filt_n_dec_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh); if(i%2==0) { // Regular. choose_filter('N',0); } else { // Flipped. choose_filter('N',-1); } lp=MFLP;Nl=Nflp; hp=MFHP;Nh=Nfhp; // Transform over columns. filt_n_dec_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh); // Copy results to trf. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) trf[p*dimi+k][q*dimj+l]=buffer[k][l]; } dimi>>=1; dimj>>=1; } } free_2d_float(buffer,Ni/2); } // Generate the complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // After this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=trf[0][k][l]+trf[3][k][l]; trf[3][k][l]=trf[0][k][l]-trf[3][k][l]; trf[0][k][l]=tmp; tmp=trf[1][k][l]+trf[2][k][l]; trf[2][k][l]=trf[1][k][l]-trf[2][k][l]; trf[1][k][l]=tmp; } } // Inverse complex wavelet packet transform. // trf is modified in intermediate computations. // Ni, Nj multiples of 2^levs // levs>=1. // // See complex_wav_pack_forw() comments. // float **complex_wav_pack_inv(float ***trf,int Ni,int Nj,int levs) { int i,j,k,l; int Nl,Nh,dimi,dimj,hdimi,hdimj; float *lp,*hp; int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE]; float tmp; float **im,**buffer; int p,q,rcnt,ccnt; // Clean the filter shift arrays. for(i=0;i<1<<levs;i++) shift_arr_r=shift_arr_c=0; // Get back the intermediate wavelet coefficients from the // complex wavelet coefficients. hdimi=Ni/(1<<levs); hdimj=Nj/(1<<levs); for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) { // Prior to this computation the complex wavelet coefficients // are given by trf[0]+jtrf[2] and trf[3]+jtrf[1]; tmp=(trf[0][k][l]+trf[3][k][l])/2; trf[3][k][l]=(trf[0][k][l]-trf[3][k][l])/2; trf[0][k][l]=tmp; tmp=(trf[1][k][l]+trf[2][k][l])/2; trf[2][k][l]=(trf[1][k][l]-trf[2][k][l])/2; trf[1][k][l]=tmp; } if(levs>1) { buffer=allocate_2d_float(Ni/2,Nj/2,0); for(i=0;i<4;i++) { dimi=Ni>>(levs-1); dimj=Nj>>(levs-1); for(j=1;j<levs;j++) { // Number of row/column bands. rcnt=Ni/dimi; ccnt=Nj/dimj; if(j>=LL_ONLY_AFTER_LEV) { // Grow wavelets only over LL band after LL_ONLY_AFTER_LEV. rcnt=ccnt=1; } for(p=0;p<rcnt;p++) for(q=0;q<ccnt;q++) { // Copy data to buffer. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) buffer[k][l]=trf[p*dimi+k][q*dimj+l]; if(i/2==0) { // Regular inverse, it just happens that orth. inverses are flipped. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over rows. ups_n_filt_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_r); if(i%2==0) { // Regular inverse. choose_filter('N',0); } else { // Flipped inverse. choose_filter('N',-1); } lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Inverse transform over columns. ups_n_filt_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-j,shift_arr_c); // Copy results to trf. for(k=0;k<dimi;k++) for(l=0;l<dimj;l++) trf[p*dimi+k][q*dimj+l]=buffer[k][l]; } dimi<<=1; dimj<<=1; } } free_2d_float(buffer,Ni/2); } // Initialize final result with zeros. im=allocate_2d_float(Ni,Nj,1); // Choose Antonini. choose_filter('A',0); lp=MILP;Nl=Nilp; hp=MIHP;Nh=Nihp; // Now the first level. for(i=0;i<4;i++) { // Set shifts as specified by Kingsbury's paper. if(i/2==0) shift_arr_r[0]=1; else shift_arr_r[0]=0; if(i%2==0) shift_arr_c[0]=1; else shift_arr_c[0]=0; // Inverse transform. wav2d_inpl(trf,Ni,Nj,1,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c); // Accumulate the results of 4x redundancy. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]+=trf[k][l]; } // Normalize results by 4 to account for 4x redundancy accumulation. for(k=0;k<Ni;k++) for(l=0;l<Nj;l++) im[k][l]/=4; return(im); }
作者: loveless 时间: 17.6.2006 22:25
我实在看得不行了，这个是不是正宗的Wavelet Transformation? 那个很多的，是不是很多多余的？

欢迎光临德国开元华人社区开元周游 (https://bbs.kaiyuan.cn/) Powered by Discuz! X3.2

德国开元华人社区 开元周游

德国开元华人社区开元周游