离散傅里叶变换DFT详解及应用

时间：2023-07-27

　　(本人为本文原作者,转载请标明出处)

　　上文请转:

　　DBinary：傅里叶变换推导详解其中,式3.51及3.52被称为傅里叶变换对,其中3.52为正变换,表示离散时间信号变换为频域复信号,式3.51位逆变换,表示由频域复信号还原为离散时间信号,在本文之后的论述中, 将有限长度且离散的有限长离散时域傅里叶变换的分析公式统一简称为离散傅里叶变换（Discrete Fourier Transform）或者其英文缩写DFT，将它的逆变换（Inverse Discrete Fourier Transform）也就是逆变换公式简称为IDFT[3]。

　　自然信号若要保留为数字离散信号，首先要做的就是对该信号进行采样，一般情况下，使用的采样方式都是等间距采样，即为样本两两之间的采样间隔相同，在数学表达上，常常使用单位冲击采样来表示这一过程。现定义单位冲击函数

　　 $\delta(t)$ ,这个函数仅在t=0时有值,其它区间的取值都是0,该函数的的定义是

　　在实际的使用中,常常定义

　　假设一连续时间信号函数 $f(t)$ ,需要对t=2时间进行采样,仅需要对 $\delta(t)$ 进行时移就可以了,定义一个连续时间信号的采样样本x[n]为式4.2

　　 $x\left\lbrack n \right\rbrack = f\left( t \right)\delta\left( t - n \right)\ \ \ \ \ (4.2)$

　　因此假设我们需要对一个信号多点进行采样,可以使用式4.2的方式进行表达:

　　 $x\left\lbrack n \right\rbrack = \sum_{n = - \infty}^{+ \infty}{f\left( t \right)\delta(t - n)}\text{}(4.2)$

　　在《信号与系统》[2] 第7.1.1章节推论如下,设对一信号以T间隔取样，那么就有如下数学式

　　 $x_{p}\left( t \right) = \sum_{n = - \infty}^{+ \infty}{x(nT)\delta(t - nT)}\text{ }(4.3)$

　　根据卷积性质时域内的相乘对应于频域内的卷积，其中 $\ P(j\omega)$ 为冲击串频域函数

　　 $X_{p}\left( \text{jω} \right) = \frac{1}{2\pi}\int_{- \infty}^{+ \infty}{x(j\omega)P(j(\omega - \theta))d\theta}\ \ (4.4)$

　　根据公式4.7,4.8:

　　 $P\left( \text{jω} \right) = \frac{2\pi}{T}\sum_{n = - \infty}^{+ \infty}{\delta(\omega - k\omega_{s})}\text{ }(4.5)$

　　 $X_{p}\left( \text{jω} \right)*{\ δ}\left( \omega - \omega_{0} \right) = X\left( j\left( \omega - \omega_{0} \right) \right)\text{ }(4.6)$

　　可得出:

　　 $X_{p}\left( \text{jω} \right) = \frac{1}{T}\sum_{n = - \infty}^{+ \infty}{X(j(\omega - k\omega_{s}))}(4.7)$

　　其中 $\omega_{s}$ 为基波频率或称之为频域间距,上式可知 $X_{p}\left( \text{jω} \right)$ 是频率为 $\omega$ 的周期函数，由一组移位的 $X\left( \text{jω} \right)$ 叠加而成，在幅度标以 $\frac{1}{T}$ 的变化，当 $\omega_{s}$ 大于频带宽度的一半时频带不会发生堆叠，反之会发生堆叠。

　　因此可以得出一个结论，信号的采样频率至少需要该信号最高频率的2倍，否则无法还原出原信号导致失真。

　　同时依据第二章第六节的推论，有DFT正变换公式如下:

　　 $X_{k} = \sum_{n = 0}^{N - 1}{x_{n}e^{\frac{- 2\pi j}{N}\text{kn}}}\text{}(4.8)$

　　其中，N表示信号的采样点数，n的范围是0到 $N-1，x_{0}$ 表示信号的第一个点， $x_{1}$ 表示信号的第二个点，依此类推，k的取值范围是0到N-1，由 $f = \frac{\omega}{2\pi}$ 可知，其分辨率等于:

