离散傅里叶变换(DFT)

Jaihk662

已于 2023-10-10 21:40:35 修改

阅读量4.1k

点赞数 3

分类专栏： # 数论文章标签：数论

于 2022-11-21 16:43:26 首次发布

本文链接：https://blog.csdn.net/Jaihk662/article/details/127961899

版权

数论专栏收录该内容

71 篇文章 1 订阅

订阅专栏

接上文：傅里叶级数与傅里叶变换

二、“虚拟世界”中的傅里叶变换

真实世界是连续的，可是计算机永远只能描述离散的点，采集离散的信号

现在的你，已经完全理解了傅里叶变换：

$F(\omega)=\int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t$ （其中 $w$ 为可变频率 $2\pi f$ ）

以及傅里叶变换的逆变换 $f(t)=\frac{1}{2 \pi} \int_{-\infty}^{+\infty} F(\omega) e^{i \omega t} d \omega$

好了，你有一台电脑，你想要让电脑也理解傅里叶变换：当然了：电脑只能描述和采集离散的信息，就别说指数空间了，电脑对积分也是完全不懂的，也永远不可能懂

当然这也并不意味没有办法了，前面提到过，傅里叶的目的实际上是计算信号幅值在频域中的分布密度，也可以说是为了求得频谱密度函数，同时在描述傅里叶变换时：积分上下无限同时意味着要对所有时间上的值进行考虑，即时间间隔为 $\infty$ 时，其极限是多少

那计算机不明白也没关系，因为只要想办法让计算机去得到/接近我们想要结果，过程其实可能没有那么重要

2.1 采样与冲激函数(Dirac δ 函数)

计算机解决积分的一个非常暴力的方式就是：将积分范围内所有可以取到的值，一个一个丢进去算出结果，最后再加在一起求个平均，只不过让我们取无数个数进去算结果必然是不可能的，因此我们可以每隔一段距离去算一个结果，以做到得到近似的答案

而这个计算过程就有采样的影子，用更通用的话来讲，采样就是在离散世界里描述连续的图像信息的手段，想想纹理采样的过程？是不是就是这个意思

采样的专业解释

想要用数学方法描述采样，需要先引入冲激函数：

当 $x \neq 0$ 时， $\delta(x)=0$ ，且满足 $\int_{-\infty}^{\infty} \delta(x) d x=1$ ，由于不是复杂的理学论文，你其实可以简单的把这个函数理解成只有当 x = 0 时存在值 1

根据他的定义，可以得到 $\int_{-\infty}^{+\infty} \delta\left(t-t_0\right) f(t) d t=f\left(t_0\right)$ ，即能够筛选出连续函数在 $t_0$ 处的取值，起到采样的作用

不过事实上我们要的采样是每隔一段距离周期性的去取一个值，而非只取 $t_0$ 处的函数结果，因此最后更应该是 $f_s=\sum_{n=-\infty}^{\infty} f(t) \delta\left(t-n T_s\right)$

其中 $T_s$ 为采样周期： $\delta_s(t)=\sum_{n=-\infty}^{\infty} \delta\left(t-n T_s\right)$ 的图像更像是这样：是由一个个脉冲序列组成

好了，讲完了，这节用一句话说就是对于连续函数 $f(t)$ ，计算机能计算的其实是 $f_s=\sum_{n=-\infty}^{\infty} f(t) \delta\left(t-n T_s\right)$ ，其中 $\delta(x)$ 为冲激函数(Dirac δ 函数)

2.2 时域与频域的离散化计算

Games101 里面对于采样也有个生动形象的例子，而前面讲述的内容正好对应着下图中左侧时域的离散采样，左下图 e 对应的 S(t) 对应上面公式中的 $f(t) \delta\left(t-n T_s\right)$

然后就是把 $f(t) \delta\left(t-n T_s\right)$ ，代入到傅里叶变换中，有

$X(\omega)=\int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t =\sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} f(t) \delta\left(t-n T_s\right) e^{-i w t} d t$

很明显，只有当 $t=n T_s$ 时，右边的 $f(t) \delta\left(t-n T_s\right) e^{-i w t}$ 才对积分结果有贡献，因此

