matlab循环卷积函数,一维卷积与循环卷积的使用(离散数据+具体例子)

前言

本人在对野外仪器数据采集的信号做一些列数字信号处理的过程中，经常遇到卷积与循环卷积。现实中，野外仪器记录到的数据都是离散的！所以我们在对数据进行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信号循环卷积结果(时域)