Admin: Created page with "理解离散傅立叶变换（二）
&n..."

2011-05-30T15:34:26Z

Created page with "<span style="font-size: 15pt;">理解离散傅立叶变换（二）</span> <div><span> &n..."

New page

<span style="font-size: 15pt;">理解离散傅立叶变换（二）</span>
<div><span>                                ------</span>实数形式离散傅立叶变换（Real DFT）</div>

<div><span>       </span>上一节我们看到了一个实数形式离散傅立叶变换的例子，通过这个例子能够让我们先对傅立叶变换有一个较为形象的感性认识，现在就让我们来看看实数形式离散傅立叶变换的正向和逆向是怎么进行变换的。在此，我们先来看一下频率的多种表示方法。</div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: -21pt;"><strong><span style="font-size: 12pt;">一、<span style="font: 7pt 'Times New Roman';">          </span></span></strong><strong><span style="font-size: 12pt;">频域中关于频率的四种表示方法</span></strong></div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 39pt; text-indent: -18pt;"><span>1、</span>序号表示方法，根据时域中信号的样本数取0 ~ N/2，用这种方法在程序中使用起来可以更直接地取得每种频率的幅度值，因为频率值跟数组的序号是一一对应的: X[k]，取值范围是0 ~ N/2；</div>
<div style="margin: 0cm 0cm 0pt 39pt; text-indent: -18pt;"><span>2、</span>分数表示方法，根据时域中信号的样本数的比例值取0 ~ 0.5: X[ƒ]，ƒ = k/N，取值范围是0 ~ 1/2；</div>
<div style="margin: 0cm 0cm 0pt 39pt; text-indent: -18pt;"><span>3、</span>用弧度值来表示，把<span>ƒ乘以一个2π得到一个弧度值，这种表示方法叫做自然频率(natural frequency)：X[ω]，ω = 2πƒ = 2πk/N，</span>取值范围是0 ~ π；</div>

<div style="margin: 0cm 0cm 0pt 39pt; text-indent: -18pt;"><span>4、</span>以赫兹(Hz)为单位来表示，这个一般是应用于一些特殊应用，如取样率为10 kHz表示每秒有10,000个样本数：取值范围是<span>0到取样率的一半。</span></div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: -21pt;"><strong><span style="font-size: 12pt;">二、<span style="font: 7pt 'Times New Roman';">          </span></span></strong><strong><span style="font-size: 12pt;">DFT基本函数</span></strong></div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 84pt;"><span style="font-size: 14pt;">c<sub>k</sub></span>[i] = cos(2πki/N)</div>
<div style="margin: 0cm 0cm 0pt 84pt;"><span style="font-size: 14pt;">s<sub>k</sub></span>[i] = sin(2πki/N)</div>

<div><span>    其中k表示每个正余弦波的频率，如为2表示在0到N长度中存在两个完整的周期，10即有10个周期，如下图：</span></div>
<div></div>
<div><span>[[File:fft21.jpg]]</span></div>
<div><span>       </span>上图中至于每个波的振幅(amplitude)值(Re X[k],Im X[k])是怎么算出来的,这个是DFT的核心，也是最难理解的部分，我们先来看看如何把分解出来的正余弦波合成原始信号(Inverse DFT)。</div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: -21pt;"><strong><span style="font-size: 14pt;">三、<span style="font: 7pt 'Times New Roman';">        </span></span></strong><strong><span style="font-size: 14pt;">合成运算方法</span></strong><strong><span style="font-size: 14pt;">(Real Inverse DFT)</span></strong></div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 21pt;">DFT合成等式：</div>

<div style="margin: 0cm 0cm 0pt 21pt;">[[File:fft22.jpg]]</div>
<div style="margin: 0cm 0cm 0pt 21pt;">如果有学过傅立叶级数，对这个等式就会有似曾相识的感觉，不错！这个等式跟傅立叶级数是非常相似的：</div>
<div style="margin: 0cm 0cm 0pt 21pt;">[[File:fft23.jpg]]</div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>           </span>当然，差别是肯定是存在的，因为这两个等式是在两个不同条件下运用的，至于怎么证明DFT合成公式，这个我想需要非常强的高等数学理论知识了，这是研究数学的人的工作，对于普通应用者就不需要如此的追根究底了，但是傅立叶级数是好理解的，我们起码可以从傅立叶级数公式中看出DFT合成公式的合理性。</div>
<div style="margin: 0cm 0cm 0pt 21pt;"> </div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       DFT</span>合成等式中的Im <span style="position: relative; top: 2pt;">[k]</span>和Re <span style="position: relative; top: 2pt;">[k]</span>跟Im X[k]和Re X[k]是不一样的，下面是转换方法:</div>