$X_s(\omega)=\sum_{n=-\infty}^{\infty} f\left(n T_s\right) e^{-i w n T_s}$

不过到这采样的规则还没有定，也就是说还要敲定①采样间隔应该是多大，可以不可以量化这个间隔；②不可能进行无穷次采样，我们要有个采样范围，否则你计算得到的频域图像依然是连续的，而且你时域采样个无数次也不现实

这也很简单，一个采样思路是这样的：

不管你是不是周期函数，我就限死一个范围，只在这个范围按照一定间隔采样
对于范围外的函数图像，我们可以假设其为采样范围内图像的延拓（覆盖掉原先的图像），在这种假设下，原函数就会成为一个周期函数，可以为我们后面的求解提供很大的便利，并且同时可以解决频率连续的问题

例如上图，如果我们只对 [-3, 3] 范围内的函数进行采样，可以直接将原函数转成像右图那样的周期函数

这样假设我们采样间隔为 $T_s$ ，采样 N 次，采样周期 $T_0=N T_s$ ，那么在一段频率内，离散信号的表达式就为：

$f_s=\sum_{n=-\infty}^{\infty} f(t) \delta\left(t-n T_s\right)$ → $x_s(t)=\sum_{n=0}^{N-1} f(t) \delta\left(t-n T_s\right)$

其傅里叶级数

$X\left(\omega\right)=\frac{1}{T_0} \int_0^T\left(\sum_{n=0}^{N-1} f(t) \delta\left(t-n T_s\right)\right) e^{-i w t} d t$

和上面计算方式一样，可以交换积分次序： $X\left(\omega\right)=\frac{1}{T_0} \sum_{n=0}^{N-1} \int_0^T f(t) \delta\left(t-n T_s\right) e^{-i \omega t} d t$

然后只有当 $t=n T_s$ 时，右边的积分才对结果有贡献： $X\left(\omega\right)=\frac{1}{T_0} \sum_{n=0}^{N-1} f(n T_s) e^{-i\omega n T_s}$

再根据 $T_0=N T_s$ ， $\omega=\frac{2 \pi}{T_0}$ ，代入进去得到： $X\left(\omega\right)=\frac{1}{NT_s} \sum_{n=0}^{N-1} f(n T_s) e^{-i\frac{2\pi}{N} n}$

当然因为是离散采样，用 $X\left(\omega\right)$ 描述最终的结果不是很准确，回到前面的那张图，我们是要采样 N 次，每次采样间隔 $T_s$ ，那么第想要得到第 $k(0{\leq}k{<}N)$ 次采样到第 $k+1(0{\leq}k{<}N)$ 次采样中间这段连续的结果，就同等于我们需求的是 $X[k] = X\left(k\omega\right) T_s$ ，同理，我们可以规定 $x[n]=f\left(n T_s\right)$ ，这样就能得到最终的离散傅里叶形式：

$X[k]= \sum_{n=0}^{N-1} x[n] e^{-i\frac{2\pi}{N}kn}$ ， $k(0{\leq}k{{<}}N)$

这即是计算机求解傅里叶变换的标准方式

三、快速傅里叶变换(FFT)

3.1 预备性质与数学原理

很明显，计算其中一项的 $X[k]$ ，复杂度为 O(N)，但是我们往往需要计算所有的 $X[k]$ ， $k(0{\leq}k{<}N)$ ，复杂度就为 O(N²)，而快速傅里叶变换(FFT)能在 O(NlogN) 的时间复杂度内得出这堆结果

3.1.1 单位复 N 次方根

把 DFT 指数部分抽出来： $e^{-i\frac{2\pi}{N}kn}$ ，为了方便后续令 $W_N=e^{-i\frac{2{\pi}}{N}}$ ， $W_N^n$ 在复数上的表示如下：

对于图像中的例子 N = 16，相邻两点的夹角 $\text { angle }=\frac{2 \pi}{N}$ ，不用任何公式推导，就能发现只要 N 为偶数，那么图中的点就可以通过与轴心连线两两配对：即得到 $W_N^{n+\frac{N}{2}}=-W_N^{n}$

