Skip to content

Latest commit

 

History

History
503 lines (412 loc) · 20.4 KB

File metadata and controls

503 lines (412 loc) · 20.4 KB

3.1.1 一维傅立叶(1D-FT)与一维离散傅立叶变换(1D-DFT)

信号学中,将沿平面分布的信号称为一维信号(1D Signal),例如音频信号。

一维傅里叶变换,能够将一组满足狄利克雷条件(Dirichlet Theorem)的一维信号分解到周期性复指数波(Complex Exponential Wave)构成的二维向量空间。

从傅里叶级数(FS)到傅里叶变换(FT)

狄利克雷条件 最初被用作傅里叶级数(FS [Fourier Series])在三角函数域上进行分解的充分不必要条件 [2] [7]。在狄利克雷条件描述中,如果选定分析的周期信号 同时满足:

  • 【单周期内,连续或存在有限个第一类间断点】;
  • 【单周期内,存在有限数目的极大值与极小值】;
  • 【单周期内,绝对可积】;

则,此周期信号就一定存在傅里叶三角级数的分解表示。

如果记周期信号函数 $$s(t)$$ 的波长(周期)为 $$T$$ ,角频率(角速度)为 $$\tfrac{2\pi}{T}$$ 。则以信号函数波长 $$T$$ 做可变 $$n \in [0, \ N]$$ 等分(即步长 $$Step = \tfrac{1}{N}$$ )选取分离函数。有分离函数(周期)为 $$\tfrac{T}{n}$$ ,角频率(角速度)为 $${\omega_n} = \tfrac{2\pi n}{T}$$ 。原周期信号函数 $$s(t)$$ 就可以被分解为:

$$ {\displaystyle \begin{aligned} s(t) &= \frac{1}{N} \sum_{n =0}^{N} a_n \cdot cos(\tfrac{2\pi n}{T}t)\ \ \ \ \ \ +\ \ \ \ \ \frac{1}{N} \sum_{n =0}^{N} b_n \cdot sin(\tfrac{2\pi n}{T}t) \\ a_n &= \int_{-\tfrac{T}{2}}^{+\tfrac{T}{2}} s(t) \cdot cos(\omega_n t) \ dt \ \ \ \ \\ b_n = \int_{-\tfrac{T}{2}}^{+\tfrac{T}{2}} s(t) \cdot sin(\omega_n t) \ dt \\ \end{aligned} } $$

如果我们对函数周期进行平移,将区间从 $$(-\tfrac{T}{2},\ +\tfrac{T}{2})$$ 偏移 $$+\tfrac{T}{2}$$ ,即变换到 $$(0,\ T)$$ ,使原周期信号函数 $$s(t)$$ 偏移为奇函数(即 $$s(-t) = - s(t)$$ ),而奇函数式可证明是不需要余弦函数项的。此时,就可以进一步化简 $$s(t)$$ 为存粹正弦函数表示:

$$ {\displaystyle \begin{aligned} s(t) &= \frac{1}{N} \sum_{n =0}^{N} b_n \cdot sin(\tfrac{2\pi n}{\lambda}t) = \frac{1}{N} \sum_{n =0}^{N} b_n \cdot sin(\omega_n t) \\ b_n &= \int_{0}^{T} s(t) \cdot sin(\omega_n t) \ dt \\ \end{aligned} } $$

简化表示 $${\omega_n}$$$${\omega}$$ ,当我们将傅里叶级数从三角函数域,扩展到复变函数域时,基底函数由正余弦函数变为了以 $${\displaystyle \begin{aligned} \lambda = \tfrac{2 \pi}{\omega} = \tfrac{T}{n}\ \end{aligned} }$$ 为周期(波长)的复指数函数 $${\displaystyle \begin{aligned} {\mathcal {S}}_{\omega}(t) = e^{i\omega t}\ \end{aligned} }$$ 。信号函数 $$s(t)$$ 的分解函数就可以表示为:

