0
$\begingroup$

I'm trying to manually implement a convolution using FFTs and it isn't working as expected. I know I'm probably missing some subtlety with padding, shifting, or conjugation, (all of which I've tried playing with) but can't find it. Here's what I've got:

import numpy as np                                                                                                                                                                
import scipy.signal as sig                                                                                                                                                        
import matplotlib.pyplot as plt                                                                                                                                                   

# Both of these work as expected, and give the same result                                                                                                                                                                                  
def myconvolve(f, g) :                                                                                                                                                         
    #return sig.fftconvolve(f, g, mode='same')                                                                                                                                  
    return sig.convolve2d(f, g, 'same', boundary='fill', fillvalue=0)                                                                                                          

# Here's my version, which doesn't work                                                                                                                                                                              
def myconvolve2(f, g) :                                                                                                                                                        
    p1, p2 = (np.r_[f.shape]-g.shape).astype(int)//2                                                                                                                           
    g = np.pad(g,((p1,p1),(p2,p2)),mode='edge')                                                                                                                                   
    F = np.fft.fft2(f)                                                                                                                                                        
    G = np.fft.fft2(g)                           
    FG = F*np.conj(G)                                                                                                                                
    res = np.real(np.fft.ifft2(FG))                                                                                                                                               
    return res  

# Generate some test functions                                                                                                                                                    
X,Y = np.mgrid[-1:1:20j,-1:1:20j]                                                                                                                                                 
g = np.exp( -X**2 -Y**2 )**2                                                                                                                                                      
np.random.seed(1)                                                                                                                                                                 
f = np.random.rand(100,70)                                                                                                                                                        
f[f<.998] = 0                                                                                                                                                                     

# Convolve                                                                                                                                                                        
fc = myconvolve2(f,g)                                                                                                                                                             

# Plot                                                                                                                                                                            
plt.figure(1,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(f.T)                                                                                                               
plt.figure(2,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(g.T)                                                                                                               
plt.figure(3,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(fc.T)   

Output: f g fg fg_bad

left to right: f, g, f*g (correct), f*g (my version)

$\endgroup$
  • $\begingroup$ only had a cursory glance but I see inside your function definition you never actually do anything with F1 and F2? $\endgroup$ – DamienBradley Jan 11 at 2:52
  • $\begingroup$ You're right. In trimming things down for SE, I accidentally removed that line! Fixed. $\endgroup$ – argentum2f Jan 12 at 16:56
  • $\begingroup$ Also, does dsp not support syntax highlighting like other SE sites do? $\endgroup$ – argentum2f Jan 13 at 14:54
1
$\begingroup$

I figured out my problem. The kernel needs to be shifted so the 'center' is on the corner of the image (which acts as the origin in an FFT). The built-in ifftshift function works just fine for this. (Note, there are some subtleties here depending on whether you have even or odd shapes or differences in shape that can result in one row or column shifts. I picked ifftshift in this case because it matches the result of the comparison function.) There's also a very faint but visible difference between my implementation and the library function (a faint streaking) that could just be numerical errors or something.

import numpy as np                                                                                                                                                                
import scipy.signal as sig                                                                                                                                                        
import matplotlib.pyplot as plt                                                                                                                                                   

def myconvolve(f, g) :                                                                                                                                                            
    return sig.fftconvolve(f, g, mode='same')                                                                                                                                     
    return sig.convolve2d(f, g, 'same', boundary='fill', fillvalue=0)                                                                                                             

def myconvolve2(f, g) :                                                                                                                                                           
    # Pad g to equal size of f. This assumes f is smaller in both dimensions                                                                                                      
    # and that the difference between f and g dimensions are even                                                                                                                 
    p1, p2 = (np.r_[f.shape]-g.shape).astype(int)//2                                                                                                                              
    gpad = np.pad(g,((p1,p1),(p2,p2)),mode='edge')                                                                                                                                

    # Shift g to 'center' on top left corner (the 'origin' in an fft)                                                                                                             
    gpad = np.fft.ifftshift(gpad)                                                                                                                                                  

    # Multiply spectra                                                                                                                                                            
    FG = np.fft.fft2(f) * np.conj(np.fft.fft2(gpad))                                                                                                                              
    return np.real(np.fft.ifft2(FG)), gpad                                                                                                                                        

# Generate some test functions                                                                                                                                                    
X,Y = np.mgrid[-1:1:21j,-1:1:21j]                                                                                                                                                 
g = np.exp( -X**2*2 -Y**2*2 )**2                                                                                                                                                  
np.random.seed(1)                                                                                                                                                                 
f = np.random.rand(101,71)                                                                                                                                                        
f[f<.998] = 0                                                                                                                                                                     

# Convolve                                                                                                                                                                        
fg, gpad = myconvolve2(f,g)                                                                                                                                                       
fg += f/4                                                                                                                                                                         

# Plot                                                                                                                                                                            
plt.figure(1,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(f.T)                                                                                                               
plt.figure(2,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(gpad.T)                                                                                                            
plt.figure(3,figsize = (10,15), dpi=80); plt.clf(); plt.imshow(fg.T)                                                                                                              

plt.imsave('f.png',f.T)                                                                                                                                                           
plt.imsave('g.png',g.T)                                                                                                                                                           
plt.imsave('fg_correct.png',fg.T) 

f gpad fg

Results left to right: f, gpad, fg

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.