<div style="margin: 0cm 0cm 0pt 21pt;">[[File:fft24.jpg]]</div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span></div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span>但k等于0和N/2时,实数部分的计算要用下面的等式:</div>
<div style="margin: 0cm 0cm 0pt 21pt;">[[File:fft25.jpg]]</div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>              </span></div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span>上面四个式中的N是时域中点的总数，k是从0到N/2的序号。</div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span>为什么要这样进行转换呢？这个可以从频谱密度(spectral density)得到理解，如下图就是个频谱图：</div>
<div style="margin: 0cm 0cm 0pt 21pt;">[[File:fft26.jpg]]</div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span></div>

<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span>这是一个频谱图，横坐标表示频率大小，纵坐标表示振幅大小，原始信号长度为N（这里是32），经DFT转换后得到的17个频率的频谱，频谱密度表示每单位带宽中为多大的振幅，那么带宽是怎么计算出来的呢？看上图，除了头尾两个，其余点的所占的宽度是2/N，这个宽度便是每个点的带宽，头尾两个点的带宽是1/N,而Im X[k]和Re X[k]表示的是频谱密度，即每一个单位带宽的振幅大小，但Im <span style="position: relative; top: 2pt;">[k]</span>和Re <span style="position: relative; top: 2pt;">[k]</span>表示2/N（或1/N）带宽的振幅大小，所以Im <span style="position: relative; top: 2pt;">[k]</span>和Re <span style="position: relative; top: 2pt;">[k]</span>分别应当是Im X[k]和Re X[k]的2/N（或1/N）。</div>
<div style="margin: 0cm 0cm 0pt 21pt;"> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: 21pt;">频谱密度就象物理中物质密度，原始信号中的每一个点就象是一个混合物，这个混合物是由不同密度的物质组成的，混合物中含有的每种物质的质量是一样的，除了最大和最小两个密度的物质外，这样我们只要把每种物质的密度加起来就可以得到该混合物的密度了，又该混合物的质量是单位质量，所以得到的密度值跟该混合物的质量值是一样的。</div>
<div style="margin: 0cm 0cm 0pt 21pt;"> </div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span>至于为什么虚数部分是负数，这是为了跟复数DFT保持一致，这个我们将在后面会知道这是数学计算上的需要（Im X[k]计算时加上了一个负号，Im <span style="position: relative; top: 2pt;">[k]</span>再加上负号，结果便是正的，等于没有变化）。</div>

<div style="margin: 0cm 0cm 0pt 21pt;"> </div>
<div style="margin: 0cm 0cm 0pt 21pt;"><span>       </span>如果已经得到了DFT结果，这时要进行逆转换，即合成原始信号，则可按如下步骤进行转换：</div>
<div style="margin: 0cm 0cm 0pt 60pt; text-indent: -18pt;"><span>1、</span>先根据上面四个式子计算得出Im <span style="position: relative; top: 2pt;">[k]</span>和Re <span style="position: relative; top: 2pt;">[k]</span>的值；</div>
<div style="margin: 0cm 0cm 0pt 60pt; text-indent: -18pt;"><span>2、</span>再根据DFT合成等式得到原始信号数据。</div>
<div style="margin: 0cm 0cm 0pt 42pt;">下面是用BASIC语言来实现的转换源代码：</div>
<div style="margin: 0cm 0cm 0pt 42pt;">100 ‘DFT逆转换方法</div>

<div style="margin: 0cm 0cm 0pt 42pt;">110 ‘/XX[]数组存储计算结果（时域中的原始信号）</div>
<div style="margin: 0cm 0cm 0pt 42pt;">120 ‘/REX[]数组存储频域中的实数分量，IMX[]为虚分量</div>
<div style="margin: 0cm 0cm 0pt 42pt;">130 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">140 DIM XX[511]</div>
<div style="margin: 0cm 0cm 0pt 42pt;">150 DIM REX[256]</div>
<div style="margin: 0cm 0cm 0pt 42pt;">160 DIM IMX[256]</div>
<div style="margin: 0cm 0cm 0pt 42pt;">170 ‘</div>