除此之外，我们对图中所有点的构成点集 $\left[1, W_N^1, W_N^2, \ldots, W_N^{n-1}\right]$ ，全部做平方，得到的点集

$\left[1, W_N^2, W_N^4, \ldots, W_{N}^{2(n-1)}\right]$ 图像就如下：

也一样不用公式推导，就有两个性质跃然纸上：

点集中后半部分的点会和前半部分的点完全重叠，也就是说新的点集中少了刚好一半的元素： $\left[1, W^2, W^4, \ldots, W^{2(\frac{n}{2}-1)}\right]$ ，并且其等价于 $\left[1, W_\frac{N}{2}^1, W_\frac{N}{2}^2, \ldots, W_\frac{N}{2}^{n-1}\right]$
新的点集大小如果还是2的倍数（也就是原先的点集大小 N 是4的倍数），那么依然满足前面提到的对称性 $W_N^{nk+\frac{N}{2}}=-W_N^{nk}$

如果原先的点集数量 N 满足 $N=2^p$ ，即2的指数幂，那么图像就可以一直递归下去，每次的分量都会少一半，在此算法时间复杂度中 O(nlogn) 的 "log" 也浮出水面

3.1.2 FFT 蝶形算法

还记得我们要求的东西嘛： $X[k]= \sum_{n=0}^{N-1} x[n] e^{-i\frac{2\pi}{N}kn}$ ， $k(0{\leq}k{{<}}N)$

为了能够用上前面的性质，我们可以假设 N 一定是二次幂，如果真实的 N 不是二次幂，那就强行补齐，并将后面要补的值 x[n] 全部当成0（关于这个方案是否是可行的当然也需要严谨证明，不过考虑到篇幅，并且文章面对的读者也不是理学研究生，这里就也做省略，有兴趣的可以参考这篇文档中关于 zero padding 的介绍部分）

然后把 DFT 根据 n 的奇偶分成两部分相加：其中 2r = n 表示偶数，2r + 1 = n + 1 表示奇数：

$X[k]=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_N^{2 r k}+\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r+1) W_N^{(2 r+1) k} \\ =\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_N^{2 r k}+W_N^{ k}\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r+1) W_N^{2 r k}$

再根据前面 3.1.1 的结论①： $\left[1, W^2, W^4, \ldots, W^{2(\frac{n}{2}-1)}\right]$ 等价于 $\left[1, W_\frac{N}{2}^1, W_\frac{N}{2}^2, \ldots, W_\frac{N}{2}^{n-1}\right]$ ，有

$W_n^2 = W_\frac{N}{2}^1$ ，上式可变为：

$X(k)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_{\frac{N}{ 2}}^{r k}+W_N^k \sum_{r=0}^{\frac{N}{ 2}-1} x(2 r+1) W_{\frac{N}{ 2}}^{r k}$

再令 $G(k)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_{\frac{N}{ 2}}^{r k}$ ， $H(k)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r+1) W_{\frac{N}{ 2}}^{r k}$

到此为止，一个采样 N 次的 DFT，就可以拆成两个采样 N/2 次的 DFT，分别为偶数点采样 G(k)，和奇数点采样 H(k)

但是还没有结束，前面 3.1.1 的结论②： $W^{nk+\frac{N}{2}}=-W^{nk}$ 还没有用上，显然

$G(k+\frac{N}{2})=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_{\frac{N}{ 2}}^{r (k+\frac{N}{2})}=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_{\frac{N}{ 2}}^{r k}W_{\frac{N}{ 2}}^{\frac{N}{ 2}}=G(k)$

同理： $H(k+\frac{N}{2})=H(k)$ 也成立

这样最终结论就出来了：

对于 $X(k)=G(k)+W_N^k H(k)$ ，有

$X(k)=G(k)+W_N^k H(k)$ ， $0{\leq}k{<}\frac{N}{2}$

$X(k)=G(k-\frac{N}{2})-W_N^{k-\frac{N}{2}} H(k-\frac{N}{2})$ ， $\frac{N}{2}{\leq}k{<}N$

也就是说，我们只需要求出前一半的 G(k) 和 H(k)，就能得出全部的 X(k)，而对于 G(k) 和 H(k)，它们本质上也是 DFT，但是规模都只有 X(k) 的一半，因此可以继续递归，直到 N = 1 结束

