前言

本人在对野外仪器数据采集的信号做一些列数字信号处理的过程中，经常遇到卷积与循环卷积。现实中，野外仪器记录到的数据都是离散的！所以我们在对数据进行FFT变换时，都默认使用的是DFT(离散傅里叶变换)。

本文目标并非介绍卷积与循环卷积的详细内容与公式，而是侧重如何编写自己的程序实现这两种操作，并简要介绍两者在使用和结果上的区别。

卷积与循环卷积区别

明确说明一点：虽然这两种操作都是对"两个信号"的处理，但两者不是同一操作！
区别1：卷积可以处理连续信号和离散信号，循环卷积只用于处理离散信号(现实中都可用)；
区别2：卷积对a和b这两个信号的长度没有要求，循环卷积要求a和b信号长度一致；
区别3：设操作后的信号为c，卷积操作后c的长度为a+b-1，循环卷积长度为a(a=b)。

卷积如何使用

卷积的操作就是"循环乘积与加和"。下面结合一个具体的例子进行说明：
信号a：a = [1 3 4 9 8 7]
信号b：b = [2 4 6 3]
卷积：a*b

离散卷积计算流程图.png

对应的Matlab语句是conv(a,b)。但是用edit conv命令查看conv的源代码时候，感觉写的很复杂，不如根据上面计算的流程图(计算原理)写自己的代码：

a = [1 3 4 9 8 7];
b = [2 4 6 3];          % a*b  a在前b在后
a_length = length(a);
b_length = length(b);
r = zeros(1,a_length+b_length-1); % 记录卷积结果，总长度是a+b-1

for k = 1:a_length    % a做循环总长 因为a在前所以是a在移动!
    c = a(k)*b;       % a所有的数都会和b数组所有元素乘一遍，但各自在r中起作用的位置不同
    
    d = r(k:k+b_length-1);  % r(k:k+wb-1)很妙！提取当前a(k)元素会影响的r的区间
    d = d+c;                % 把当前a(k)影响的范围"对应"加进去
    
    r(k:k+b_length-1) = d;  % 把r当前受影响的区域更新(加进去)
end

fprintf('卷积结果为:');
r

注意：进行卷积计算的两个信号谁在前谁在后结果都是一样的，即a*b = b*a；但是上面的程序还是要区分一下"谁在前"的问题(小细节，注意一下即可)。

循环卷积如何使用

循环卷积的操作一般用"傅里叶逆变换"来计算。下面结合一个具体的例子进行说明：
信号a：a = [2,1,2,1,6,2,8,9]
信号b：b = [3,4,2,4,3,5,1,8]
循环卷积：a*b
注意：两个信号的长度一样

用自己编的DFT与IDFT进行循环卷积的Matlab程序：

clc; clear;

% 手动输入的信号:
fprintf('原始信号为:\n');
a = [2,1,2,1,6,2,8,9]
b = [3,4,2,4,3,5,1,8]
N = length(a);

WN = exp(-i*2*pi/N);    % 常数
WN_nk_a = zeros(N)+WN;  % a的WN_kn
WN_nk_b = zeros(N)+WN;  % b的WN_kn
xk_a = a';     % a时域信号振幅(列矩阵)
xk_b = b';     % b时域信号振幅(列矩阵)
E_a = zeros(N);   % a的辅助的E(WN_kn的幂,单独拿出来算)
E_b = zeros(N);   % a的辅助的E(WN_kn的幂,单独拿出来算)

%%% a,b信号傅里叶正变换即结果 %%%
for row = 0:N-1
    for cow = 0:N-1
        E_a(row+1,cow+1) = row*cow;
        E_b(row+1,cow+1) = row*cow;
        WN_nk_a(row+1,cow+1) = WN_nk_a(row+1,cow+1)^E_a(row+1,cow+1);
        WN_nk_b(row+1,cow+1) = WN_nk_b(row+1,cow+1)^E_b(row+1,cow+1);
    end
end

Xk_a = WN_nk_a * xk_a;
Xk_b = WN_nk_b * xk_b;

Xk_t = Xk_a.*Xk_b;

%%% T信号(循环卷积)傅里叶逆变换即结果 %%%
WN_nk_t = zeros(N)+WN;  % t的WN_kn
E_t = zeros(N);   % t的辅助的E(WN_kn的幂,单独拿出来算)

for row = 0:N-1
    for cow = 0:N-1
        E_t(row+1,cow+1) = -row*cow;
        WN_nk_t(row+1,cow+1) = WN_nk_t(row+1,cow+1)^E_t(row+1,cow+1);
    end
end

fprintf('a,b信号的循环卷积结果为:\n');
xk_t = real((WN_nk_t * Xk_t)/N)'

优势1：信号可用是任意长度(不用像FFT那种必须是2^N的长度)；
优势2：完全理解其内涵。

当然，全部用Matlab自带函数写就是下面：

fprintf('原始信号为:\n');
a = [2,1,2,1,6,2,8,9]
b = [3,4,2,4,3,5,1,8]

a_tmp = fft(a);
b_tmp = fft(b);  % 用FFT计算a与b信号的频域信号

r_tmp = a_tmp.*b_tmp;  % 对应元素相乘(点乘)
r = ifft(r_tmp)        % 反变换回去，就是a,b信号循环卷积结果(时域)