title: High Pass Filter using FFTW in C
body: <p>I have a question regarding FFT. I already manage to do FFT forward and backward using FFTW in C. Now, I want to apply high pass filter for edge detection, some of my source said that just zeroing the centre of the magnitude.</p>
  
  <p>This is my input image
  <a href="http://i62.tinypic.com/2wnxvfl.jpg" rel="nofollow">http://i62.tinypic.com/2wnxvfl.jpg</a></p>
  
  <p>Basically what I do are :</p>
  
  <ol>
  <li>Forward FFT</li>
  <li>Convert the output to 2D array</li>
  <li>Do forward FFT shifting</li>
  <li>Make the real and imag value to 0 when the distance from the centre is 25% of the height</li>
  <li>Generate the magnitude</li>
  <li>Do backward FFT shifting</li>
  <li>Convert into 1D array</li>
  <li>Do Backward FFT.</li>
  </ol>
  
  <p>This is the original magnitude, the processed magnitude, and the result</p>
  
  <p><a href="http://i58.tinypic.com/aysx9s.png" rel="nofollow">http://i58.tinypic.com/aysx9s.png</a></p>
  
  <p>can someone help me, to tell me which part is wrong and how to do the high pass filtering using FFTW in C.</p>
  
  <p>Thank You.</p>
  
  <p>The Source Code:</p>
  
  <pre><code>unsigned char **FFT2(int width,int height, unsigned char **pixel, char line1[100],char line2[100], char line3[100],char filename[100])
  {
    fftw_complex* in, * dft, * idft, * dft2;
  
    //fftw_complex tmp1,tmp2;
    fftw_plan plan_f,plan_i;
    int i,j,k,w,h,N,w2,h2;
  
    w = width;
    h = height;
    N = w*h;
  
    unsigned char **pixel_out;
    pixel_out = malloc(h*sizeof(unsigned char*));
    for(i = 0 ; i&lt;h;i++)
      pixel_out[i]=malloc(w*sizeof(unsigned char));
  
  
  
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
    dft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
    dft2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
    idft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
  
  
  
    /*run forward FFT*/
  
    plan_f = fftw_plan_dft_2d(w,h,in,dft,FFTW_FORWARD,FFTW_ESTIMATE);
    for(i = 0,k = 0 ; i &lt; h ; i++)
    {
      for(j = 0 ; j &lt; w ; j++,k++)
      {
        in[k][0] = pixel[i][j];
        in[k][1] = 0.0;
      }
    }
    fftw_execute(plan_f);
  
  
    double maxReal = 0.0;
    for(i = 0 ; i &lt; N ; i++)
      maxReal = dft[i][0] &gt; maxReal ? dft[i][0] : maxReal;
  
    printf("MAX REAL : %f\n",maxReal);
  
  
  
  
    /*fftshift*/
    //convert to 2d
    double ***temp1;
    temp1 = malloc(h * sizeof (double**));
    for (i = 0;i &lt; h; i++){
        temp1[i] = malloc(w*sizeof (double*));
        for (j = 0; j &lt; w; j++){
            temp1[i][j] = malloc(2*sizeof(double));
        }
    }
  
    double ***temp2;
    temp2 = malloc(h * sizeof (double**));
    for (i = 0;i &lt; h; i++){
        temp2[i] = malloc(w*sizeof (double*));
        for (j = 0; j &lt; w; j++){
            temp2[i][j] = malloc(2*sizeof(double));
        }
    }
  
  
  
  
    for (i = 0;i &lt; h; i++){
        for (j = 0; j &lt; w; j++){
            temp1[i][j][0] = dft[i*w+j][0];
            temp1[i][j][1] = dft[i*w+j][1];
        }
    }
  
  
  
  
  
    int m2 = h/2;
    int n2 = w/2;
  
  
  
    //forward shifting
    for (i = 0; i &lt; m2; i++)
    {
       for (k = 0; k &lt; n2; k++)
       {
            double tmp13[2]         = {temp1[i][k][0],temp1[i][k][1]};
            temp1[i][k][0]       = temp1[i+m2][k+n2][0];
            temp1[i][k][1]    = temp1[i+m2][k+n2][1];
            temp1[i+m2][k+n2][0] = tmp13[0];
            temp1[i+m2][k+n2][1] = tmp13[1];
  
            double tmp24[2]       = {temp1[i+m2][k][0],temp1[i+m2][k][1]};
            temp1[i+m2][k][0]    = temp1[i][k+n2][0];
            temp1[i+m2][k][1]    = temp1[i][k+n2][1];
            temp1[i][k+n2][0]   = tmp24[0];
            temp1[i][k+n2][1]    = tmp24[1];
       }
    }
  
  
  
    //process
  
    for (i = 0;i &lt; h; i++){
        for (j = 0; j &lt; w; j++){
            if(distance_to_center(i,j,m2,n2) &lt; 0.25*h)
            {
              temp1[i][j][0] = (double)0.0;
              temp1[i][j][1] = (double)0.0;
            }
        }
    }
  
  
  
    /* copy for magnitude */
    for (i = 0;i &lt; h; i++){
        for (j = 0; j &lt; w; j++){
            temp2[i][j][0] = temp1[i][j][0];
            temp2[i][j][1] = temp1[i][j][1];
        }
    }
  
  
    //backward shifting
    for (i = 0; i &lt; m2; i++)
    {
       for (k = 0; k &lt; n2; k++)
       {
            double tmp13[2]         = {temp1[i][k][0],temp1[i][k][1]};
            temp1[i][k][0]       = temp1[i+m2][k+n2][0];
            temp1[i][k][1]    = temp1[i+m2][k+n2][1];
            temp1[i+m2][k+n2][0] = tmp13[0];
            temp1[i+m2][k+n2][1] = tmp13[1];
  
            double tmp24[2]       = {temp1[i+m2][k][0],temp1[i+m2][k][1]};
            temp1[i+m2][k][0]    = temp1[i][k+n2][0];
            temp1[i+m2][k][1]    = temp1[i][k+n2][1];
            temp1[i][k+n2][0]   = tmp24[0];
            temp1[i][k+n2][1]    = tmp24[1];
       }
    }
  
  
  
    //convert back to 1d
    for (i = 0;i &lt; h; i++){
        for (j = 0; j &lt; w; j++){
            dft[i*w+j][0] = temp1[i][j][0];
  
            dft[i*w+j][1] = temp1[i][j][1];
  
  
            dft2[i*w+j][0] = temp2[i][j][0];
  
            dft2[i*w+j][1] = temp2[i][j][1];
  
        }
    }
  
  
  
  
    /* magnitude */
  
    double max = 0;
    double min = 0;
    double mag=0;
    for (i = 0, k = 1; i &lt; h; i++){
        for (j = 0; j &lt; w; j++, k++){
            mag = sqrt(pow(dft2[i*w+j][0],2) + pow(dft2[i*w+j][1],2));
            if (max &lt; mag)
                max = mag;
        }
    }
  
  
    double **magTemp;
    magTemp = malloc(h * sizeof (double*));
    for (i = 0;i &lt; h; i++){
        magTemp[i] = malloc(w*sizeof (double));
    }
  
    for(i = 0,k = 0 ; i &lt; h ; i++)
    {
      for(j = 0 ; j &lt; w ; j++,k++)
      {
  
        double mag = sqrt(pow(dft2[i*w+j][0],2) + pow(dft2[i*w+j][1],2));
        mag = 255*(mag/max);
        //magTemp[i][j] = 255-mag; //Putih
        magTemp[i][j] = mag; //Item
  
  
      }
    }
  
  
  
    /* brightening magnitude*/
  
    for(i = 0,k = 0 ; i &lt; h ; i++)
    {
      for(j = 0 ; j &lt; w ; j++,k++)
      {
        //double temp = magTemp[i][j];
        double temp = (double)(255/(log(1+255)))*log(1+magTemp[i][j]);
        pixel_out[i][j] = (unsigned char)temp;
  
      }
    }
  
    generateImage(width,height,pixel_out,line1,line2,line3,filename,"magnitude");
  
  
    /* backward fft */
    plan_i = fftw_plan_dft_2d(w,h,dft,idft,FFTW_BACKWARD,FFTW_ESTIMATE);
    fftw_execute(plan_i);
    for(i = 0,k = 0 ; i &lt; h ; i++)
    {
      for(j = 0 ; j &lt; w ; j++,k++)
      {
        double temp = idft[i*w+j][0]/N;
        pixel_out[i][j] = (unsigned char)temp; //+ pixel[i][j];
  
      }
    }
    generateImage(width,height,pixel_out,line1,line2,line3,filename,"backward");
  
    return pixel_out;
  }
  </code></pre>
  
  <p><strong>EDIT new source code</strong></p>
  
  <p>I add this part before the forward shifting, the result is as expected also.</p>
  
  <pre><code>//proses
    //create filter
    unsigned char **pixel_filter;
    pixel_filter = malloc(h*sizeof(unsigned char*));
    for(i = 0 ; i&lt;h;i++)
      pixel_filter[i]=malloc(w*sizeof(unsigned char));
    for (i = 0;i &lt; h; i++){
        for (j = 0; j &lt; w; j++){
            if(distance_to_center(i,j,m2,n2) &lt; 20)
            {
              pixel_filter[i][j] = 0;
            }
            else
            {
              pixel_filter[i][j] = 255;
            }
        }
    }
    generateImage(width,height,pixel_filter,line1,line2,line3,filename,"filter1");
    for (i = 0; i &lt; m2; i++)
    {
       for (k = 0; k &lt; n2; k++)
       {
            unsigned char tmp13         = pixel_filter[i][k];
            pixel_filter[i][k]       = pixel_filter[i+m2][k+n2];
            pixel_filter[i+m2][k+n2] = tmp13;
  
            unsigned char tmp24       = pixel_filter[i+m2][k];
            pixel_filter[i+m2][k]    = pixel_filter[i][k+n2];
            pixel_filter[i][k+n2]   = tmp24;
       }
    }
    generateImage(width,height,pixel_filter,line1,line2,line3,filename,"filter2");
    for (i = 0;i &lt; h; i++){
        for (j = 0; j &lt; w; j++){
            temp1[i][j][0] *= pixel_filter[i][j]; 
            temp1[i][j][1] *= pixel_filter[i][j];
        }
    }
  </code></pre>