　　 $\left( \frac{k + 1}{N} - \frac{k}{N} \right)f_{s} = \frac{f_{s}}{N}\ \ (4.9)$

　　因此得出另一个结论,当对一个信号以频率 $f_{s}$ 进行采样时,对该信号的傅里叶变换后的频域分辨率是采样率除以其样本点数。

　　若某一函数在坐标系上以x=N轴对称，那么其可以写成式4.10的形式

　　 $f\left( x \right) = f\left( N - x \right)\ \ \ \ (4.10)$

　　如图4-2-1，其中的函数 $f\left( x \right)及g\left( x \right)$ 都是以x=0的轴线对称。

　　图4-2-1 两个轴对称函数图像

　　在复数中，两个实部相等，虚部互为相反数的复数互为共轭复数(conjugate complex number)，假设一复数z，则其共轭写作 $\overset{\overline{}}{z}$ ,其定义满足若

　　 $z = a + ib\ \ \ \ (4.11)$

　　则其共轭为

　　 $\overset{\overline{}}{z} = a - \text{ib }(4.12)$

　　现定义一个离散序列x[n]若x[n]满足式4.13

　　 $x\left\lbrack n \right\rbrack = \overset{\overline{}}{x}\left\lbrack - n \right\rbrack\ \ (4.13)$

　　则称x[n]为共轭对称序列。

　　假设：

　　 $W_{N} = e^{- \frac{2\text{πi}}{N}}\text{ }(4.14)$

　　依据欧拉公式，可以得出

　　 $W_{N} = \cos\left( \frac{2\pi}{N} \right) - isin\left( \frac{2\pi}{N} \right)\$

　　那么，上式可根据4.8写成

　　 $X_{k} = \sum_{n = 0}^{N - 1}{x_{n}W_{N}^{\text{kn}}}\text{}(4.15)$

　　设 $k = N - k$ ，并将它带入 $W_{N}^{\text{kn}}$ 中，那么就有

　　 $W_{N}^{(N - k)n} = e^{- i\frac{2\pi}{N}(N - k)n} = \cos\left( \frac{2\pi}{N}(N - k)n \right) - isin\left( \frac{2\pi}{N}\left( N - k \right)n \right) = \cos\left( \frac{2\pi}{N}\text{kn} \right) + isin\left( \frac{2\pi}{N}\text{kn} \right)\text{ }(4.16)$

　　 $W_{N}^{\text{kn}} = \cos\left( \frac{2\pi}{N}\text{kn} \right) - isin\left( \frac{2\pi}{N}\text{kn} \right) = {\overset{\overline{}}{W}}_{N}^{(N - k)n}\text{ }(4.17)$

　　由式4.17可知,傅里叶变换结果在 $k \in \lbrack 1,N - 1\rbrack$ 范围内呈共轭对称,而当k=0时有式4.18

　　 $X_{0} = \sum_{n = 0}^{N - 1}{x_{n}\cos\left( 0 \right)} - i\sum_{n = 0}^{N - 1}{x_{n}\sin\left( 0 \right)} = \sum_{n = 0}^{N - 1}x_{n}$

　　也即说明变换的第一个结果为 $x_{n}$ 的累加和,即频率为0部分的累加和,在变换结果中称之为信号的直流分量。

　　离散傅里叶的正变换的结果为离散复数序列,设某一结果为

　　 $a,bi$

　　在实际的使用中,常常使用 $\sqrt{a^{2} + b^{2}}$ 表示对应频率的能量,该能量并非物理意义上的能量,它用来表示该频率信号对原始信号的贡献,同时使用 $\frac{b}{a}$ 来表示对应频率的相位,描述的是该频率的位置信息。

　　DFT在声学的分析中被广泛的使用，WavFreq是一款能够加载Wav音乐（波形文件）并进行解析的软件。下图显示的是该软件对一段音乐进行分析的频谱图：

　　图4-3-1 音乐频谱图像其中，该音频信号的采样频率为44100HZ，每个样本以16位短整形记录，这也就意味着声波信号被量化在范围[-32768,32767]内,从图中的频谱图中可以看出,大部分能量较高的频率被集中在了低频当中,在正常的音乐中,人声大多被集中在了低频范围。女声的频率也会相对的比男声的频率高出一些。