<div style="margin: 0cm 0cm 0pt 42pt;">180 PI = 3.14159265</div>
<div style="margin: 0cm 0cm 0pt 42pt;">190 N% = 512</div>
<div style="margin: 0cm 0cm 0pt 42pt;">200 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">210 GOSUB XXXX ‘转到子函数去获取REX[]和IMX[]数据</div>
<div style="margin: 0cm 0cm 0pt 42pt;">220 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">230 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">240 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">250 FOR K% = 0 TO 256</div>

<div style="margin: 0cm 0cm 0pt 42pt;">260<span>   REX[K%] = REX[K%] / (N%/2)</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">270<span>   IMX[K%] = -IMX[K%] / (N%/2)</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">280 NEXT k%</div>
<div style="margin: 0cm 0cm 0pt 42pt;">290 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">300 REX[0] = REX[0] / N</div>
<div style="margin: 0cm 0cm 0pt 42pt;">310 REX[256] = REX[256] / N</div>
<div style="margin: 0cm 0cm 0pt 42pt;">320 ‘</div>

<div style="margin: 0cm 0cm 0pt 42pt;">330 ‘ 初始化XX[]数组</div>
<div style="margin: 0cm 0cm 0pt 42pt;">340 FOR I% = 0 TO 511</div>
<div style="margin: 0cm 0cm 0pt 42pt;">350<span>   XX[I%] = 0</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">360 NEXT I%</div>
<div style="margin: 0cm 0cm 0pt 42pt;">370 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">380 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">390 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">400 ‘</div>

<div style="margin: 0cm 0cm 0pt 42pt;">410 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">420 FOR K% =0 TO 256</div>
<div style="margin: 0cm 0cm 0pt 42pt;">430<span>   FOR I%=0 TO 511</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">440 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">450<span>      XX[I%] = XX[I%] + REX[K%] * COS(2 * PI * K% * I% / N%) </span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">460<span>      XX[I%] = XX[I%] + IMX[K%] * SIN(2 * PI * K% * I% / N%)</span></div>

<div style="margin: 0cm 0cm 0pt 42pt;">470 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">480<span>   NEXT I%</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">490 NEXT K%</div>
<div style="margin: 0cm 0cm 0pt 42pt;">500 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">510 END</div>
<div style="margin: 0cm 0cm 0pt 42pt;"> </div>
<div style="margin: 0cm 0cm 0pt 42pt;">上面代码中420至490换成如下形式也许更好理解，但结果都是一样的：</div>
<div style="margin: 0cm 0cm 0pt 42pt;">420 FOR I% =0 TO 511</div>

<div style="margin: 0cm 0cm 0pt 42pt;">430<span>   FOR K%=0 TO 256</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">440 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">450<span>      XX[I%] = XX[I%] + REX[K%] * COS(2 * PI * K% * I% / N%) </span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">460<span>      XX[I%] = XX[I%] + IMX[K%] * SIN(2 * PI * K% * I% / N%)</span></div>
<div style="margin: 0cm 0cm 0pt 42pt;">470 ‘</div>
<div style="margin: 0cm 0cm 0pt 42pt;">480<span>   NEXT I%</span></div>

