注：本文为《The Scientist and Engineer‘s Guide to Digital Signal Processing》的读书笔记系列之第八章。
Chapter 8. The Discrete Fourier Transform
- 1. The Family of Fourier Transform
- 2. Notation and Format of the Real DFT
- 3. The Frequency Domain's Independent Variable
- 4. DFT Basis Functions
- 5. Synthesis, Calculating the Inverse DFT
本章主要讲述实数表示的 DFT
1. The Family of Fourier Transform
Jean Baptiste Joseph Fourier对热的传播很感兴趣，他发表的一篇论文中提出了一个有争议的观点：任何连续周期信号可以表示为一系列正弦信号的和。
图1展示了如何将一个信号分解成正弦信号和余弦信号。
图 1-a 时域信号
图 1-b 1-a 被分解后得到的余弦和正弦信号
图1-a展示了16点的时域信号，图1-b展示了对该时域信号进行傅立叶分解得到的9个余弦信号和9个正弦信号，每一个都具有不同的频率和幅度。
为什么使用正弦波而不是方波或者三角波？其实，有很多种方式都能对信号进行分解，分解的目标是得到比原始信号更容易处理的东西。正弦波和余弦波比原始信号更简单，因为它们具有原始信号不具有的特性：sinusoidal fidelity， 如第五章所讨论的，系统对一个正弦波的响应还是一个正弦波。只有信号的幅度和相位可以改变，频率和波形必须保持一致。正弦波是唯一具有此属性的波形。
通用术语Fourier transform可以被细分为四种类型，这是由于有四种不同类型的信号。信号可以被分为连续的和离散的，也可以被分为周期的和非周期的，总共会产生四种不同类型的信号，如图2所示。
图 2 四种不同类型的信号及其对应的变换名称
连续非周期
这种信号的Fourier transform被简单地称为the Fourier transform。
连续周期
这种信号的Fourier transform被称为the Fourier Series。
离散非周期
这种信号的Fourier transform被称为the Discrete Time Fourier transform。
离散周期
这种信号的Fourier transform有时被称为the Discrete Fourier Series，但最常被称为the Discrete Fourier transform。
这四种类型的信号均被设定是从负无穷延伸到正无穷，即无限长。假如有一个1024点的有限长信号，那么是不是就没有适合其Fourier transform的版本？正弦波和余弦波都是无限长的信号，不能使用一组无限长的信号合成有限长的信号。解决这个难题的方法是让有限长信号look like无限长信号。这是通过想象在有限长信号的左右两侧具有无限采样点来完成的。如果这些想象出来的无限采样点全部是零，这时候信号看起来是离散非周期的，可以对其应用the Discrete Time Fourier transform；也可以将想象出来的采样点认为是真实的1024个采样点的拷贝。这时候信号看起来是离散周期的，可以对其应用the Discrete Fourier transform。
事实证明，合成一个非周期的信号需要无数个正弦波。计算机算法是不可能计算the Discrete Time Fourier transform的。可以在DSP中使用的唯一的Fourier transform类型是DFT。也就是说，数字计算机只能使用离散且长度有限的信息。当在理论推导或者家庭作业中思考一些数学问题时，你可能会发现你使用的是Fourier transform家族中的前三种类型，而当在使用计算机时，将仅使用DFT，本章主要讲解该中类型的变换。
回顾图1的例子，表面上看，似乎是一个16点的信号被分解成为18个正弦波，每个正弦波包含16个点。用更正式的术语说，图1-a中所示的16点的信号必须被认为是无限长周期信号的某一周期。同样的，图1-b中所示的18个正弦信号中的每个都代表无限长正弦曲线中的16个点。如果将其视为由16个点正弦波合成的16个点的信号，或者视为由无限长的正弦波合成的无限长的周期信号，这真的有关系吗？答案是：usually no, but sometimes, yes。
数学术语transform经常在数字信号处理中出现，那么什么是transform呢？回忆一下function，function是指将一个数或者一些数变换为另一个数的算法或者过程。transform是对function的扩展，它允许输入输出都是一些数。a transform is any fixed procedure that changes one chunk of data into another chunk of data。
2. Notation and Format of the Real DFT
如图3所示，离散傅立叶变换将一个 N N N点的输入信号变换成两个 N / 2 + 1 N/2+1 N/2+1点的输出信号。这两个输出信号包含了正弦波和余弦波的幅度。输入的信号一般称为时域，输出的信号一般称为频域。术语频域被用来描述正弦波和余弦波的幅度。
图 3 DTF 示例
频域包含了和时域完全一样的信息，只不过是用不同的形式表示出来。标准的DSP表示法使用小写字母表示时域信号，如： x [ ] , y [ ] x[],y[] x[],y[]和 z [ ] z[] z[]。其对应的频域使用大写字母表示，如 X [ ] , Y [ ] X[],Y[] X[],Y[]和 Z [ ] Z[] Z[]。为了说明，假设 N N N点的时域信号表示为 x [ ] x[] x[]，其频域表示为 X [ ] X[] X[]，包含了实部和虚部两部分，每部分有 N / 2 + 1 N/2+1 N/2+1个点。 X [ ] X[] X[]的实部写作 R e X [ ] Re X[] ReX[]，虚部写作 I m X [ ] Im X[] ImX[]，其实部包含了余弦波的幅度，虚部包含了正弦波的幅度。
3. The Frequency Domain’s Independent Variable
图4演示了 N = 128 N=128 N=128点的DFT。
图 4-a 时域信号
图 4-b、 图 4-c 频域信号
频域的横轴可以有四种不同的解释。第一种，横轴被标记为从0到64，使用字母 k ∈ [ 0 , N / 2 ] k\in[0, N/2] k∈[0,N/2]表示该索引，程序员最喜欢使用这种方式。图4-b采用的这种方式。第二种，图4-c所采用的方式，横轴被标记为fraction of the sampling rate，这表示横轴的取值从0到0.5，因为离散数据只能包含介于DC和一半采样率之间的频率。使用字母 f f f表示该索引，且 f = ( k / N ) ∈ [ 0 , 0.5 ] f=(k/N)\in[0, 0.5] f=(k/N)∈[0,0.5]。第三种，和第二种类似，只不过横轴上乘了 2 π 2\pi 2π。使用字母 ω = ( 2 π f ) ∈ [ 0 , π ] \omega=(2\pi f)\in[0, \pi] ω=(2πf)∈[0,π]表示该索引，其表示natural frequency，单位为弧度。第四种，使用模拟频率Hz表示。
4. DFT Basis Functions
DFT中所使用的正弦波和余弦波通常称为DFT的基函数。换句话说，DFT的输出是一组表示幅度的数字。基函数是一组具有单位幅度的余弦波和正弦波。DFT的基函数如下：
c k [ i ] = c o s ( 2 π k i / N ) c_k[i]=cos(2\pi ki/N) ck[i]=cos(2πki/N)
s k [ i ] = s i n ( 2 π k i / N ) (1) s_k[i]=sin(2\pi ki/N)\tag{1} sk[i]=sin(2πki/N)(1)
其中， c k [ ] c_k[] ck[]是余弦波，其幅度存储在 R e X [ k ] ReX[k] ReX[k]中， s k [ ] s_k[] sk[]是正弦波，其幅度存储在 I m X [ k ] ImX[k] ImX[k]中。 k k k决定了其频率。图5展示了 N = 32 N=32 N=32点的DFT使用17个正弦波和17个余弦波作为其基函数。
图 5 DFT 的基函数
由于这些信号加起来共同构成了输入信号，因此它们必须具有和输入信号同样的长度。 k k k决定了每个信号的频率。比如， c 1 [ ] c_1[] c1[]是在N点上完成一个周期的余弦波， c 5 [ ] c_5[] c5[]是在 N N N点上完成五个周期的余弦波等等。这是理解基函数的一个重要概念，参数 k k k等于在 N N N点上完成的整周期数。
图5-a展示了频率为0的余弦波 c 0 [ ] c_0[] c0[]，其值全为1。这表示 R e X [ 0 ] ReX[0] ReX[0]存储时域信号所有点的平均值，也称 R e X [ 0 ] ReX[0] ReX[0]存储DC offset(直流偏移)。图5-b展示了频率为0的正弦波 s 0 [ ] s_0[] s0[]，其值全为0。 由于其不会影响时域信号的合成，因此 I m X [ 0 ] ImX [0] ImX[0]的值无关紧要，并且始终设置为零。
图5-c和图5-d分别展示了 c 2 [ ] c_2[] c2[]和 s 2 [ ] s_2[] s2[]，它两均在 N N N点上完成了两个整周期，分别对应于 R e X [ 2 ] ReX[2] ReX[2]和 I m X [ 2 ] ImX[2] ImX[2]。图5-e和图5-f与之类似，不过要是图中的离散点之间没有曲线链接，肉眼将很难看出它们是正弦波和余弦波。
图5-g和图5-h分别展示了最高频率的基函数 c 16 [ ] c_{16}[] c16[]和 s 16 [ ] s_{16}[] s16[]。离散余弦波的值在1和-1两者之间交替，这可以解释为在峰值处采样连续余弦波。相比之下，离散正弦波的值全为0，这是由于在连续正弦波的零点处进行采样。
有个问题：如果对 N N N个点进行DFT，而输出了 N + 2 N+2 N+2个点，那么额外的信息从何而来？ 答案是：输出中的两个点不包含任何信息，从而使其他 N N N个点完全独立。不携带任何信息的点是 I m X [ 0 ] ImX [0] ImX[0]和 I m X [ N / 2 ] ImX [N / 2] ImX[N/2]，其值始终为零。
5. Synthesis, Calculating the Inverse DFT
可以将合成方程写为：
x [ i ] = ∑ k = 0 N / 2 R e X ‾ [ k ] c o s ( 2 π k i / N ) + ∑ k = 0 N / 2 I m X ‾ [ k ] s i n ( 2 π k i / N ) (2) x[i]=\sum_{k=0}^{N/2}Re\overline{X}[k]cos(2\pi ki/N)+\sum_{k=0}^{N/2}Im\overline{X}[k]sin(2\pi ki/N)\tag{2} x[i]=k=0∑N/2ReX[k]cos(2πki/N)+k=0∑N/2ImX[k]sin(2πki/N)(2)
R e X ‾ [ k ] = R e X [ k ] N / 2 Re\overline{X}[k]=\frac{ReX[k]}{N/2} ReX[k]=N/2ReX[k]
I m X ‾ [ k ] = − I m X [ k ] N / 2 Im\overline{X}[k]=-\frac{ImX[k]}{N/2} ImX[k]=−N/2ImX[k]
e x p e c t f o r t w o s p e c i a l c a s e s : expect \, for \, two \, special \, cases: expectfortwospecialcases:
R e X ‾ [ 0 ] = R e X [ 0 ] N Re\overline{X}[0]=\frac{ReX[0]}{N} ReX[0]=NReX[0]
R e X ‾ [ N / 2 ] = R e X [ N / 2 ] N (3) Re\overline{X}[N/2]=\frac{ReX[N/2]}{N}\tag{3} ReX[N/2]=NReX[N/2](3)
其中， N N N是时域信号的采样点数， x [ i ] x[i] x[i]是合成的信号，且 i ∈ [ 0 , N − 1 ] i\in [0,N-1] i∈[0,N−1]。 R e X ‾ [ k ] Re\overline{X}[k] ReX[k]和 I m X ‾ [ k ] Im\overline{X}[k] ImX[k]分别表示合成时域信号所需的余弦波和正弦波的幅度， R e X [ k ] Re{X}[k] ReX[k]和 I m X [ k ] Im{X}[k] ImX[k]分别表频域的实部和虚部，记住它们的关系！！！ k ∈ [ 0 , N / 2 ] k\in [0,N/2] k∈[0,N/2]。式 2 2 2和 3 3 3一起定义了逆DFT。
表1展示了两种实现逆DFT的编码方法。
100 /*THE INVERSE DISCRETE FOURIER TRANSFORM
110 The time domain signal, held in XX[], is calculated from the frequency domain signals,
120 held i REX[] and IMX[].*/
130
140 float XX[512]; // XX[] holds the time domain signal
150 float REX[257]; // REX[] holds the real part of the frequency domain
160 float IMX[257]; // IMX[] holds the imaginary part of the frequency domain
170
180 const float PI = 3.14159265; // Set the constant, PI
190 const int N = 512 // N is the number of points in XX[]
200
210 // Initing REX[] and IMX[]
220
230
240 // Find the cosine and sine wave amplitudes using Eq. 3
250 for (int k=0; k <257; k++){
260 REX[k] = REX[k] / (N/2);
270 IMX[k] = -IMX[k] / (N/2);
280 }
290
300 REX[0] = REX[0] / 2;
310 REX[256] = REX[256] / 2;
320
330
340 for (int i=0; i<512; i++){ // Zero XX[] so it can be used as an accumulator
350 XX[i] = 0;
360 }
370
380 /*Method 1: Loop through each frequency generating the entire length of the sine and cosine waves, and add them to the accumulator signal, XX[] */
390 for (int k=0; k<257; k++){
400 for (int i=0; i<512; i++){
410 XX[i] = XX[i] + REX[k] * cos(2*PI*k*i/N);
420 XX[i] = XX[i] + IMX[k] * sin(2*PI*k*i/N);
430 }
440 }
450
460
470
480
490
500 /*Method 2: Loop through each sample in the time domain, and sum the corresponding samples from each cosine and sine wave*/
510 for (int i=0; i<512; i++){
520 for (int k=0; k<257; k++){
530 XX[i] = XX[i] + REX[k] * cos(2*PI*k*i/N);
540 XX[i] = XX[i] + IMX[k] * sin(2*PI*k*i/N);
550 }
560 }
图6展示了逆DFT的操作以及频域信号和合成时域信号所需的波形的幅度差异。
图 6 逆 DFT
图6-a展示了一个待合成的时域信号，是一个冲击信号，其第一个采样值为32，其余采样值均为0。图6-b展示了该信号的频域表示，其实部是值为32的常量，其虚部均为0（这儿没有展示出来）。正如下一章节讨论的那样，这是一个重要的DFT变换对：时域的冲击信号对应于频域的常量。现在要注意的是，图6-b是图6-a的DFT，图6-a是图6-b的逆DFT。式 3 3 3将频域信号图6-b转换为余弦波的幅度，即图6-c。可以看出，除了第一个和最后一个余弦波的幅度为1之外，其余余弦波的幅度均为2。此处没有展示正弦波的幅度，因为，它们均为0值。式 2 2 2将余弦波的幅度图6-b转换为时域信号图6-a。
上述部分描述了频域信号和正弦信号的幅度有何不同，但并没有解释为什么会不同。它们两者之间的差异是因为频域信号是定义在spectral density上的。图7展示了这是如何工作的。它展示了一个32点的时域信号所对应的频域信号的实部。spectral density describes how much signal(amplitude) is preset per unit of bandwidth。为了将正弦信号的幅度转化为spectral density，divide each amplitude by the bandwidth represented by each amplitude。那么，问题来了：如何确定频域每个离散频率所表示的带宽？
持续更新中。。。。。。