$$ {\displaystyle \begin{aligned} s(t) &= \frac{1}{N} \sum_{n = 0}^{N} \hat{s}(\tfrac{2\pi n}{T}) \cdot e^{i \tfrac{2\pi n}{T}t}
= \frac{1}{N} \sum_{\omega = 0}^{\omega_N} \hat{s}(\omega) \cdot e^{i \omega t} \ &= \frac{1}{N}\sum_{n = 0}^{N} \hat{s}(\omega) \cdot {\mathcal {S}}{\omega}(t) \ \hat{s}(\omega) &= \int{0}^{T} s(t) \cdot e^{-i \omega t} \ dt \ \end{aligned} } $$

根据 欧拉公式(Euler's Formula) 可知 $${\displaystyle \begin{aligned} e^{ix} = cos(x) + i \cdot sin(x) \end{aligned} }$$ , 带入上式有:

$$ {\displaystyle \begin{aligned} s(t) &= \frac{1}{N}\sum_{n = 0}^{N} \hat{a}{\omega} \cdot cos(\omega t) + i \cdot \hat{b}{\omega} \cdot sin(\omega t)\ \hat{a}{\omega} &= \hat{s}(-\omega) + \hat{s}(\omega) \quad \quad \hat{b}{\omega} = \tfrac{1}{i} \cdot (\hat{s}(-\omega)-\hat{s}(\omega)) \end{aligned} } $$

转换到欧氏空间下的三角函数表示 $${\mathcal {S}}{\omega}(t)$$ ,记构成原信号函数 $$s(t)$$ 的复指数函数 $${\mathcal {S}}{\omega}(t)$$ 的初相为 $$\angle\phi_{\omega}$$ ,振幅为 $$A_{\omega}$$ ,则:

$$ {\displaystyle \begin{aligned} {\mathcal {S}}{\omega}(t) : \quad \angle\phi{\omega} = \arctan(\tfrac{\hat{a}{\omega}}{\hat{b}{\omega}}) \quad A_{\omega} = \sqrt{ (\hat{a}{\omega}) ^2 + (\hat{b}{\omega}) ^2 } \ \end{aligned} } $$

同三角函数域的情况,复变函数域下的傅里叶级数仍然可以进一步精简。我们仍然需要对原函数 $$s(t)$$ 平移 $$+\tfrac{\lambda}{2}$$ 并将周期变换到 $$(0,\ \lambda)$$ ,使 $$s(t)$$ 表现为奇函数。由于原信号函数 $$s(t)$$ 必为实函数的特性,会使得 $$a_{\omega}$$$$b_{\omega}$$ 互为共轭复数。因此在奇函数条件下, $$a_{\omega}$$$$b_{\omega}$$ 表现为符号相反的纯虚数,此时:

$$ {\displaystyle \begin{aligned} \hat{a}{\omega} &= 1 \cdot [\hat{s}(-\omega) + \hat{s}(\omega)] = 0 \ \ \ \ \ \hat{b}{\omega} = \tfrac{1}{i} \cdot [\hat{s}(-\omega)-\hat{s}(\omega)] = \tfrac{2}{i} \cdot \hat{s}(-\omega) \ s(t) &= \frac{1}{N} \sum_{\omega =0}^{\omega_N} \ \ \ \ 0 \cdot cos(\omega t) \ \ \ \ \ + \ \ \ \ i \cdot (\tfrac{2}{i} \cdot \hat{s}(-\omega)) \cdot sin(\omega t) \ &= \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{1}{N}\sum_{n = 0}^{N} \hat{s}(-\omega) \cdot sin(\omega t) \ \end{aligned} } $$

如果我们将 $$\hat{s}(-\omega)$$ 的负号划入公式,并将离散级数扩展到原信号函数 $$s(t)$$ 的连续实数空间上以积分形式表示。则 $$s(t)$$$$\hat{s}(-\omega)$$ 的关系就展现为:

$$ {\displaystyle \begin{aligned} s(t) &= \frac{1}{N}\int_{0}^{N} \hat{s}(\omega) \cdot sin(\omega t) \ d{n} \\ \hat{s}(\omega) &= \int_{0}^{T} s(t) \cdot sin(-\omega t) \ dt \\ \end{aligned} } $$

这就是傅里叶变换的奇函数表达式,也被称为 正弦傅里叶变换(SFT [Sine Fourier Transform])

同理,如果我们取偶函数,有 $$a_{\omega}$$$$b_{\omega}$$ 表现为符号相同的纯实数。即:

$$ {\displaystyle \begin{aligned} \hat{a}{\omega} &= 1 \cdot [\hat{s}(-\omega) + \hat{s}(\omega)] = 2 \cdot \hat{s}(\omega) \ \ \ \ \ \hat{b}{\omega} = \tfrac{1}{i} \cdot [\hat{s}(-\omega)-\hat{s}(\omega)] = 0 \ s(t) &= \frac{1}{N} \sum_{\omega =0}^{\omega_N} \ \ \ \ {2 \cdot \hat{s}(\omega)} \cdot cos(\omega t) \ \ \ \ \ + \ \ \ \ i \cdot 0 \cdot sin(\omega t) \ &= \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{1}{N}\sum_{n = 0}^{N} \hat{s}(\omega) \cdot cos(\omega t) \ \end{aligned} } $$

采用相同处理,有余 弦傅里叶变换(CFT [Cosine Fourier Transform]) 结果如下:

$$ {\displaystyle \begin{aligned} s(t) &= \frac{1}{N} \int_{0}^{N} \hat{s}(\omega) \cdot cos(\omega t) \ d{n} \\ \hat{s}(\omega) &= \int_{-\tfrac{T}{2}}^{+\tfrac{T}{2}} s(t) \cdot cos(-\omega t) \ dt \\ \end{aligned} } $$

然而工程中的信号并不存在有限周期且并不都能判定奇偶性,这是否意味着我们无法对其进行分解和化简?

答案是否定的。首先来看,针对周期性需要进行的操作。

解构一维信号 - 时频分离(Time-Frequency Separation)

如果我们换个角度就会发现,不存在有限周期只不过是因为周期太长,以至函数周期等于信号完整时长或着趋近无穷而导致的。所以我们分解原函数到对应的复指数函数和,所选择基底复指数函数也趋近于无穷,并使其对应频率从 $$0$$$$\infty$$ 而周期从极大到极小即可。不过在计算上就需要利用傅立叶变化的空间特征了。

结合上文,记被分解的原信号函数为 $$f(t)$$ 。根据傅立叶基的正交特性,如果存在 $${\mathcal {F}}(t)$$ 为当前 $$f(t)$$ 的解函数空间,则必然有 $$f(t) \cdot {\mathcal {F}}_{\omega}^{-1}(t)$$ 内积在时间 $$t$$ 范围为 $$(0,\ \infty)$$ 有固定值 $$\hat{f}(\omega)$$ ,使得:

$$ {\displaystyle \begin{aligned} \hat{f}(\omega) &= \int_{0}^{\infty} f(t) \cdot {\mathcal {F}}{\omega}^{-1}(t) \ dt = \int{0}^{\infty} f(t) \cdot e^{-i \omega t}\ dt \ \end{aligned} } $$

以函数空间角度排除 $$f(t)$$ 周期干扰。而复指数波的波函数,顾名思义就是复指数函数,有:

$$ {\displaystyle \begin{aligned} \hat{f}(\omega) &= \int_{-\infty}^{+\infty} a_{\omega} \cdot cos(\omega t) + i \cdot b_{\omega} \cdot sin(\omega t) \ dt\\ \end{aligned} } $$

使 $$b_{\omega}$$ 可取复数域,就可以转换为:

$$ {\displaystyle \begin{aligned} \hat{f}(\omega) &= \int_{-\infty}^{+\infty} a_{\omega} \cdot cos(\omega t) + b_{\omega} \cdot sin(\omega t) \ dt\\ \end{aligned} } $$

由于实际信号并不能严格确定奇偶性,不过对于小于四维的情况下,大多数条件都能保证其本身为实函数(即函数只有实数域取值),因而构成原信号的分离基底函数是存在不同强度和初项的。我们沿用前文中对初相和振幅的定义,记 $${\mathcal {F}}{\omega}(t)$$ 初相为 $$\angle\phi{\omega}$$ ,振幅为 $$A_{\omega}$$ ,则有:

$$ {\displaystyle \begin{aligned} {\mathcal {F}}{\omega}(t) : \quad \angle\phi{\omega} = \arctan(\tfrac{\hat{a}{\omega}}{\hat{b}{\omega}}) \ \ \ \ A_{\omega} = \sqrt{ (\hat{a}{\omega}) ^2 + (\hat{b}{\omega}) ^2 } \ \end{aligned} } $$

根据 帕西瓦尔定理(Parseval’s Theorem) 转复数空间,我们会发现 $$A_{\omega}$$ 就是 $$\hat{f}(\omega)$$$$2$$ 范数后的结果,而初项其实就是 $$\hat{f}(\omega)$$$$t = 0$$ 时,自身相位在复数空间上与实轴的夹角。即:

$$ {\displaystyle \begin{aligned} {\mathcal {F}}{\omega}(t) &: \ \angle\phi{\omega} &= \angle{\vert \hat{f}(t) \vert} \ = \arctan(\tfrac{\hat{a}{\omega}}{\hat{b}{\omega}}) \ A_{\omega} &=\ \ \Vert \hat{f}(t) \Vert 2 =\sqrt{ (\hat{a}{\omega}) ^2 + (\hat{b}_{\omega}) ^2 } \ \end{aligned} } $$

进而有:

$$ {\displaystyle \begin{aligned} {\mathcal {F}}{\omega}(t) &= A{\omega} \cdot sin(\omega t -\angle\phi_{\omega}) = A_{\omega} \cdot cos(\omega t +\angle\phi_{\omega}) \ &= {\Vert \hat{f}(t) \Vert 2} \cdot sin(\omega t -\angle{\vert \hat{f}(t) \vert}) = {\Vert \hat{f}(t) \Vert 2} \cdot cos(\omega t +\angle{\vert \hat{f}(t) \vert}) \ \hat{f}(\omega) &= \int{0}^{\infty} f(t) \cdot e^{-i \omega t}\ dt \ \ \ \ \ \Leftrightarrow \ \ \ \ \ f(t) = \frac{1}{N} \int{-\infty}^{+\infty} \hat{f}(\omega) \cdot {\mathcal {F}}_{\omega}(t) \ d \omega \ \end{aligned} } $$

显然,大部分信号都是有限时间下的,且基本都能满足无穷区间的狄利克雷条件,也因此可以使用傅里叶变换分解。

如果频率范围在 $$\omega \in [\omega_{0},\ \omega_{1}]$$ ,对于选定的时间点 $$t = t_c$$ ,有频率 $$\omega$$ 、原函数 $$f(t)$$$$t = t_c$$ 时的取值 $$f(t_c)$$ 、基底函数族 $${\mathcal {F}}{\omega}(t)$$ 锁定时间 $$t = t_c$$ 的变体 $${\mathcal {F}}{t_c}(\omega)$$ ,构成该频率范围的 频域投影(FDP [Frequency Domain Projection])

反之,如果时间范围在 $$t\in [\ t_0,\ \ t_1]$$ ,对于频率范围 $$\omega \in [\omega_{0},\ \omega_{1}]$$ ,有时间 $$t$$ 、原函数 $$f(t)$$ 、基底函数族 $${\mathcal {F}}_{\omega}(t) $$,就构成了原函数在该时间范围的 时域投影(TDP [Time Domain Projection])

两者的区别仅在于观察角度的不同:

$$ {\displaystyle \begin{aligned} {Frequency\ Domain\ Projection:} &\ \ (\ \ \omega\ ,\ \ f(t_c)\ ,\ \ {\mathcal {F}}{t_c}(\omega) \ \ ) \ {Time\ Domain\ Projection:} &\ \ (\ \ t\ \ ,\ \ f(t)\ \ ,\ \ {\mathcal {F}}{\omega}(t) \ \ \ ) \ {\omega \in [\omega_0,\ \omega_n]} \ \ \ \ & \ \ {\ t\ \in [\ t_0,\ \ t_n\ ]} \ \end{aligned} } $$

周期的问题解决了,现在我们能够拿到时频分离(Time-Frequency Separation)的原信号函数信息并可以依此还原信号本身。但积分对于计算机来说任务有些繁重。同时,由于计算机只能处理离散化后的数字信号,因此离散化的算法才能够被计算机有效使用。

所以还需要在此基础上,找到更为便捷的算法实现。

精简运算过程 - 一维离散傅立叶变换(1D-DFT)

如果将积分重新转换为级数形式积化和差表示,并在允许误差范围内取有限子集。那么就能够化解掉大部分运算量,从而得到一个相对理论而言的低时间复杂度算法。这种想法促成了基于计算机运算的一维离散傅立叶(1D-DFT)的诞生。

一维离散傅立叶(1D-DFT [1D-Discrete Fourier Transform])本质上包含了两部分离散化作业,即对时域的离散化(TDD [Time Domain Discrete])和对频域的离散化(FDD [Frequency Domain Discrete])。

时域离散化(TDD) 方面,一维离散傅立叶采用了离散时间傅立叶变化(DTFT [Discrete Time Fourier Transform])中,对时域信号间隔采样的操作。即将:

$$ {\displaystyle \begin{aligned} \hat{f}(\omega) &= \int_{0}^{\infty} f(t) \cdot e^{-i \omega t}\ dt \\ \end{aligned} } $$

以时间采样(切片)数量为 $${n_1}$$ ,转为级数形式:

$$ {\displaystyle \begin{aligned} \hat{f}(\omega) &= \sum_{t = t_0}^{t_{n_1}} f(t) \cdot e^{-i \omega t} \\ \end{aligned} } $$

打破时间上的连续性。需要注意的是,此时频域仍然是连续的。

频域离散化(FDD) 方面,离散傅立叶做的操作就更为直观了。如果在频率采样时就以离散化的方式采样数据,那得到的频域信息天然就是离散的。同样,从某个时刻 $$t = t_c$$ 离散化的频域信息上还原当前实际频率,则也是一个线性求和的过程。因此有:

$$ {\displaystyle \begin{aligned} f(t) = \frac{1}{N} \int_{-\infty}^{+\infty} \hat{f}(\omega) \cdot {\mathcal {F}}_{\omega}(t) \ d \omega \\ \end{aligned} } $$

以频率采样(切片)数量为 $${n_2}$$ ,转为级数形式:

$$ {\displaystyle \begin{aligned} f(t) = \frac{1}{n_2} \sum_{\omega = \omega_0}^{\omega_{n_2}} \hat{f}(\omega) \cdot {\mathcal {F}}_{\omega}(t) \\ \end{aligned} } $$

而随着有限采样,基底函数族 $${\mathcal {F}}_{\omega}(t) $$$ 构成的解函数空间也是有限维的,即:

$$ {\displaystyle \begin{aligned} {\mathcal {F}}{\omega} = [{\mathcal {F}}{\omega_1},{\mathcal {F}}{\omega_2},\ ...\ ,{\mathcal {F}}{\omega_{n_2}}] \ \end{aligned} } $$

至此,由时域离散化(TDD)与频域离散化(FDD)共同构成离散傅立叶(DFT)的完整表达如下所示:

$$ {\displaystyle \begin{aligned} {\mathcal {F}}{\omega} = [{\mathcal {F}}{\omega_1},&{\mathcal {F}}{\omega_2},\ ...\ ,{\mathcal {F}}{\omega_{n_2}}] \ \hat{f}(\omega) = \sum_{t = t_0}^{t_{n_1}} f(t) \cdot e^{-i \omega t} \ \ \ \ \ &\Leftrightarrow \ \ \ \ \ f(t) = \frac{1}{n_2} \sum_{\omega = \omega_0}^{\omega_{n_2}} \hat{f}(\omega) \cdot {\mathcal {F}}_{\omega}(t) \ \end{aligned} } $$

经过离散化后的有限采样更适合计算机有限的算力,因此才能被程序化。不过由于并没有保存连续的完整信息,经过离散傅里叶变换后再还原的数据,相对于采样自然源的原始数据终归还是会有一定损失的。但是由于变换与逆变换,并不会导致解构再还原后的数据存在差异。所以离散傅里叶变换被归类为 有损采样(Lossy Sampling)的无损算法(Lossless Algorithm)

一维离散傅立叶变换(1D-DFT)的 C 语言实现

既然需要做程序化,那么首先需要将离散傅里叶变换的过程抽象化。理清逻辑思路的同时,方便构造迭代器和代码的处理流水线。这里我们使用伪码表示:

/**
 * 1D-DFT [Discrete Fourier Transform]
 * [How to Use]
 *   <case:>
 *     Fo[T] = {...};
 *     Fn[N] = {};
 *     dft_1d(&Fo, &Fn, T, N);
 * [theorem::definitions]
 *   Fo meanings Original Function
 *   Fn meanings Fourier Basis at [n]
 *   pi meanings π
 *   T  meanings Periodic of Fo
 *   N  meanings Slice of Frequency
 *   Wn meanings Angular Frequency of Basis Fn is Wn = ((2*pi*n)/T)
 * [theorem::formula]
 *   Fo[t] = sum_{n=0}^{N-1} x Fn[t] * exp(  i * ((2*pi*n)/T) * t), 0<=n<N
 *   Fn[t] = sum_{t=0}^{T-1} x Fo[t] * exp( -i * ((2*pi*n)/T) * t), 0<=t<T
 *         = sum_{t=0}^{T-1} x Fo[t] * (An * cos(Wn * t) + Bn * sin(Wn * t)), 0<=t<T
 *         = All_in_one(An) * cos(Wn * t) + All_in_one(Bn) * sin(Wn * t), 0<=t<T
 *   So:
 *     An = Re(Fn[t]);  Bn = Im(Fn[t])
 * [logistic]
 * {
 *   result = []; // as byte array
 *   // do TDD:
 *   for n in range(N) {
 *     An = 0; Bn = 0;
 *     // do FDD:
 *     for t in Range(T) {
 *       Wn = (2 * pi * n) / T;
 *       An = Re += Cos(Wn * t) * Fo(t);
 *       Bn = Im += Sin(Wn * t) * Fo(t);
 *     }
 *     result[n] = Fn.to_complex_angular_form(An, Bn)
 *   }
 *   return result;
 * }
 * @param Fo Original Function input array
 *           (already sampled by T as count-of-time-slice)
 * @param Fn Fourier Basis
 *           (already sampled by N as count-of-frequency-slice)
 * @param T  Periodic of Fo, Number of Time Slice
 * @param N  Number of Frequency Slice
 */

同时,我们还需要提供离散傅里叶变换的逆变换(IDFT [Inverse Discrete Fourier Transform])来使得电脑能够还原信息:

/**
 * 1D-IDFT [Inverse Discrete Fourier Transform]
 * [How to Use]
 *   <case:>
 *     Fo[T] = {};
 *     Fn[N] = {...};
 *     dft_1d(&Fo, &Fn, T, N);
 * [theorem::definitions]
 *   Fo meanings Original Function
 *   Fn meanings Fourier Basis at [n]
 *   pi meanings π
 *   T  meanings Periodic of Fo
 *   N  meanings Slice of Frequency
 *   Wn meanings Angular Frequency of Basis Fn is Wn = ((2*pi*n)/T)
 * [theorem::formula]
 *   Fo[t] = sum_{n=0}^{N-1} x Fn[t], 0<=n<N
 *   Fn[t] = sum_{t=0}^{T-1} x Fo[t] * exp( -i * ((2*pi*n)/T) * t), 0<=t<T
 *         = sum_{t=0}^{T-1} x Fo[t] * (An * cos(Wn * t) + Bn * sin(Wn * t)), 0<=t<T
 *         = All_in_one(An) * cos(Wn * t) + All_in_one(Bn) * sin(Wn * t), 0<=t<T
 *   So:
 *     An = Re(Fn[t]);  Bn = Im(Fn[t])
 * [logistic]
 * {
 *   result = []; // as byte array
 *   // do TDD:
 *   for t in range(T) {
 *     Re = 0; Im = 0;
 *     // do FDD:
 *     for n in Range(N) {
 *       Wn = (2 * pi * n) / T;
 *       An = Re(Fn[n]);
 *       Bn = Im(Fn[n]);
 *       result[t] += Fn[n].to_value(Wn, An, Bn);
 *     }
 *   }
 *   return result;
 * }
 * @param Fo Original Function input array
 *           (already sampled by T as count-of-time-slice)
 * @param Fn Fourier Basis
 *           (already sampled by N as count-of-frequency-slice)
 * @param T  Periodic of Fo, Number of Time Slice
 * @param N  Number of Frequency Slice
 */

现在思路有了,只需要以代码实现即可:

#include "stdio.h"
#include "math.h"

#define PI 3.1415926f

typedef struct FBasis {
    double re_;
    double im_;
    double w_;
} FBasis;

void dft_1d(double *Fo, FBasis *Fn, size_t T, size_t N) {
    for (int n = 0; n < N; ++n) {
        double An = 0;
        double Bn = 0;
        double Wn = (2 * PI * n) / T;
        for (int t = 0; t < T; ++t) {
            An += cos(Wn * t) * Fo[t];
            Bn += sin(Wn * t) * Fo[t];
        }
        FBasis e_ = {An, Bn, Wn};
        Fn[n] = e_;
    }
}

void idft_1d(double *Fo, FBasis *Fn, size_t T, size_t N) {
    for (int t = 0; t < T; ++t) {
        for (int n = 0; n < N; ++n) {
            FBasis e_ = Fn[n];
            Fo[t] += (e_.re_ * cos(e_.w_ * t) + e_.im_ * sin(e_.w_ * t)) / N;
        }
    }
}

写完后简单测试一下:

int main(void) {
    FBasis Fn[6] = {};
    double Fo[6] = {1, 2, 3, 4, 5, 6};
    double iFo[6] = {};
    size_t T = sizeof(Fo) / sizeof(double);
    size_t N = sizeof(Fn) / sizeof(FBasis);
    printf("\n Original_data: \n");
    for (int t = 0; t < T; ++t) {
        printf("%f  ", Fo[t]);
    }

    printf("\n DFT_result: \n");
    dft_1d(Fo, Fn, T, N);
    for (int n = 0; n < N; ++n) {
        printf("%f + i %f \n", Fn[n].re_, Fn[n].im_);
    }

    printf("\n IDFT_result: \n");
    idft_1d(iFo, Fn, T, N);
    for (int t = 0; t < T; ++t) {
        printf("%f  ", iFo[t]);
    }

    return 0;
}

得到结果和标准几近相同:

 Original data: 
1.000000  2.000000  3.000000  4.000000  5.000000  6.000000  

 DFT_result: 
21.000000 + i 0.000000 
-3.000003 + i -5.196152 
-3.000002 + i -1.732048 
-3.000000 + i -0.000002 
-2.999996 + i 1.732057 
-2.999979 + i 5.196158 

 IDFT_result: 
1.000003  2.000000  2.999999  3.999999  4.999999  6.000000   

运行结束。

到这里,我们已经基本掌握了傅里叶变换原理和最基础的应用。

如果拓展傅里叶变换到相对复杂的二维情况,那么和一维时有哪些不同呢?