你可能会感受到一阵眼花，这是参透新事物的重要一环

没事，还有例子，如果 N = 4，那么就有：

$X(0)=G(0)+W_4^0 H(0) \quad X(1)=G(1)+W_4^1 H(1) \\ X(2)=G(0)-W_4^0 H(0) \quad X(3)=G(1)-W_4^1 H(1)$

其中 G 和 H 就是两个点的 DFT，并且此时我们要求的就是 G(0)，G(1)，H(0)，H(1)

此时 G 和 H 对应的 N = 2，对于 G(0) 和 H(0) 可以直接带入前面的公式：

$G(k)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r) W_{\frac{N}{ 2}}^{r k}$ ， $H(k)=\sum_{r=0}^{\frac{N}{ 2}-1} x(2 r+1) W_{\frac{N}{ 2}}^{r k}$ ，有

$G(0) = x(0)+W_1^0x(2) \quad H(0) = x(1)+W_1^0x(3)$

一样根据对称性：

$G(1) = x(0)-W_1^0x(2) \quad H(0) = x(1)-W_1^0x(3)$

搞定！

3.2 速解多项式乘法：从另一个角度看 FFT

对多项式 $f(x)=\sum_{i=0}^n a_i x^i$ ， $g(x)=\sum_{j=0}^m b_j x^j$ ，定义其乘积 fg 为

$(f g)(x)=\left(\sum_{i=0}^n a_i x^i\right)\left(\sum_{j=0}^m b_j x^j\right)$ ，现在给你 n, m，以及 f 和 g 的每一项系数，要求你求解它们乘积的多项式的每一项系数，当然 n 和 m 都是 100000 级别，因此时间复杂度要控制在 nlogn 级别，请问该如何编写相应的程序？

对于视 FFT 为基本操作的 ACMer 而言，这种题目真的再简单不过了，当然对于刚了解 FFT 的初学者，可能还是很难发觉这题和 FFT 或是 DFT 有什么关系

3.2.1 何为卷积

在此之前，先介绍一种函数运算关系：卷积

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

不过由于我们不做信号分析，研究的内容更偏向于图像处理，因此我们主要关心的应该是卷积方程的离散形式： $(f * g)(n)=\sum_{\tau=-\infty}^{\infty} f(\tau) g(n-\tau)$

网上一个经典的例子就是抛骰子：连续投掷两枚骰子，它们点数之和为4的概率是多少：

设 f(x) 为投掷第一枚骰子点数为 x 的概率，g(x) 为投掷第二枚骰子点数为 x 的概率，那么最后概率就是： $(f * g)(4)=\sum_{m=1}^3 f(4-m) g(m)=f(1) g(3)+f(2) g(2)+f(3) g(1)$

就这么简单：你甚至不必去分析这个公式，只需要把“卷积”理解为是一种特殊的计算方式就好了：即用一组函数进行一种特定的线性组合算法（加权叠加）得到一个新的目标函数

除此之外，卷积还有一个性质：时域的卷积可以转换为频域的乘积，反之亦然（因篇幅问题，关于这个定理的证明从略）：

$F[f(t) * g(t)] =F[f(\tau)] F[g(t-\tau)] =F_1(w) F_2(w)$

$2 \pi f(t) g(t) \longleftrightarrow F_1( \omega) * F_2( \omega)$

好了，然后再看回前面的多项式问题：

首先是对于任意 n+1 次多项式 $f(x)=\sum_{i=0}^n a_i x^i$ ，只需知道每一项的系数就可以唯一确定 $f(x)$ ，故 $f(x)$ 可以被写作 $\left(a_0, a_1, a_2, \ldots, a_n\right)$ 的形式，即系数的有序排列，这就是多项式的系数表示法

其次对于多项式 $f(x)=\sum_{i=0}^n a_i x^i$ ， $g(x)=\sum_{j=0}^m b_j x^j$ 的最终乘积结果 $(f g)(x)$ ，它的第 n 项 $c_nx^n$ 的系数就等于 $a_0b_n+a_1b_{n-1}+ \ldots+a_nb_0$ ，也就是说，其实本质上计算多项式乘法，就是在计算多项式系数表示法对应的离散函数的卷积