<div style="margin: 0cm 0cm 0pt 42pt;">490 NEXT K%</div>
<div> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: -21pt;"><strong><span style="font-size: 14pt;">四、<span style="font: 7pt 'Times New Roman';">        </span></span></strong><strong><span style="font-size: 14pt;">分解运算方法（</span></strong><strong><span style="font-size: 14pt;">DFT</span></strong><strong><span style="font-size: 14pt;">）</span></strong></div>
<div style="margin: 0cm 0cm 0pt 21pt;"> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: 21pt;">有三种完全不同的方法进行DFT：一种方法是通过联立方程进行求解, 从代数的角度看，要从N个已知值求N个未知值，需要N个联立方程，且N个联立方程必须是线性独立的，但这是这种方法计算量非常的大且极其复杂，所以很少被采用；第二种方法是利用信号的相关性（correlation）进行计算，这个是我们后面将要介绍的方法；第三种方法是快速傅立叶变换（FFT），这是一个非常具有创造性和革命性的的方法，因为它大大提高了运算速度，使得傅立叶变换能够在计算机中被广泛应用，但这种算法是根据复数形式的傅立叶变换来实现的，它把N个点的信号分解成长度为N的频域，这个跟我们现在所进行的实域DFT变换不一样，而且这种方法也较难理解，这里我们先不去理解，等先理解了复数DFT后，再来看一下FFT。有一点很重要，那就是这三种方法所得的变换结果是一样的，经过实践证明，当频域长度为32时，利用相关性方法进行计算效率最好，否则FFT算法效率较高。现在就让我们来看一下相关性算法。</div>
<div style="margin: 0cm 0cm 0pt 42pt;"> </div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: 21pt;">利用信号的相关性(correlation)可以从噪声背景中检测出已知的信号，我们也可以利用这个方法检测信号波中是否含有某个频率的信号波：把一个待检测信号波乘以另一个信号波，得到一个新的信号波，再把这个新的信号波所有的点进行相加，从相加的结果就可以判断出这两个信号的相似程度。如下图：</div>
<div style="margin: 0cm 0cm 0pt 21pt; text-indent: 21pt;">[[File:fft27.jpg]]</div>
<div><span>        </span>上面a和 b两个图是待检测信号波，图a很明显可以看出是个3个周期的正弦信号波，图b的信号波则看不出是否含有正弦或余弦信号，图c和d都是个3个周期的正弦信号波，图e和f分别是a、b两图跟c、d两图相乘后的结果，图e所有点的平均值是0.5，说明信号a含有振幅为1的正弦信号c，但图f所有点的平均值是0，则说明信号b不含有信号d。这个就是通过信号相关性来检测是否含有某个信号的方法。</div>

<div> </div>
<div><span>       </span>相应地，我也可以通过把输入信号和每一种频率的正余弦信号进行相乘（关联操作），从而得到原始信号与每种频率的关联程度（即总和大小），这个结果便是我们所要的傅立叶变换结果，下面两个等式便是我们所要的计算方法：</div>
<div>[[File:fft28.jpg]]</div>
<div><span>       </span></div>
<div><span>       </span>第二个式子中加了个负号，是为了保持复数形式的一致，前面我们知道在计算Im <span style="position: relative; top: 2pt;">[k]</span>时又加了个负号，所以这只是个形式的问题，并没有实际意义，你也可以把负号去掉，并在计算Im <span style="position: relative; top: 2pt;">[k]</span>时也不加负号。</div>
<div> </div>
<div><span>       </span>这里有一点必须明白一个正交的概念：两个函数相乘，如果结果中的每个点的总和为0，则可认为这两个函数为正交函数。要确保关联性算法是正确的，则必须使得跟原始信号相乘的信号的函数形式是正交的，我们知道所有的正弦或余弦函数是正交的，这一点我们可以通过简单的高数知识就可以证明它，所以我们可以通过关联的方法把原始信号分离出正余弦信号。当然，其它的正交函数也是存在的，如：方波、三角波等形式的脉冲信号，所以原始信号也可被分解成这些信号，但这只是说可以这样做，却是没有用的。</div>

<div><span>       </span>下面是实域傅立叶变换的BASIC语言代码：</div>
<div> [[File:fft29.jpg]]</div>
<div><span>       </span>到此为止，我们对傅立叶变换便有了感性的认识了吧。但要记住，这只是在实域上的离散傅立叶变换，其中虽然也用到了复数的形式，但那只是个替代的形式，并无实际意义，现实中一般使用的是复数形式的离散傅立叶变换，且快速傅立叶变换是根据复数离散傅立叶变换来设计算法的，在后面我们先来复习一下有关复数的内容，然后再在理解实域离散傅立叶变换的基础上来理解复数形式的离散傅立叶变换。</div>

理解离散傅立叶变换（二） - Revision history

Admin: Created page with "理解离散傅立叶变换（二） &n..."

Admin: Created page with "理解离散傅立叶变换（二）
&n..."