Understanding in convolutions

gistfile1.py

import numpy as np

# http://www.songho.ca/dsp/convolution/convolution.html#separable_convolution
# https://www.mathworks.com/help/signal/ug/linear-and-circular-convolution.html



# make circular convolution equivalent to linear convolution
# M = 13
# N = 7
# a = np.random.random(M)
# b = np.random.random(N)
#
# t1 = np.zeros(M+N-1)
# for i in xrange(M+N-1):
#     s = 0
#     for j in xrange(i+1):
#         if j < M and i-j < N:
#             s += a[j] * b[i-j]
#     t1[i] = s
#
#
# t2 = np.fft.irfft(np.fft.rfft(a, n=M+N-1) * np.fft.rfft(b, n=M+N-1), n=M+N-1)
#
# print "t1=", t1
# print "t2=", t2
#
# print t1 - t2
# print np.all(t1 == t2)




# check 2d circular convolution equivalent to 2d linear convolution
# M = 5
# N = 5
# k = 3
# a = np.random.random([M, N])
# b = np.random.random([k, k])
#
# s = [M+k-1, M+k-1]
# t1 = np.zeros(s)
# for m in xrange(M+k-1):
#     for n in xrange(M+k-1):
#         for i in xrange(m+1):
#             for j in xrange(n+1):
#                 if i < M and j < N and m-i < k and n-j < k:
#                     t1[m, n] += a[i, j] * b[m-i, n-j]
#
#
#
# t2 = np.fft.irfft2(np.fft.rfft2(a, s) * np.fft.rfft2(b, s), s)
# print "t1=", t1
# print "t2=", t2
#
# print t1 - t2
# print np.all(t1 == t2)



# check column vector times row vector is equivalent to their linear convolution
# M = 5
# N = 5
# a = np.random.random([M, 1])
# b = np.random.random([1, N])
#
# s = [M+1-1, M+1-1]
# t1 = np.zeros(s)
# for m in xrange(M+1-1):
#     for n in xrange(1+N-1):
#         for i in xrange(m+1):
#             for j in xrange(n+1):
#                 if i < M and j < 1 and m-i < 1 and n-j < N:
#                     t1[m, n] += a[i, j] * b[m-i, n-j]
#
# t2 = np.fft.irfft2(np.fft.rfft2(a, s) * np.fft.rfft2(b, s), s)
#
# t3 = np.dot(a, b)
#
# print "t1=", t1
# print "t2=", t2
# print "t3=", t3
#
#
# print t1 - t2
# print np.all(t1 == t2)
#
# print t1 - t3
# print np.all(t1 == t3)






# check that linear convolution is associative
def conv2d(a, b):
    M, N = a.shape
    P, Q = b.shape
    s = [M + P - 1, N + Q - 1]
    t = np.zeros(s)
    for m in xrange(M +P - 1):
        for n in xrange(N + Q - 1):
            for i in xrange(m + 1):
                for j in xrange(n + 1):
                    if i < M and j < N and m - i < P and n - j < Q:
                        t[m, n] += a[i, j] * b[m - i, n - j]
    return t


def conv2d_fft(a, b):
    M, N = a.shape
    P, Q = b.shape
    s = [M + P - 1, N + Q - 1]
    t = np.fft.irfft2(np.fft.rfft2(a, s) * np.fft.rfft2(b, s), s)
    return np.real(t)
#
#
#
# a = np.random.random([37, 37])
# b = np.random.random([7, 1])
# c = np.random.random([1, 5])
#
#
#
#
#
#
# ft1 = conv2d_fft(conv2d_fft(a, b), c)
# ft2 = conv2d_fft(a, conv2d_fft(b, c))
# ft3 = conv2d_fft(a, np.dot(b, c))
#
# #print ft1 - ft2
# print np.all(ft1 - ft2 < 1e-12)
#
# #print ft1 - ft3
# print np.all(ft1 - ft3 < 1e-12)
#
#
#
#
#
# t1 = conv2d(conv2d(a, b), c)
# t2 = conv2d(a, conv2d(b, c))
# t3 = conv2d(a, np.dot(b, c))
#
#
# # print "t1=", t1
# # print "t2=", t2
# # print "t3=", t3
#
# #print t1 - t2
# print np.all(t1 - t2 < 1e-14)
#
# ##print t1 - t3
# print np.all(t1 - t3 < 1e-14)
#
#
# print np.all(ft1 - t1 < 1e-12)




# replace 2d convolution with 1d convolution
def conv1d(a, b, M, N): # M for a, N for b
    t = np.zeros(M+N-1)
    for i in xrange(M+N-1):
        s = 0
        for j in xrange(i+1):
            if j < M and i-j < N:
                s += a[j] * b[i-j]
        t[i] = s
    return t


def conv1d_fft(a, b, M, N): # M for a, N for b
    t = np.fft.irfft(np.fft.rfft(a, n=M+N-1) * np.fft.rfft(b, n=M+N-1), n=M+N-1)
    return t

M = 3
N = 5
a = np.random.random([31, 31])
b = np.random.random(M) # [M, 1]
c = np.random.random(N) # [1, N]

# v = np.apply_along_axis(conv1d, 0, a, b, a.shape[1], M) # slice according to dim 0
# #print "v.shape", v.shape
# h = np.apply_along_axis(conv1d, 1, v, c, a.shape[0], N) # slice according to dim 0
# #print "h.shape", h.shape
#
# t = conv2d(a, np.dot(b.reshape([M, 1]), c.reshape([1, N])))
#
# print np.all(h - t < 1e-12)


# v = np.apply_along_axis(conv1d_fft, 0, a, b, a.shape[1], M) # slice according to dim 0
# h = np.apply_along_axis(conv1d_fft, 1, v, c, a.shape[0], N) # slice according to dim 0
# t = conv2d(a, np.dot(b.reshape([M, 1]), c.reshape([1, N])))
# print np.all(h - t < 1e-12)




v = np.apply_along_axis(conv1d_fft, 0, a, b, a.shape[1], M) # slice according to dim 0
h = np.apply_along_axis(conv1d_fft, 1, v, c, a.shape[0], N) # slice according to dim 0
t = conv2d_fft(a, np.dot(b.reshape([M, 1]), c.reshape([1, N])))
print np.all(h - t < 1e-12)