3.2.2 nlogn 的卷积计算

回到题目上，先忘掉 FFT 相关的内容，考虑多项式本身的性质：对于 n 次多项式上 n+1 个不同的点能唯一确定这个多项式，也就是说 $f(x)$ 也可以被写作

$\tau(f) =\left\{\left(x_0, f\left(x_0\right)\right),\left(x_1, f\left(x_1\right)\right), \ldots,\left(x_n, f\left(x_n\right)\right)\right\}$ 点值表示

如果再令 $\tau(g)=\left\{\left(x_0, g\left(x_0\right)\right),\left(x_1, g\left(x_1\right)\right), \ldots,\left(x_m, g\left(x_m\right)\right)\right\}$ ，则有

$\tau(f g)=\left\{\left(x_0, f\left(x_0\right) \cdot g\left(x_0\right)\right),\left(x_1, f\left(x_1\right) \cdot g\left(x_1\right)\right), \cdots,\left(x_{n+m}, f\left(x_{n+m}\right)\right)\right\}$

这样的话，我们就可以大致得出解题的步骤了（为了简化问题，我们假设 $f(x)$ 和 $g(x)$ 的项数是相同的，都为 N，且 N 为一个二次幂数）：

求出 $f(x)$ 和 $g(x)$ 的点值表示法（暴力算法复杂度 O(N²)）
对于得到的每个点值，x 坐标不变，y 坐标相乘，这样就可以得到一个新的点值序列，这个新的序列正式我们要求的多项式乘积结果 $(f * g)(n)$ 的点值表示 $\tau(f g)$ （复杂度O(N)）
再根据性质，我们可以根据 $\tau(f g)$ 逆推出 $(f * g)(n)$ ，得出答案（未知算法??）

然后就是设计步骤①和步骤③的算法

在这里再推荐一个视频：【官方双语】奥数级别的数数问题，虽然内容和 FFT 无关，但是也能给你带来不少启发，感受到数学之美

不过在此之前，还是留意一下步骤②，不知道你发现了什么？步骤②是个乘积算法，而前面我们提到过：计算多项式乘法，就是在计算多项式系数表示法对应的离散函数的卷积，而这正印证了前面的定理：时域的卷积可以转换为频域的乘积……

没错对了，不可能这么巧合，我们是不是就可以说：

多项式系数表示法 $f(x)=\left(a_0, a_1, a_2, \ldots, a_n\right)$ 对应着频域
而多项式点值表示法 $tau(f)=\left\{\left(x_0, f\left(x_0\right)\right),\left(x_1, f\left(x_1\right)\right), \ldots,\left(x_n, f\left(x_n\right)\right)\right\}$ 正是前者的时域表示

这样目标就很明确了，求出点值表示法的每一项 $\left(x_i, f\left(x_i\right)\right)$

还等什么，直接对 $f(x)=\left(a_0, a_1, a_2, \ldots, a_n\right)$ 上 DFT $a_k= \sum_{n=0}^{N-1} f[n] e^{-i\frac{2\pi}{N}kn}$ 啊

令 $w_N = e^{-i\frac{2\pi}{N}}$ ，DFT 转换成矩阵形式就是：

$\left[\begin{array}{ccccc} 1 & 1 & 1 & \cdots & 1 \\ 1 & \left(\omega_N^1\right)^1 & \left(\omega_N^1\right)^2 & \cdots & \left(\omega_N^1\right)^{n-1} \\ 1 & \left(\omega_N^2\right)^1 & \left(\omega_N^2\right)^2 & \cdots & \left(\omega_N^2\right)^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \left(\omega_N^{n-1}\right)^1 & \left(\omega_N^{n-1}\right)^2 & \cdots & \left(\omega_N^{n-1}\right)^{n-1} \end{array}\right]\left[\begin{array}{c} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{array}\right]=\left[\begin{array}{c} f\left(\omega_N^0\right) \\ f\left(\omega_N^1\right) \\ f\left(\omega_N^2\right) \\ \vdots \\ f\left(\omega_N^{n-1}\right) \end{array}\right]$