　　第一节 C语言与环境

　　本文所述的程序部分，全部都使用C语言来完成，同时，本文使用开发环境为windows 10 64位系统，集成开发环境使用微软官方发布的Visual Studio 2010（编译器MSVC 2010），C语言依据C99标准编写。因此，程序运行需要安装VC++2010的支持库。

　　为了保证代码的可移植性，程序部分仅使用C语言提供的数学标准库math.h,该标准库在众多的编译器中都被支持。

　　第二节使用C语言编写DFT代码

　　因为涉及到复数运算，因此，本文定义一个结构体来描述一个复数

　　其中该结构体名为complex，当中有两个float类型的成员分别代表复数的实部和虚部，其中re意味really，表示复数的实部，im表示Imaginary表示虚部。同时，实现三个函数用于对虚部进行运算。

　　其中 complex complexBuild(float re,float im);用于构建并初始化一个复数结构体,re表示构建复数的实部,im表示构建复数的虚部

　　complex complexAdd(complex a,complex b)用于实现两个复数的相加运算,即虚部加虚部,实部加实部.

　　complex complexMult(complex a,complex b)用于实现两个复数的乘法运算.依据公式

　　 $\left( a_{1} + ib_{1} \right)\left( a_{2} + ib_{2} \right) = a_{1}a_{2} + ia_{1}b_{2} + ia_{2}b_{1} - b_{1}b_{2}$

　　根据第一章节第六节的推导公式3.52

　　 $X\left\lbrack n \right\rbrack = \sum_{t = 0}^{T - 1}{x\left\lbrack t \right\rbrack e^{- jn\frac{2\pi}{T}t}}\text{}(3.52)$

　　根据欧拉公式上式可以写成

　　 $X\left\lbrack k \right\rbrack = \sum_{n = 0}^{N - 1}{x\left\lbrack n \right\rbrack}\left( \cos\left( \frac{2\text{πkn}}{N} \right) - \text{jsin}\left( \frac{2\text{πkn}}{N} \right) \right)\ (5.1)$

　　那么,DFT的公式就可以写成

　　其中x[]表示输入的复信号序列,实信号可以当做虚数为0的复数信号进行输入,X[]表示输出的复数序列,其分配的内存大小至少应该为输入序列的大小,N为输入序列的复信号个数。

　　在项目 DFT_DEMO中,对一个长度为10的序列进行DFT变换,其中,该序列为1,2,3,4,5,6,7,8,9,10十个数,如图4-2-1

　　图4-2-1 DFT变换的结果

　　代码如下

　　将上图结果与matlab中的DFT变换结果对比计算结果一致,可知代码编写正确。

　　第三节使用C语言编写IDFT代码

　　根据第六节傅里叶变换对的逆变换公式3.46

　　 $x\left\lbrack t \right\rbrack = \frac{1}{2\pi}\int_{0}^{2\pi}{X\left( e^{\text{jω}} \right)e^{\text{jωnt}}\text{dω}}\text{ }(3.46)$

　　可得到适合代码编写的公式5.2

　　 $x\left\lbrack k \right\rbrack = \frac{1}{N}\sum_{n = 0}^{N - 1}{X\left\lbrack n \right\rbrack}\left( \cos\left( \frac{2\text{πkn}}{N} \right) + i\sin\left( \frac{2\text{πkn}}{N} \right) \right)\ (5.2)$

　　通过第三节的准备工作,可以直接套用第三节给出的结构体声明,直接在源代码当中写出逆变换的代码。

　　定义IDFT变换函数为

　　实现代码如下所示

　　其中X[]表示输入的复信号序列,x[]表示输出的信号序列,N表示输入复信号序列的长度，如果逆变换编写正确，那么将第二节的输出作为逆变换的输入，就可以得到正变换前的离散信号序列，为了验证这一结果，将main函数更改如下：

　　运行结果如图4-3-1

　　图4-3-1 傅里叶变换对程序执行结果从结果来看IDFT的程序成功将DFT变换的结果还原成原信号，逆变换的结果有极小的偏差，这是因为计算过程中浮点运算存在极小的精度损失，导致变换的结果有极小的误差。