Bonjour à tous,

J'aimerai écrire un programme qui effectue le filtrage d'une fonction par un noyau gaussien.
Pour se faire j'ai codé une fonction que me permet de faire une convolution de ma fonction (signal) avec le noyau gaussien. Mon souci c'est que j'ai des problèmes aux bords. Pour y palier j'ai considéré que ma fonction est périodique, mais les résultats obtenus ne sont pas ceux escomptés. Pouvez vous me dire ce qui ne va pas dans mon programme (il est écrit en C)?


/* noyau gaussien */
int noyauGaussien1D(double **noyauGaussien, double sigma)
{

int rayon=(int)round(sigma*3.); /*rayon*2+1 est la taille de la fenetre*/
short i;
double *noyau;
double constante;

noyau=malloc(((rayon*2)+1)*sizeof(double));

constante=1/(sigma*sqrt(2.*M_PI));

for(i=-rayon;i<=rayon;++i)
                       noyau[rayon+i]=constante*exp(-(i*i)/(2.*sigma*sigma));
  
                                                                        
*noyauGaussien=noyau;
            
return(rayon);

}

/*application du filtre a  l'aide de la convolution*/
void filtre_gaussien1D(int n, double *signal,double sigma,double *out)
{
  short i,j,k;
  int rayon;
  double *noyau;
  double *signalPeriodique;
  double temp;
  
  /* ajout d'elements dans le tableau d'entree pour eviter les effets de bords (on considere que le signal est periodique) */
  rayon=noyauGaussien1D(&noyau,sigma);
  signalPeriodique=malloc((n+(2.*rayon))*sizeof(double));
  
  for(i=0; i<rayon; ++i){
   signalPeriodique[i]=signal[n-rayon+i];
   signalPeriodique[n+rayon+i]=signal[i];
                       }
                      
  for(i=0; i<n; ++i) signalPeriodique[i+rayon]=signal[i];
      
    
   /* application de la convolution*/              
   for(i = rayon,k = 0; (i < n+rayon) & (k < n); ++i,++k){   temp=0.;
       for(j=-rayon; j<=rayon; ++j)
                             temp+=(signalPeriodique[i+j]*noyau[rayon+j]);
                                
                                            out[k]=temp;
                        }
                        
      free(noyau);
      free(signalPeriodique);                  
  return;
  
}