使用overlap-add方法计算两个信号的卷积示例

假设信号数组如下,使用overlap-add方法计算线性卷积

a = [1, 3, -1, 2, 5, 3, -2, -4, -2, 1]
h = [1, 2, -1]

step1： 分组

根据h长度划分a
h长度为3，每3个一分，将a分为4组，如下：

[
[1, 3, -1],
[2, 5, 3],
[-2, -4, -2],
[1,0,0]
]

最后一个数组不够补零处理。

step2： 补零

h的长度为3，M=3
a被划分的长度为3， L=3
所以每个block补零后的长度应该为
$$N=L+M-1$$
N = 5

补零后

[
[1, 3, -1, 0, 0],    # block 0
[2, 5, 3, 0, 0],     # block 1
[-2, -4, -2, 0, 0],  # block 2
[1,0,0, 0, 0]        # block 3
]

step 3: 计算线性卷积

这里有个tricky小技巧，将h转换为用于计算卷积的矩阵：

将h补零到长度5， 即 [1, 2, -1, 0, 0]，然后每次向右移动一位。
即

[
[1, 2, -1, 0, 0], 
[0, 1, 2, -1, 0],
[ 0, 0, 1, 2, -1],
[-1, 0,0,1, 2],
[2,-1, 0,0,1]
]

转置一下，得到用于计算卷积的矩阵h_matrix

[
[  1  0  0  -1  2]
[  2  1  0  0  -1]
[ -1  2  1  0  0]
[ 0  -1  2  1  0]
[ 0  0  -1  2  1]
]

然后用次矩阵每次乘block
所以

y[0] =  h_matrix * block 0 = [ 1, 5, 4, -5, 1]
y[1] =  h_matrix * block 1 = [ 2, 9, 11, 1, -3]
y[2] =  h_matrix * block 2 = [ -2, -8, -8, 0, 2]
y[3] =  h_matrix * block 3 = [ 1, 2, -1, 0, 0]

step 4: Overlap-add计算

M-1部分做overlap-add计算

1, 5, 4, -5, 1
          2, 9, 11, 1, -3
                   -2, -8, -8, 0, 2
                               1, 2, -1, 0, 0
--------------------------------------------------------
1, 5, 4, -3, 10, 11,-1,-11,-8, 1, 4, -1, 0, 0

计算的结果长度 = 10（a的长度）+ 2 （M-1） =12

示例代码如下

import numpy as np


def get_blocks(signal, block_length):
    a_length = len(signal)
    padding_length = block_length - 1
    segments = -(-a_length // block_length)
    result = np.zeros((segments, block_length + padding_length))
    for i in range(segments):
        b = signal[i * block_length: (i + 1) * block_length]
        p = block_length + padding_length - len(b)
        b = np.pad(b, (0, p), 'constant', constant_values=(0, 0))
        result[i] = b
    return result


def generate_convolve_matrix(h, length):
    m = np.zeros((length, length))
    if len(h) < length:
        h = np.pad(h, (0, length - len(h)), 'constant', constant_values=(0, 0))
    for i in range(length):
        m[i] = h
        h = np.roll(h, 1)
    m = np.transpose(m)
    return m


def convolve_blocks(blocks, h_matrix):
    result = np.zeros((blocks.shape[0], blocks.shape[1]))
    for i in range(blocks.shape[0]):
        result[i] = np.dot(h_matrix, blocks[i])
    return result


def convolve_overlap_add(a, h):
    a_length = len(a)
    h_length = len(h)
    o_length = h_length-1
    blocks = get_blocks(a, h_length)
    length = blocks.shape[1]
    hm = generate_convolve_matrix(h, length)
    cb = convolve_blocks(blocks, hm)

    result = np.zeros(length * length)
    read_offset = 0
    overlap_add = np.zeros(length)
    for i in range(cb.shape[0]):
        if i == 0:
            result[0:length] = cb[i]
            read_offset = h_length
        else:
            result[read_offset:read_offset + length] = overlap_add + cb[i]
            read_offset = read_offset + h_length
        overlap_add = cb[i][h_length:]
        overlap_add = np.pad(overlap_add, (0, length - len(overlap_add)), 'constant', constant_values=(0, 0))
    return result[:a_length+o_length]


a = np.asarray([1, 3, -1, 2, 5, 3, -2, -4, -2, 1])
h = np.asarray([1, 2, -1])


c1 = convolve_overlap_add(a, h)
print('convolve_overlap_add:',c1)


c2 = np.convolve(a, h)
print('np.convolve:', c2)