Python实现循环卷积算法

本文介绍循环卷积算法的基本原理，并使用Python实现循环卷积算法。

1 基本原理

循环卷积的定义式是
$$y(k)=\sum _{n=0}^{N-1}x(n)h(k-n);k=1,\mathrm{2...},N-1$$
如果k-n<0，则h的序号加N。对于长度不足N的序列，在其末尾补0，上式展开为矩阵的形式为
$$\left[\begin{array}{c}y(0)\\ y(1)\\ y(2)\\ y(3)\end{array}\right]=\left[\begin{array}{cccc}h(0)& h(3)& h(2)& h(1)\\ h(1)& h(0)& h(3)& h(2)\\ h(2)& h(1)& h(0)& h(3)\\ h(3)& h(2)& h(1)& h(0)\end{array}\right]\left[\begin{array}{c}x(0)\\ x(1)\\ x(2)\\ x(3)\end{array}\right]$$
根据卷积定理，函数卷积的傅里叶变换是函数傅里叶变换的乘积，所以循环卷积可以转换为先算序列x和h的傅里叶变换，然后将傅里叶变换的结果相乘，最后对乘积后的序列进行傅里叶逆变换，如下式所示
y ( k ) = I F F T { F F T [ x ( n ) ] × F F T [ h ( n ) ] } y(k) = IFFT\{FFT[x(n)]×FFT[h(n)]\} y(k)=IFFT{FFT[x(n)]×FFT[h(n)]}

2 程序实现

import numpy as np
import cmath
from numpy import pi

def calc1(x, h, y): #直接使用卷积定理计算
  length = len(x)
  for i in range(length):
    for j in range(length):
      idx = i-j
      if (idx < 0):
        idx = idx+length
      y[i] = y[i] + x[j]*h[idx]

def calc2(x, h, y): #使用快速傅里叶变换计算
  x_fft = np.fft.fft(x)
  h_fft = np.fft.fft(h)
  length = len(x)
  xh = np.zeros(length, complex)
  for i in range(length):
    xh[i] = x_fft[i] * h_fft[i]
  xh_ifft = np.zeros(length, complex)
  xh_ifft = np.fft.ifft(xh)
  #print(xh_ifft)
  for i in range(length):
    y[i] = complex(xh_ifft[i].real,xh_ifft[i].imag)

########################################
N = 5
x = np.zeros(N, complex)
h = np.zeros(N, complex)
for i in range(N):
  x[i] = i*100 + i*150j
  h[i] = cmath.exp(1j*pi*(i-N)**2/N)
length = len(x)
y = np.zeros(length, complex)
calc1(x, h, y)
print('y = ', end='')
print(y)
calc2(x, h, y)
print('y = ', end='')
print(y)

3 运行结果