[编程练习] 卷积和实现demo——对阶跃函数进行一维高斯模糊

#include 
   
   
    
    
#include 
    
    
     
     
#include 
     
     
      
      
#include 
      
      
       
       
#define _PI_ 3.14159265358979323846264
#define _E_  2.7182818284590452353602874713527

template 
       
       
         struct Array{ T * data; int length; Array(){data=0;length=0;}; Array(int len){ printf("allocing %d\n",len); data = new T[len]; length = len; }; Array(Array 
        
          & a){data = new T[a->length];length = a->length;}; }; template 
         
           struct pair{ T a,b; }; double normal(double x, double u, double sigma){ return exp((-1*pow(x-u,2))/(2*sigma*sigma))/(sqrtf(2*_PI_)*sigma); } Array< pair 
          
            > * normalCurve(double u, double sigma, int count, float width = 3.0){ double step = sigma * width * 2.0 /(double)count; double ptr = u-sigma*width; Array< pair 
           
             > * curve = new Array< pair 
            
              >(count); for(int i=0;i 
             
               data[i].a = ptr; curve->data[i].b = normal(ptr,u,sigma); ptr += step; } puts("normalCurve"); return curve; } Array 
              
                * curveTemplate(Array< pair 
               
                 > * a){ Array 
                
                  * tpl = new Array 
                 
                   (a->length); for(int i=0;i 
                  
                    length;i++) tpl->data[i] = a->data[i].b; puts("curveTemplate"); return tpl; } void curvePrt( Array< pair 
                   
                     > * curve, FILE * f = stdout){ for(int i=0;i 
                    
                      length;i++) fprintf(f,"%lf,%lf\n",curve->data[i].a,curve->data[i].b); puts("curvePrt"); return; } void curvePrt( Array< pair 
                     
                       > * curve, const char * filename){ FILE * f = fopen(filename,"w"); for(int i=0;i 
                      
                        length;i++) fprintf(f,"%lf,%lf\n",curve->data[i].a,curve->data[i].b); puts("curvePrt"); fclose(f); return; } Array 
                       
                         * convolution( Array 
                        
                          * a, Array 
                         
                           * b ){ int bstart; int bend; int blenh = b->length / 2; Array 
                          
                            * conv = new Array 
                           
                             (a->length); memset(conv->data,0,a->length*sizeof(double)); for(int i=0;i 
                            
                              length;i++){ bstart = (i-blenh)<0 ? (blenh-i) : 0; bend = (i+blenh)>a->length ? (a->length-i+blenh) : b->length; for(int j=bstart;j 
                             
                               data[i+j-blenh] += a->data[i]*b->data[j]; } } puts("convolution"); return conv; } Array< pair 
                              
                                > * makeCurve( Array 
                               
                                 * a, double start, double step ){ Array< pair 
                                
                                  > * curve = new Array< pair 
                                 
                                   >(a->length); double prt = start; for(int i=0;i 
                                  
                                    length;i++){ curve->data[i].a = prt; curve->data[i].b = a->data[i]; prt += step; } puts("makeCurve"); return curve; } int main(){ int stepSize = 2000; Array 
                                   
                                     * stepFunc = new Array 
                                    
                                      (stepSize); memset(stepFunc->data,0,stepSize*sizeof(double)); for(int i=stepSize/2-40;i 
                                     
                                       data[i] = 1; } puts("start"); curvePrt( normalCurve(0,0.4,100) , "norm.csv" ); curvePrt( makeCurve( stepFunc,-1*stepSize/2/10.0,1/10.0 ) , "step.csv" ); curvePrt( makeCurve( convolution( stepFunc, curveTemplate( normalCurve( 0, 0.4, 100 )) ), -1*stepSize/2/10.0, 1/10.0 ) , "stdout.csv" ); }