Python实现循环卷积算法
- 1 基本原理
- 2 程序实现
- 3 运行结果
本文介绍循环卷积算法的基本原理,并使用Python实现循环卷积算法。
1 基本原理
循环卷积的定义式是
y ( k ) = ∑ n = 0 N − 1 x ( n ) h ( k − n ) ; k = 1 , 2... , N − 1 y(k) = \sum_{n=0}^{N-1}x(n)h(k-n); k=1,2...,N-1 y(k)=n=0∑N−1x(n)h(k−n);k=1,2...,N−1
如果k-n<0,则h的序号加N。对于长度不足N的序列,在其末尾补0,上式展开为矩阵的形式为
[ y ( 0 ) y ( 1 ) y ( 2 ) y ( 3 ) ] = [ 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 ) ] [ x ( 0 ) x ( 1 ) x ( 2 ) x ( 3 ) ] \begin{bmatrix} y(0)\\y(1)\\y(2)\\y(3) \end{bmatrix}= \begin{bmatrix} 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{bmatrix}\begin{bmatrix} x(0)\\x(1)\\x(2)\\x(3) \end{bmatrix} y(0)y(1)y(2)y(3) = h(0)h(1)h(2)h(3)h(3)h(0)h(1)h(2)h(2)h(3)h(0)h(1)h(1)h(2)h(3)h(0) x(0)x(1)x(2)x(3)
根据卷积定理,函数卷积的傅里叶变换是函数傅里叶变换的乘积,所以循环卷积可以转换为先算序列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 pidef calc1(x, h, y): #直接使用卷积定理计算length = len(x)for i in range(length):for j in range(length):idx = i-jif (idx < 0):idx = idx+lengthy[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*150jh[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)