离散卷积算法与傅里叶|MATLAB-based

前言

熟悉卷积在相关领域的“妙用”之后，这里再从数学与计算机中讲讲离散卷积的相关算法以及离散傅里叶变换📕

离散卷积

与数学中常用的运算符（如： $+$ 、 $-$ 、 $\times$ 、 $\div$ ）一样，卷积从数学角度出发，也是一种运算方式

卷积的数学定义如下：

$f(t)*g(t)=(f*g)(t)=\int_{-\infty}^{+\infty}f(\tau)·g(t-\tau)d\tau$

当然，上式是连续函数（或者说是连续时间信号）的卷积公式，下面给出与之对应的离散版本：

$x[n]*h[n]=(x*h)[n]=\sum_{k=-\infty}^{+\infty}x[k]·h[n-k]$

在《信号分析》相关知识中，我们定义上式中 $h[n]$ 为一单位脉冲响应，而任意输入序列 $x[n]$ 的零状态响应我们就可以用二者的卷积来表示其结果 $y[n]$ 。

下面我们将详细讨论离散序列卷积的算法与程序实现。

手工计算

根据已给出的数学公式，对于有限长的序列我们可以很轻易地、按部就班地得出对应的卷积结果。

（有限长序列可以看成除有数值外其他索引下的值都为0的无限长序列）

如：给定序列：

$x[n]=[1,-1,3,2,4,3]\\ h[n]=[1,2,3,4,5,6]$

默认有： $n=0,1,...,5$

其手工计算过程如下：

🔵Step1:计算 $y[0]=\sum_{k=0}^5x[k]h[-k]$

可见，只有k取0时，有 $h[-k]=h[0]=1$ 有值，其他情况下h[-k]=0

故结果是： $y[0]=1$

🔵Step2:计算 $y[1]=\sum_{k=0}^5x[k]h[1-k]$

仅有 $k=0,1时，即h[0]=1,h[1]=2$ 被用上时乘积不为0.

因此结果为： $y[1]=1$

🔵Skip…

由卷积的数学表达式，我们可以看出在计算时我们以k作为“自变量”。

而式子 $h[n-k]$ 就相当于：将 $h[k]$ 翻转得到 $h[-k]$ 然后再向右平移n个单位 $h[-(k-n)]$

综上，我们得到图形化的直观理解手工算法的要点图示，如下：

🔵Output

结果为： $y[n]=[1,1,4,9,18,30,35,53,44,39,18]$

结论一：若x[n]长度为N，h[n]长度为M，则y[n]=(x*h)[n]的长度为：(N+M-1)

编程计算

结合数学式，我们可以发现离散卷积的计算就是简单的 “乘积+求和”的循环过程，因此，设计一个计算函数也是较为简单的。

下面以matlab为例，给出具体代码.

注意：1.matlab中索引是从1开始，而不是0开始，因此代码需要做出相应的调整

2.矩阵索引只能为正整数，不能超出设定好的长度，默认值为0的部分需要自己过滤

P.S. 此处贴出的代码是初级代码，经过算法优化的代码在后文

%matlab函数文件 myConv.m

function y = myConv(x,h)
    % 用于计算x和h的卷积结果
    N = length(x);
    M = length(h);
    y = zeros(1,N+M-1); %运用结论一，定义一个全0序列
                        %也就是一个1×(N+M-1)的矩阵
    for n = 1:N+M-1
        ans = 0;
        for k = 1:N
        	if n-k+1>0 && n-k+1 <= M  	%过滤+调整
            	ans = ans + x(k)*h((n+1)-k);	%调整
            end
        end
        y(1,n) = ans;
   	end
end

%matlab脚本文件

x = [1,-1,3,2,4,3];
h = [1,2,3,4,5,6];
y1 = myConv(x,h)  %调用自定义的函数进行处理
y2 = conv(x,h)   %调用matlab自带的卷积函数计算

为因为数据选取的避免特殊性，采用isequal()函数多次进行对比：

1
2
3

a = randi(100,[1,100]);   %随机生成长度为100、数值在1~100中选取的序列
b = randi(100,[1,2000]);
isequal(myConv(a,b),conv(a,b))  %判断二者结果是否相等

结果展示

✅：结果正确

事实上，除了上述较为直接明了的算法外，还要其他类似的计算方法。将在后文贴出并作适当解释。

什么？你以为自己被骗了？明明matlab自身就有conv()函数，为什么自己还要编写？
下面我们来讨论讨论conv函数与我们的区别！

我们先利用时间函数计算两个函数运行时长

% 为了说明问题，将序列长度取得足够长
a = randi(100,[1,1e4]);
b = randi(100,[1,2e4]);
timeit(@() myConv(a,b))
% ans = 0.9240
timeit(@() conv(a,b))
% ans = 0.0068

自定义的函数在速度上几乎被conv函数碾压！

当然我们明显能看出，自制的计算法要用到两次for循环，时间复杂度可想而知

而做算法的同学当然对此则更加敏感，比如在思考有没有什么方法能够优化这个算法呢？

The Answer Is Yes 👾

问：conv函数为什么那么快？

答：因为conv函数并没有直接根据卷积的定义进行计算，而是用到了一个数学上的恒等式(~~果然最后还是比数学~~)

$f(t)*g(t)=F(s)\times G(s)$

也就是说，时域卷积等于频域相乘

毕竟做单独的乘法肯定比一个个乘然后再累加快得多嘛~😵