这就相当于我们把 $w_N^k = e^{-i\frac{2\pi}{N}k}$ 代进原 $f(x)$ 以求得其点值表示法

好了，别忘了这个过程是可以 FFT 加速的，也就是说我们能够能在 nlogn 的时间复杂度内求得上述步骤①对应的两个函数的点值表示

因此，最终的解体步骤应如下：

对于 $f(x)$ 和 $g(x)$ 的系数表示法，经 DFT 可得到其点值表示法（FFT 复杂度 O(NlogN)）
对于得到的每个点值，x 坐标不变，y 坐标相乘，这样就可以得到一个新的点值序列，这个新的序列正式我们要求的多项式乘积结果 $(f * g)(n)$ 的点值表示 $\tau(f g)$ （复杂度O(N)）
将最终的结果 $\tau(f g)$ 进行 IDFT 逆离散傅里叶变换，即可求得最终的 $(f * g)(n)$ ，得出答案（FFT 复杂度 O(NlogN)）

对于步骤③，IDFT 的矩阵可以由 DFT 的矩阵求逆得出，也可以直接套公式：

$\mathbf{C}=\frac{1}{n}\left[\begin{array}{ccccc} 1 & 1 & 1 & \cdots & 1 \\ 1 & \left(\omega_N^{n-1}\right)^1 & \left(\omega_N^{n-1}\right)^2 & \cdots & \left(\omega_N^{n-1}\right)^{n-1} \\ 1 & \left(\omega_N^{n-2}\right)^1 & \left(\omega_N^{n-2}\right)^2 & \cdots & \left(\omega_N^{n-2}\right)^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \left(\omega_N^1\right)^1 & \left(\omega_N^1\right)^2 & \cdots & \left(\omega_N^1\right)^{n-1} \end{array}\right]$

搞定！其实这个思考解题流程，是有要求对傅里叶变换及其意义有一定了解的，你也可以看出来，我们在解题的过程中有些先入为主。如果你之前没有了解过傅里叶变换，但是却又急切的要解出这道题，那么可以直接从另一个角度思考，在此继续推荐视频：快速傅里叶变换(FFT)——有史以来最巧妙的算法？

3.3 FFT 的程序实现

https://www.luogu.com.cn/problem/P3803

#include<stdio.h>
#include<algorithm>
#include<complex>
using namespace std;
complex<double> a[1 << 22], b[1 << 22];

//长度为 n 的 DFT
void FFT(complex<double>* a, int n, int inv)
{
    if (n == 1)
        return;
    int mid = n / 2;
    complex<double> *A1, *A2;
    A1 = new complex<double>[mid];
    A2 = new complex<double>[mid];
    //采样 N 次的 DFT，可以拆成两个采样 N/2 次的 DFT，分别为偶数点采样 A1(k)，和奇数点采样 A2(k)
    for (int i = 0; i < n; i += 2)
    {
        A1[i / 2] = a[i];
        A2[i / 2] = a[i + 1];
    }
    FFT(A1, mid, inv);
    FFT(A2, mid, inv);
    complex<double> w(1, 0), w1(cos(2 * acos(-1.0) / n), inv * sin(2 * acos(-1.0) / n));
    for (int i = 0; i < mid; i++)
    {
        a[i] = A1[i] + w * A2[i];
        a[i + n / 2] = A1[i] - w * A2[i];
        w *= w1;
    }
    delete []A1;
    delete []A2;
}

int main()
{
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i++)
    {
        double x;
        scanf("%lf", &x);
        a[i].real(x);
    }
    for (int i = 0; i <= m; i++)
    {
        double x;
        scanf("%lf", &x);
        b[i].real(x);
    }
    int len = 1 << max((int)ceil(log2(n + m)), 1);
    FFT(a, len, 1);         //DFT
    FFT(b, len, 1);
    for (int i = 0; i <= len; ++i)
        a[i] = a[i] * b[i];
    FFT(a, len, -1);        //DTFT
    for (int i = 0; i <= n + m; ++i)
        printf("%.0f ", a[i].real() / len + 0.000001f); 
    return 0;
}