问：如何让一个信号从时域变成频域？

答：对该信号做傅里叶变换

离散傅里叶变换

离散傅里叶变换（Discrete Fourier Transform，DFT）是《信号分析》中提供的一种信号分析的方法。可以将时域信号转变为频域信号进行信号研究。这里仅给出表达式，不做过多推导。

$X[k]=\sum_{n=0}^{N-1}x[n]W_N^{kn}$

其中， $X[k]$ 即是时域离散序列 $x[n]$ 的频域表示， $W_N=e^{-j\frac{2\pi}{N}}=exp\{-j\frac{2\pi}{N}\}$

不难发现，我们想要计算频域信号，只需一次循环即可。

即：

%此处仅给出伪代码
for n = 0:N-1
	nk = n*k
	ans = x(n)*wn^nk
end

关于Matrix

matlab作为一款出色的数学编程工具，最有意义的一件事那就是引入了矩阵的概念。

因此，其实前面我们用到的类似C语言的“数组”/序列其实都是一个 $（1\times n）$ 的矩阵

用之前的代码来实现相关功能，大材小用了（雾）

请回忆矩阵乘法的基本概念，我们就可以发现矩阵乘法也是一个乘积求和的过程，因此我们可以借助matlab自带的矩阵乘法机制，优化我们的代码，使得代码更加简洁、快捷和有效

这也是我仅给出一般的求DFT伪代码的原因之一。（~~绝对不是因为懒~~）

下面，我将借助Martix给出matlab的DFT实现代码。

DFT的MATLAB实现

复习一下矩阵乘法：

$A_{m\times n}\times B_{n\times m}= \begin{bmatrix} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \cdots & a_{mn} \end{bmatrix} \times \begin{bmatrix} b_{11} & \cdots & b_{1m} \\ \vdots & \ddots & \vdots \\ b_{n1} & \cdots & b_{nm} \end{bmatrix} \\=\begin{bmatrix} \sum_{j=1}^na_{1j}b_{j1}&\cdots &\sum_{j=1}^na_{1j}b_{jm}\\ \vdots & \ddots &\vdots\\ \sum_{j=1}^na_{mj}b_{j1}&\cdots &\sum_{j=1}^na_{mj}b_{jm} \end{bmatrix}\\ =C_{m\times m}$

也就是说，对于C的某一行i某一列j有：

$C_{(i,j)}=\sum_{t=1}^na_{it}b_{tj}$

这样一来，我们就可以将 $X[k]$ 的表达式与之对应。

$X$ 就是此处的 $C$ ，且为一维行向量， $X[k]$ 就是第一行第k列的值。

同理， $x[n]$ 也是一维行向量。为了计算，我们需要将 $W_N^{kn}$ 作为方阵处理，由此得到：

$\textbf{X}=\textbf{x}\times\textbf{W}_N$

其中，有：

$\textbf{x}=\begin{bmatrix} x_0,x_1,...,x_n \end{bmatrix} \\\textbf{W}_N=\begin{bmatrix} W_N^{0\times 0}&\cdots&W_N^{0\times N-1} \\\vdots&\ddots&\vdots \\W_N^{N-1\times 0}&\cdots&W_N^{N-1\times N-1} \end{bmatrix}$

（注：此处都是以起始为0计算的）

可见，W矩阵的指数部分其实就是由0~(N-1)的数组合而成，可以使用如下方式表示之：

$\textbf{n}=[0,1,...,N-1] \\\textbf{n}'=\begin{bmatrix} 0 \\1 \\2 \\\vdots \\N-1 \end{bmatrix} \\\textbf{ans}=\textbf{n}\times\textbf{n}'$

于是，我们就可以得到如下的计算方法了。

function Xk = mydft(xn, N)
    % DFT 给定离散序列xn和参数N，求其傅里叶变换
    n = [0:1:N-1];
    k = n;
    WN= exp(-j*2*pi/N);
    nk = n' * k;            
    Xk = xn(1:N) * WN.^nk;
end

注：n'表示n的转置；.^表示每一项都进行指数运算

反变换(iDFT)

结合逆变换的表达式，我们可以很轻易得到逆变换的函数。

$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}$

function xn = myidft(Xk, N)
    % iDFT 给定离散序列Xk和参数N，求其傅里叶反变换
    n = [0:1:N-1];
    k = n;
    WN= exp(-j*2*pi/N);
    nk = n' * k;            
    xn = (Xk(1:N) * WN.^(-nk))/N;
end

事实上，前面的卷积表达式也可以利用矩阵的方式简化代码。此处不再赘述，其他类型的代码会放在后文

再谈离散卷积

起点不为0情况

其他的卷积算法

利用FFT的卷积计算(首推)

问：什么是FFT？

答：快速傅里叶变换(Fast Fourier Transform)的简称。是计算机领域内基于DFT设计并优化的更好的算法

关于FFT将在新文章中继续探讨~

之前我们已经介绍了可以利用公式：

$f(t)*g(t)=F(s)\times G(s)$

计算卷积。

因此我们可以直接调用matlab内自带的函数fft()来进行计算。

对，FFT也自带得有（笑）

function y = myNewConv(x,h)
    %借助fft实现快速离散卷积
    
    %先将二者填充0使得矩阵维度相同
    xn = [x,zeros([1,length(h)-1])];
    hn = [h,zeros([1,length(x)-1])];
    y = ifft(fft(xn).*fft(hn));
end