一天征服傅里叶变换
如果你对信号处理感兴趣,无疑会说这个标题是太夸张了我赞同这点当然,没有反覆实践和钻研数学,
您无法在一天里学会傅立叶变换的方方面面无论如何,这个在线课程将提供给您怎样进行傅立叶变换运
算的基本知识能有效和能非常简单地领会的原因是我们使用了一种不太传统的逼近重要的是你将学习
傅立叶变换的要素而完全不用超过加法和乘法的数学计算!我将设法在不超过以下六节里解释在对音像信
号处理中傅立叶变换的实际应用
步骤1:一些简单的前提
在下面,您需要理解以下四件最基本的事情:加法,乘除法什么是正弦,余弦和正弦信号明显
地,我将跳第一二件事和将解释位最后一个您大概还记得您在学校学过的“三角函数”[1],它神秘地用于
与角度一起从它们的内角计算它们的边长,反之亦然我们这里不需要所有这些事,我们只需要知道二个
最重要的三角函数,"正弦"和"余弦"的外表特征这相当简单:他们看起来象是以峰顶和谷组成的从观
察点向左右无限伸展的非常简单的波浪
(附图一)
Thesinewave
Thecosinewave
如同你所知道的,这两种波形是周期性的,这意味着在一定的时间周期之后,它们看起来再次一样
两种波形看起来也很象,但当正弦波在零点开始时余弦波开始出现在最大值在实践中,我们如何判定我
们在一个给定时间所观测到的波形是开始在它的最大值或在零?问的好:我们不能实践上没有办法区分
正弦波和余弦波,因此看起来象正弦或余弦波的我们统称为正弦波,在希腊语中译作"正弦类"正弦波的
一个重要性质是"频率"它告诉我们在一个给定的时间内有多少个波峰和波谷高频意味许多波峰和波谷,
低频率意味少量波峰和波谷:
(附图二)
Lowfrequencysinusoid
Middlefrequencysinusoid
Highfrequencysinusoid
步骤2:了解傅立叶定理
Jean-BaptisteJosephFourier是孩子们中让父母感到骄傲和惭愧的的一个,因为他十四岁时就开始对他们
说非常复杂的数学用语他的一生中做了很多重要工作,但最重大的发现可能是解决了材料热传导问题
他推导出了描述热在某一媒介中如何传导的公式,即用三角函数的无穷级数来解决这个问题(就是我们在
上面讨论过的正弦余弦函数)主要和我们话题有关的是:傅里叶的发现总结成一般规律就是任意复杂
的信号都能由一个个混合在一起的正弦函数的和来表示
这是一个例子:
Thisisouroriginal
Onesine
Twosines
Foursines
Sevensines
Fourteensines
(附图三)
在这里你看到的是一个原始的信号,以及如何按某一确定的关系(“配方”)混合在一起的正弦函数混
合物(我们称它们为分量)所逼近我们将简略地谈论一下那份配方如你所知,我们用的正弦函数愈多
其结果就愈精确地接近我们的原始信号波形在“现实”世界中,在信号连续的地方,即你能以无穷小的间
隔来测量它们,精度仅受你的测试设备限制,你需要无限多的正弦函数才能完美地建立任意一个给定的信
号幸运地是,和数字信号处理者们一样,我们不是生活在那样的世界相反,我们将处理仅以有限精度
每隔一定间隔被测量的现实世界的采样信号因而,我们不需要无限多地正弦函数,我们只需要非常多
稍后我们也将讨论这个“非常多”是多少目前重要的一点是你能够想象,任意一个在你计算机上的信号,
都能用简单正弦波按配方组成
步骤3:“非常多”是多少
正如我们所知道的,复杂形状的波形能由混合在一起的正弦波所建立我们也许要问需要多少正弦波来
构造任意一个在计算机上给定的信号当然,倘若我们知道正在处理的信号是如何组成的,这可能至少是
一个单个正弦波在许多情况下,我们处理的现实世界的信号可能有非常复杂的结构,以至于我们不能深
入知道实际上有多少“分量”波存在在这种情况下,即使我们无法知道原始的信号是由多少个正弦波来构
成的,肯定存在一个我们将需要多少正弦波的上限尽管如此,这实际上没解决有多少的问题让我们试
着来直观地逼近它:假设一个信号我们有1000个样采,可能存在的最短周期正弦波(即多数波峰波谷在其
中)以交替的波峰波谷分布在每个采样内因此,最高频率的正弦波将有500个波峰和500个波谷在我们
的1000个采样中,且每隔一个采样是波峰下图中的黑点表示我们的采样,所以,最高频率的正弦波以
看起来象这样:
Thehighestfrequencysinewave
(附图四)
现在让我们来看一下最低频率正弦波可能多么低如果我们只给一个单独的采样点,我们将如何能测
量穿过这点的正弦波的峰顶和谷?我们做不到,因为有许多不同周期正弦波穿过这点
(附图五)
所以,一个单独数据点不足以告诉我们关于频率的任何事现在,如果我们有两个采样,那么穿过这
两点的正弦波的最低频率是什么?在这种情况下它很常简单只有一个穿过这两点的非常低频率的正弦
波它看起来向这样:
Thelowestfrequencysinewave
(附图六)
想象最左面的两个点是二个钉子和一个跨越它们之间的弦(图六描述三个数据点是因为正弦波的周期
性,但我们实际上只需要最左面的两点说明它的频率)我们能领会的最低的频率是来回地摆动在二个钉
子之间的弦,象是图六中左边两点之间的正弦波所做的如果我们有1000个采样,那么两个“钉子”相当
于第一个或最后的采样,比如1号采样和1000号采样从对乐器的体验我们知道,当长度增加时弦的频
率将下降所以我们可以想象,当我们将两个钉子向彼此远离的方向移开时,最小正弦波的频率将变得更
低例如,如果我们选择2000个采样,因为我们的“钉子”,在1号或2000号采样间,所以最低正弦波将
更低事实上,它将低两倍,因为因为我们的钉子比1000个采样时远两倍这样,如果我们有更多的采
样,我们将能辨别出一个更低频率的正弦波,因为它们的零交叉点(我们的“钉子”)将移动得更远这对
了解下面的解释是非常重要的
象我们看到的那样,在两个“钉子”之后,我们的波形开始重复上升斜坡(第一个钉子和第三个钉子同
样)这意味着任意两个相邻的钉子准确地包含完整正弦波的一半,换句话说一个峰或一个谷或半个周期
概括一下我们刚学过的东西,我们知道,一个采样正弦波的上限频率是所有其它有一个波峰和波谷的
采样,并且,低频下限是我们看到的正好匹配采样数的正弦波的周期的一半但等一下,这难到不意味着,
当上限频率保持固定时,当有很多采样时最低频率可以降低?确实如此!结果我们将从一个低频开始增加
更多正弦波来组成一个较长的未知内容的信号
一切都清楚了,但我们仍然不知道我们最终需要多少正弦波就象我们现在知道任意正弦波分量所能
具的有上下限频率一样,我们能计算出适合这两个极限之间的量是多少既然我们沿着最左到最右的采样
固定了最低正弦波分量,我们要所有其它正弦波最好也使用这些钉子(为什么我们要不同地对待他们?所有
的正弦波被同等的创建!)假设正弦波束是系在吉他上二个固定点的弦它们能只摇摆在二个钉子之间(除
非他们断了),就象我们下图的正弦波这导致如下关系,我们的最低分量(1)以1/2周期装配,第二分
量(2)以1个周期装配,第三分量以11/2周期装配以此类推直到我们看到的1000个采样形象地,看
起来象这:
(附图七)
现在,如果我们算一下要多少正弦波以那种方法装配我们的1000个采样,就会发现我们精确地需要
1000个正弦波叠加起来表示1000个采样实际上,我们总是发现我们需要和采样一样多的正弦波
步骤4关于烹饪食谱
在前面的段落我们看到,任一个给定的在计算机上的信号能被正弦波混合物来构造我们考虑了他们
的频率,并且考虑了需要多大的最低和最高频率的正弦波来完美地重建任一个我们所分析信号我们明白
了为确定所需最低的正弦波分量,我们考察的采样的数量是重要,但我们还未论述实际正弦波如何必需被
混合产生某一确定的结果由叠加正弦波组成任何指定的信号,我们需要测量他们的另外一个方面实际
上,频率不是我们需要知道的唯一的事我们还需要知道正弦波的幅度,也就是说每个正弦波幅度有多高
才能混合在一起产生我们需要的输入信号高度是正弦波的峰顶的高度,意即峰顶和零线之间距离幅度
赿高,我们听到的声音也就赿大所以,如果您有一个含有许多低音的信号,无疑可以预期混合体中的低
频率正弦波的分量比例比高频正弦波分量更大因此一般情况下,低音中的低频正弦波有一个比高频正弦
波更高的幅度在我们的分析中,我们将需要确定各个分量正弦波的幅度以完成我们的配方
步骤5:关于苹果和桔子
如果你一直跟着我,我们几乎完成了通向傅里叶变换的旅程我们学了需要多少正弦波,它的数量依
赖于我们查看的采样的数量有一个频率上下限界,并且不知道怎么确定单个分量的幅度以完成我们的配方
我们一直不清楚究竟如何从我们的采样来确定实际的配方直观上我们可以断定能找到正弦波的幅度,
设法把一个已知频率正弦波和采样作对比,我们测量找出它们有多么接近如果它们精确地相等,我们知
道该正弦波存在着相同的幅度,如果我们发现我们的信号与参考正弦波一点也不匹配,我们将认为这个不
存在尽管如此,我们如何高效地把一个已知的正弦波同采样信号进行比较?幸运地是,数字信号处理工
作者早已解决了如何作这些事实上,这象加法和乘法一样容易-我们取一个已知频率的单位正弦波(这
意味着它的振幅是1,可从我们的计算器或计算机中精确地获得)和我们的信号采样相乘累加乘积之后,
我们将得到我们正在观测的这个频率上正弦波分量的幅度这是个举例,一个简单的完成这些工作的C代码
片段:
Listing1.1:ThedirectrealizationoftheDiscreteSineTransform(DST):
#defineM_PI3.14159265358979323846
longbin,k;
doublearg;
for(bin=0;bin transformData[bin]=0.;
for(k=0;k arg=(float)binM_PI(float)k/(float)transformLength;
transformData[bin]+=inputData[k]sin(arg);
}
}
这段代码变换存储在inputData[0...transformLength-1]中我们测量的采样点成为一个正弦波分量的幅度
队列transformData[0...transformLength-1]根据通用术语,我们称参考正弦波的频率步长为盒(bin),这
意味着它们被认为象是一个我们放置我们估计的任意分量波的幅度的容器离散正弦变换(DST)是一个
普通程序,它假设我们无法想象我们的信号看起来象什么样,否则我们能使用一个更加高效率的方法来确
定正弦波分量的幅度(例如,我们予先知道,我们的信号是一个已知频率的正弦波我们能直接地找出它
的高度而不用计算正弦波的整个范围实现这个有效的逼近是基于傅里叶原理,它能在文献的戈策尔(Goe
rtzel)算法条目下找到)
这些就是你坚持想要的我们为什么用这样的方法计算正弦变换的一个解释:对我们用一个已知频率正弦
波的乘积来作一种非常直观逼近的理由,可以设想,这大致相当于一个固有频率的“共振”在系统内发生时
物理世界发生的事情sin(arg)项本质上是一个获得由输入信号波形激励的谐振器如果输入(信号)有在
我们正观测的频率上的分量,它的输出将是参考正弦波谐振的幅度因为我们的参考波是单位幅度的,输
出是一个在那个频率上的分量的实际幅度的一个直接测量因为谐振器只是简单的滤波器,变换(不可否
认是在稍微宽松条件下)被认为有极窄的带通滤波器组的特征,它位于我们估值的频率中心的周围这有
助于解释一个事实,为什么傅立叶变换提供了对信号进行过滤的一个高效工具
只是为了完备性:当然,上述程序是可逆的,当我们知道它的正弦波分量时,我们的信号(在数字精确度
极限内)能完全被重建,通过简单地把正弦波加起来这留下给读者做为一个练习同样程序能改变使用余
弦波做为基本函数工作-我们只需简单地改变sin(arg)条件到cos(arg)来获得离散余弦变换的直接实现(DC
T)
现在,就象在这篇文章较前面的段落中我们讨论过的那样,我们在实践中没有办法区分一个被测量的正
弦类函数象是正弦波还是余弦波做为代替我们总是测量正弦信号,且正弦和余弦变换在实践中没有太大
的用途,除了一些特殊情况(象图象压缩的地方,即毎块图象具有能用一个基本的余弦或正弦函数较好模
拟特性,例如能用余弦基本函数较好表现的相同颜色的大区域)正弦信号是一个比正弦或余弦波更一般
的片断,因为它可以开始在在它的周期中的一个任意位置我们记得,当余弦波开始于1时,正弦波总开
始在0当我们采取正弦波作为参考,余弦波开始在它的周期的最后1/4之处一般用度或弧度测量它们
的偏移量,这是两个一般与三角函数相关的单位一个完整的周期等于360°(代表"度")或2π个弧度(代
表"2π","π"发音象“pie”π是希腊字表示数3.14159265358979323846...在三角学方面有重要意义)
余弦波因而有一个90°或π/2的偏移这偏移叫正弦信号的相位,因此余弦波相对正弦信号有90°或π/2相
位
相位的事情就有这些内容因为我们一直不能限定信号在0°或90°相位开始(因为我们正观测一个我们
可能无法控制的信号),它对同时直接唯一的描述信号的频率振幅相位至关重要以正弦或余弦做变
换相位限制在0°或90°,一个具有任意相位的正弦信号将引起相邻频率出现假峰(因为它们试图“帮助”分
析,强制给被测信号加上一个0°或90°的相位作用)它有些象用一圆石头去填满一个方孔:你需要小一
些的圆石头去填充剩余的空间,并且更小的石头填好依然留出空的空间,等等我们所需要的是能处理一
般信号的变换,它能处理任意相位正弦波构成的信号
步骤6:离散傅叶变换
从正弦变换到傅里叶变换的步骤是简单的,只需用更一般的方法在正弦变换中对每个频率上的测度
使用正弦波,在傅里叶变换中正弦余弦波二者都使用就是说,对任意的当前频率,我们以同一频率的
正弦和余弦波来“比较”(或“共振”)被测信号如果我们的信号看起来很象正弦波,变换的正弦部份将有
一个大的幅值如果它看起来象余弦波,变换的余弦部份将有一个大的幅值如果看起来象反相的正弦波,
也就是说,它开始于0但下降至-1取代上升至1,它的正弦部份将有一个大的负幅值这表明用+-符
号和正弦余弦相位能表示任意给定频率的正弦信号[2]
Listing1.2:ThedirectrealizationoftheDiscreteFourierTransform[3]:
#defineM_PI3.14159265358979323846
longbin,k;
doublearg,sign=-1.;/sign=-1->FFT,1->iFFT/
for(bin=0;bin<=transformLength/2;bin++){
cosPart[bin]=(sinPart[bin]=0.);
for(k=0;k arg=2.(float)binM_PI(float)k/(float)transformLength;
sinPart[bin]+=inputData[k]signsin(arg);
cosPart[bin]+=inputData[k]cos(arg);
}
}
我们仍遗留着一个问题,就是如何获得傅里叶变换所缺乏的那些有用的东西我说过傳里叶变换的优赿性
超过正弦和余弦变换是因为用正弦信号工作但至今我们还未看到任何正弦信号,仍只有正弦和余弦好,
这需要一点附加处理步聚:
#defineM_PI3.14159265358979323846
longbin;
for(bin=0;bin<=transformLength/2;bin++){
/frequency/
frequency[bin]=(float)binsampleRate/(float)transformLength;
/magnitude/
magnitude[bin]=20.log10(2.sqrt(sinPart[bin]sinPart[bin]+cosPart[bin]c
osPart[bin])/(float)transformLength);
/phase/
phase[bin]=180.atan2(sinPart[bin],cosPart[bin])/M_PI-90.;
}
在运行清单1.3所示的关于DFT输出的代码段之后,我们结束被看作以正弦信号波的和的输入信号表示
K序正弦信号是用frequency[k],magnitude[k]和phase[k]来描述的单位是Hz(Hertz,周/秒),dB(Decibel),
和°(Degree)请注意在经过清单1.3的后加工(处理)即把正弦和余弦函数部份转换成一个单一的正弦信
号之后,我们命名K序正弦信号的振幅-DFT存贮为幅度,且它总是取相对值我们可以说一个-1.0的振
幅对应于1.0的幅度,对应于相位+或-180°在文献中,做做傅里叶变换的场合,队列magnitude[]被称作
被测信号的幅度谱,队列phase[]被称作被测信号的相位谱
如用分贝测量存贮幅度的参考,输入波也期望有一个在[-1.0,1.0]之间的采样值,相对于0dB幅度满刻
度数字做为一个DFT的有趣应用,比如清单1.3就可被用于写一个基于离散傅里叶变换的谱分析
结论
象我们已知那样,傅里叶变换和其系列的离散正弦和余弦变换,提供了把一个信号分解成一束分波的
便利工具结果有正弦或余弦之一,或正弦信号(用正弦和余弦波的组合来描述)在傅里叶变换中同时
使用正弦和余弦波的好处是我们因而能引入相位的概念,它使变换更一般化,因而我们能用它有效清楚地
分析既不是纯正弦也不是纯余弦的正弦信号,当然其它信号也一样
傅里叶变换与被考察信号无关,因而无论我们正分析的信号是一个正弦信号或是一些其它的更复杂的,
变换需要相同的操作数这就是为什么傅里叶变换被称做无参数变换的原因,这意味着它对需要的信号“智
能的”分析没有直接的帮助(在考察一个我们已知是一个信号是正弦曲线的情况下,我们更喜欢精确地获
得关于相位,频率,幅度的信息以代替一串在一些预定频率上的正弦和余弦波)
现在我们也知道了我们是在求输入信号在一组固定频率栅格上的值,输入信号实际存在的频率组在这
组栅格上可能不起作用我们在分析中利用的栅格是人为的,因为我们几乎按照关于它们的频率的尝试来
选择参考正弦余弦波说到了这些,马上清楚了一个将要很容易遇到的要点,即被测信号的频率位于变
换栅格的频率之间因此,有一个频率发生在位于两个频率栅格之间的正弦曲线,在变换中将不好被描述
包围着与输入信号频率最接近的栅格的相邻的栅格将试图‘改正’频率的背离因而,输入信号的能量将拖
尾至数个相邻的栅格这也是傅里叶变换不能迅速地分析声音返回它的基波和谐波(并且,这也是为什么
我们称正弦和余弦波为分波而不谐波和泛音)
简单的说,没有进一步的快速处理,DFT和一个狭窄的坝一样,细小并行的带通滤波器组(“通道”)
和每个通道带有附加的相位信息这对分析信号做滤波器和运用其它的技巧是有益的(改变一个信号的
音调而不改变它的速度是它们其中之一,说明在DSPdimension.com上另一篇不同的文章中),但它需要
对少量普通任务附加快速处理同样,它能被认为是使用除了正弦和余弦波基本函数的变换系列的一个特
例在这个方向上展开概念超出了这篇文章的范围
最后,重要的是要提及一个更高效的DFT工具,也就是一个被称做快速傅里叶变换的算法它最初
是由库利和图克在1969年构思的(它的根源仍然要追溯到高斯和其它人的工作)FFT只是一个高效的
算法,它比上面给出的以直接逼近计算DFT所化的时间少,它是结果完全相同的其它方法无论如何,F
FT是以库利/图克算法实施的,它需要变换长度是2的幂在实践中,对大多数应用来说这是一个可以接
受的限制有大量的以不同方法实施FFT的可利用的文献,因而,可以说有足够多不同的FFT实现,其
中一些并不需要经典FFT的2的幂的限制下面清单1.4以程序smbFft()给出了一个FFT的实现
Listing1.4:TheDiscreteFastFourierTransform(FFT):
#defineM_PI3.14159265358979323846
voidsmbFft(floatfftBuffer,longfftFrameSize,longsign)
/
FFTroutine,(C)1996S.M.Bernsee.Sign=-1isFFT,1isiFFT(inverse)
FillsfftBuffer[0...2fftFrameSize-1]withtheFouriertransformofthetimedomaindatainfftBuffer
[0...2fftFrameSize-1].TheFFTarraytakesandreturnsthecosineandsinepartsinaninterleaved
manner,ie.fftBuffer[0]=cosPart[0],fftBuffer[1]=sinPart[0],asf.fftFrameSizemustbeapowerof
2.Itexpectsacomplexinputsignal(seefootnote2),ie.whenworkingwith'common'audiosig
nalsourinputsignalhastobepassedas{in[0],0.,in[1],0.,in[2],0.,...}asf.Inthatcase,thetransform
ofthefrequenciesofinterestisinfftBuffer[0...fftFrameSize].
/
{
floatwr,wi,arg,p1,p2,temp;
floattr,ti,ur,ui,ur,ui,p1r,p1i,p2r,p2i;
longi,bitm,j,le,le2,k;
for(i=2;i<2fftFrameSize-2;i+=2){
for(bitm=2,j=0;bitm<2fftFrameSize;bitm<<=1){
if(i&bitm)j++;
j<<=1;
}
if(i p1=fftBuffer+i;p2=fftBuffer+j;
temp=p1;(p1++)=p2;
(p2++)=temp;temp=p1;
p1=p2;p2=temp;
}
}
for(k=0,le=2;k<(long)(log(fftFrameSize)/log(2.));k++){
le<<=1;
le2=le>>1;
ur=1.0;
ui=0.0;
arg=M_PI/(le2>>1);
wr=cos(arg);
wi=signsin(arg);
for(j=0;j p1r=fftBuffer+j;p1i=p1r+1;
p2r=p1r+le2;p2i=p2r+1;
for(i=j;i<2fftFrameSize;i+=le){
tr=p2rur-p2iui;
ti=p2rui+p2iur;
p2r=p1r-tr;p2i=p1i-ti;
p1r+=tr;p1i+=ti;
p1r+=le;p1i+=le;
p2r+=le;p2i+=le;
}
tr=urwr-uiwi;
ui=urwi+uiwr;
ur=tr;
}
}
}
TheDFT“àPied”:MasteringTheFourierTransforminOneDay
PostedbyBernseeonSeptember21,1999·50Comments
Ifyou’reintosignalprocessing,youwillnodoubtsaythattheheadline
isaverytallclaim.Iwouldsecondthis.Ofcourseyoucan’tlearn
allthebellsandwhistlesoftheFouriertransforminonedaywithout
practisingandrepeatingandeventuallydelvingintothemaths.However,
thisonlinecoursewillprovideyouwiththebasicknowledgeofhowthe
Fouriertransformworks,whyitworksandwhyitcanbeverysimpleto
comprehendwhenwe’reusingasomewhatunconventionalapproach.The
importantpart:youwilllearnthebasicsoftheFouriertransform
completelywithoutanymathsthatgoesbeyondaddingandmultiplying
numbers!IwilltrytoexplaintheFouriertransforminitspractical
applicationtoaudiosignalprocessinginnomorethansixparagraphs
below.
Step1:Somesimpleprerequisites
Whatyouneedtounderstandthefollowingparagraphsareessentiallyfour
things:howtoaddnumbers,howtomultiplyanddividethemandwhata
sine,acosineandasinusoidisandhowtheylook.Obviously,Iwill
skipthefirsttwothingsandjustexplainabitthelastone.You
probablyrememberfromyourdaysatschoolthe‘trigonometric
functions’thatweresomehowmysteriouslyusedinconjunctionwith
trianglestocalculatethelengthofitssidesfromitsinneranglesand
viceversa.Wedon’tneedallthesethingshere,wejustneedtoknowhow
thetwomostimportanttrigonometricfunctions,the?sine?and
?cosine?looklike.Thisisquitesimple:theylooklikeverysimple
waveswithpeaksandvalleysinthemthatstretchouttoinfinityto
theleftandtherightoftheobserver.
Asyoucansee,bothwavesareperiodic,whichmeansthatafteracertain
time,theperiod,theylookthesameagain.Also,bothwaveslookalike,
butthecosinewaveappearstostartatitsmaximum,whilethesinewave
startsatzero.Nowinpractice,howcanwetellwhetherawaveweobserve
atagiventimestartedoutatitsmaximum,oratzero?Goodquestion:
wecan’t.There’snowaytodiscernasinewaveandacosinewavein
practice,thuswecallanywavethatlookslikeasineorcosinewave
a?sinusoid?,whichisGreekandtranslatesto?sinus-like?.An
importantpropertyofsinusoidsis?frequency?,whichtellsushowmany
peaksandvalleyswecancountinagivenperiodoftime.Highfrequency
meansmanypeaksandvalleys,lowfrequencymeansfewpeaksandvalleys:
Step2:UnderstandingtheFourierTheorem
Jean-BaptisteJosephFourierwasoneofthosechildrenparentsareeither
proudorashamedof,ashestartedthrowinghighlycomplicated
mathematicaltermsatthemattheageoffourteen.Althoughhedidalot
ofimportantworkduringhislifetime,theprobablymostsignificant
thinghediscoveredhadtodowiththeconductionofheatinmaterials.
Hecameupwithanequationthatdescribedhowtheheatwouldtravelin
acertainmedium,andsolvedthisequationwithaninfiniteseriesof
trigonometricfunctions(thesinesandcosineswehavediscussedabove).
Basically,andrelatedtoourtopic,whatFourierdiscoveredboilsdown
tothegeneralrulethateverysignal,howevercomplex,canberepresented
byasumofsinusoidfunctionsthatareindividuallymixed.
Whatyouseehereisouroriginalsignal,andhowitcanbeapproximated
byamixtureofsines(wewillcallthempartials)thataremixedtogether
inacertainrelationship(a‘recipe’).Wewilltalkaboutthatrecipe
shortly.Asyoucansee,themoresinesweusethemoreaccuratelydoes
theresultresembleouroriginalwaveform.Inthe‘real’world,where
signalsarecontinuous,ie.youcanmeasurethemininfinitelysmall
intervalsatanaccuracythatisonlylimitedbyyourmeasurement
equipment,youwouldneedinfinitelymanysinestoperfectlybuildany
givensignal.Fortunately,asDSPerswe’renotlivinginsuchaworld.
Rather,wearedealingwithsamplesofsuch‘realworld’signalsthat
aremeasuredatregularintervalsandonlywithfiniteprecision.Thus,
wedon’tneedinfinitelymanysines,wejustneedalot.Wewillalso
talkaboutthat‘howmuchisalot’lateron.Forthemoment,itis
importantthatyoucanimaginethateverysignalyouhaveonyourcomputer
canbeputtogetherfromsimplesinewaves,aftersomecookingrecipe.
Step3:Howmuchis?alot?
Aswehaveseen,complexshapedwaveformscanbebuiltfromamixture
ofsinewaves.Wemightaskhowmanyofthemareneededtobuildanygiven
signalonourcomputer.Well,ofcourse,thismaybeasfewasonesingle
sinewave,providedweknowhowthesignalwearedealingwithismade
up.Inmostcases,wearedealingwithrealworldsignalsthatmighthave
averycomplexstructure,sowedonotknowinadvancehowmany‘partial’
wavesthereareactuallypresent.Inthiscase,itisveryreassuring
toknowthatifwedon’tknowhowmanysinewavesconstitutetheoriginal
signalthereisanupperlimittohowmanywewillneed.Still,thisleaves
uswiththequestionofhowmanythereactuallyare.Let’strytoapproach
thisintuitively:assumewehave1000samplesofasignal.Thesinewave
withtheshortestperiod(ie.themostpeaksandvalleysinit)thatcan
bepresenthasalternatingpeaksandvalleysforeverysample.So,the
sinewavewiththehighestfrequencyhas500peaksand500valleysin
our1000samples,witheveryothersamplebeingapeak.Theblackdots
inthefollowingdiagramdenoteoursamples,sothesinewavewiththe
highestfrequencylookslikethis:
Nowlet’slookhowlowthelowestfrequencysinewavecanbe.Ifweare
givenonlyonesinglesamplepoint,howwouldwebeabletomeasurepeaks
andvalleysofasinewavethatgoesthroughthispoint?Wecan’t,as
therearemanysinewavesofdifferentperiodsthatgothroughthispoint.
So,asingledatapointisnotenoughtotellusanythingaboutfrequency.
Now,ifweweregiventwosamples,whatwouldbethelowestfrequency
sinewavethatgoesthroughthesetwopoints?Inthiscase,itismuch
simpler.Thereisoneverylowfrequencysinewavethatgoesthroughthe
twopoints.Itlookslikethis:
Imaginethetwoleftmostpointsbeingtwonailswithastring
spannedbetweenthem(thediagramdepictsthreedatapointsasthesine
waveisperiodic,butwereallyonlyneedtheleftmosttwototellits
frequency).Thelowestfrequencywecanseeisthestringswinging
backandforthbetweenthetwonails,likeoursinewavedoesinthe
diagrambetweenthetwopointstotheleft.Ifwehave1000samples,
thetwo‘nails’wouldbethefirstandthelastsample,ie.sample
number1andsample
number1000.Weknowfromourexperiencewithmusicalinstrumentsthat
thefrequencyofastringgoesdownwhenitslengthincreases.Sowewould
expectthatourlowestsinewavegetslowerinfrequencywhenwemove
ournailsfartherawayfromeachother.Ifwechoose2000samples,for
instance,thelowestsinewavewillbemuchlowersinceour‘nails’
arenowsamplenumber1andsamplenumber2000.Infact,itwillbetwice
aslow,sinceournailsarenowtwiceasfarawayasinthe1000samples.
Thus,ifwehavemoresamples,wecandiscernsinewavesofalower
frequencysincetheirzerocrossings(our‘nails’)willmovefarther
away.Thisisveryimportanttounderstandforthefollowing
explanations.
Aswecanalsosee,aftertwo‘nails’ourwavestartstorepeatwith
theascendingslope(thefirstandthethirdnailareidentical).This
meansthatanytwoadjacentnailsembraceexactlyonehalfofthe
completesinewave,orinotherwordseitheronepeakoronevalley,or1/2
period.
Summarizingwhatwehavejustlearned,weseethattheupperfrequency
ofasampledsinewaveiseveryothersamplebeingapeakandavalley
andthelowerfrequencyboundishalfaperiodofthesinewavewhich
isjustfittinginthenumberofsampleswearelookingat.Butwait–
wouldn’tthismeanthatwhiletheupperfrequencyremainsfixed,the
lowestfrequencywoulddropwhenwehavemoresamples?Exactly!Theresult
ofthisisthatwewillneedmoresinewaveswhenwewanttoputtogether
longersignalsofunknowncontent,sincewestartoutatalower
frequency.
Allwellandgood,butstillwedon’tknowhowmanyofthesesinewaves
wefinallyneed.Aswenowknowthelowerandupperfrequencyanypartial
sinewavecanhave,wecancalculatehowmanyofthemfitinbetweenthese
twolimits.Sincewehavenailedourlowestpartialsinewavedownto
theleftmostandrightmostsamples,werequirethatallothersinewaves
usethesenailsaswell(whyshouldwetreatthemdifferently?Allsine
wavesarecreatedequal!).Justimaginethesinewaveswerestringson
aguitarattachedtotwofixedpoints.Theycanonlyswingbetweenthe
twonails(unlesstheybreak),justlikeoursinewavesbelow.Thisleads
totherelationshipthatourlowestpartial(1)fitsinwith1/2period,
thesecondpartial(2)fitsinwith1period,thethirdpartial(3)fits
inwith11/2periodasf.intothe1000sampleswearelookingat.
Graphically,thislookslikethis:
Nowifwecounthowmanysinewavesfitinour1000samplesthatway,
wewillfindthatweneedexactly1000sinewavesaddedtogetherto
representthe1000samples.Infact,wewillalwaysfindthatweneed
asmanysinewavesaswehadsamples.
Step4:Aboutcookingrecipes
Inthepreviousparagraphwehaveseenthatanygivensignalonacomputer
canbebuiltfromamixtureofsinewaves.Wehaveconsideredtheir
frequencyandwhatfrequencythelowestandhighestsinewavesneedto
havetoperfectlyreconstructanysignalweanalyze.Wehaveseenthat
thenumberofsampleswearelookingatisimportantfordeterminingthe
lowestpartialsinewavethatisneeded,butwehavenotyetdiscussed
howtheactualsinewaveshavetobemixedtoyieldacertainresult.
Tomakeupanygivensignalbyaddingsinewaves,weneedtomeasureone
additionalaspectofthem.Asamatteroffact,frequencyisnottheonly
thingweneedtoknow.Wealsoneedtoknowtheamplitudeofthesine
waves,ie.howmuchofeachsinewaveweneedtomixtogethertoreproduce
ourinputsignal.Theamplitudeistheheightofthepeaksofasinewave,
ie.thedistancebetweenthepeakandourzeroline.Thehigherthe
amplitude,thelouderitwillsoundwhenwelistentoit.So,ifyouhave
asignalthathaslotsofbassinityouwillnodoubtexpectthatthere
mustbeagreaterportionoflowerfrequencysinewavesinthemixthan
therearehigherfrequencysinewaves.So,generally,thelowfrequency
sinewavesinabassysoundwillhaveahigheramplitudethanthehigh
frequencysinewaves.Inouranalysis,wewillneedtodeterminethe
amplitudeofeachpartialsinewavetocompleteourrecipe.
Step5:Aboutapplesandoranges
Ifyouarestillwithme,wehavealmostcompletedourjourneytowards
theFouriertransform.Wehavelearnedhowmanysinewavesweneed,that
thisnumberdependsonthenumberofsampleswearelookingat,thatthere
isalowerandupperfrequencyboundaryandthatwesomehowneedto
determinetheamplitudeoftheindividualpartialwavestocompleteour
recipe.We’restillnotclear,however,onhowwecandeterminethe
actualrecipefromoursamples.Intuitively,wewouldsaythatwecould
findtheamplitudesofthesinewavessomehowbycomparingasinewave
ofknownfrequencytothesampleswehavemeasuredandfindouthow
‘equal’theyare.Iftheyareexactlyequal,weknowthatthesinewave
mustbepresentatthesameamplitude,ifwefindoursignaltonotmatch
ourreferencesinewaveatallwewouldexpectthisfrequencynottobe
present.Still,howcouldweeffectivelycompareaknownsinewavewith
oursampledsignal?Fortunately,DSPershavealreadyfiguredouthowto
dothisforyou.Infact,thisisaseasyasmultiplingandaddingnumbers
–wetakethe‘reference’sinewaveofknownfrequencyandunit
amplitude(thisjustmeansthatithasanamplitudeof1,whichisexactly
whatwegetbackfromthesin()functiononourpocketcalculatororour
computer)andmultiplyitwithoursignalsamples.Afteraddingtheresult
ofthemultiplicationtogether,wewillobtaintheamplitudeofthe
partialsinewaveatthefrequencywearelookingat.
Toillustratethis,here’sasimpleCcodefragmentthatdoesthis:
Thiscodesegmenttransformsourmeasuredsamplepointsthatarestored
ininputData[0...transformLength-1]intoanarrayofamplitudesofits
partialsinewavestransformData[0...transformLength-1].Accordingto
commonterminology,wecallthefrequencystepsofourreferencesine
wavebins,whichmeansthattheycanbethoughtofasbeing‘containers’
inwhichweputtheamplitudeofanyofthepartialwavesweevaluate.
TheDiscreteSineTransform(DST)isagenericprocedurethatassumes
wehavenoideawhatoursignallookslike,otherwisewecoulduseamore
efficientmethodfordeterminingtheamplitudesofthepartialsinewaves
(ifwe,forexample,knowbeforehandthatoursignalisasinglesine
waveofknownfrequencywecoulddirectlycheckforitsamplitudewithout
calculatingthewholerangeofsinewaves.Anefficientapproachfordoing
thisbasedontheFouriertheorycanbefoundintheliteratureunder
thenamethe?Goertzel?algorithm).
Forthoseofyouwhoinsistonanexplanationforwhywecalculatethe
sinetransformthatway:Asaveryintuitiveapproachtowhywemultiply
withasinewaveofknownfrequency,imaginethatthiscorrespondsroughly
towhatinthephysicalworldhappenswhena‘resonance’atagiven
frequencytakesplaceinasystem.Thesin(arg)termisessentiallya
resonatorthatgetsexcitedbyourinputwaveform.Iftheinputhasa
partialatthefrequencywe’relookingat,itsoutputwillbethe
amplitudeoftheresonancewiththereferencesinewave.Sinceour
referencewaveisofunitamplitude,theoutputisadirectmeasureof
theactualamplitudeofthepartialatthatfrequency.Sincearesonator
isnothingbutasimplefilter,thetransformcan(admittedlyunder
somewhatrelaxedconditions)beseenasahavingthefeaturesofabank
ofverynarrowbandpassfiltersthatarecenteredaroundthefrequencies
we’reevaluating.ThishelpsexplainingthefactwhytheFourier
transformprovidesanefficienttoolforperformingfilteringofsignals.
Justforthesakeofcompleteness:ofcourse,theaboveroutineis
invertible,oursignalcan(withinthelimitsofournumericalprecision)
beperfectlyreconstructedwhenweknowitspartialsinewaves,bysimply
addingsinewavestogether.Thisisleftasanexercisetothereader.
Thesameroutinecanbechangedtoworkwithcosinewavesasbasis
functions–wesimplyneedtochangethesin(arg)termtocos(arg)to
obtainthedirectrealizationoftheDiscreteCosineTransform(DCT).
Now,aswehavediscussedintheveryfirstparagraphofthisarticle,
inpracticewehavenowaytoclassifyameasuredsinus-likefunction
assinewaveorcosinewave.Insteadwearealwaysmeasuringsinusoids,
soboththesineandcosinetransformareofnogreatusewhenweare
applyingtheminpractice,exceptforsomespecialcases(likeimage
compressionwhereeachimagemighthavefeaturesthatarewellmodelled
byacosineorsinebasisfunction,suchaslargeareasofthesamecolor
thatarewellrepresentedbythecosinebasisfunction).Asinusoidis
abitmoregeneralthanthesineorcosinewaveinthatitcanstartat
anarbitrarypositioninitsperiod.Werememberthatthesinewavealways
startsoutatzero,whilethecosinewavestartsoutatone.Whenwetake
thesinewaveasreference,thecosinewavestartsout1/4thlaterin
itsperiod.Itiscommontomeasurethisoffsetindegreeorradians,
whicharetwounitscommonlyusedinconjunctionwithtrigonometric
functions.Onecompleteperiodequals360°(pron.?degree?)or2π
radian(pron.?twopi?with?pi?pronouncedliketheword?pie?.
πisaGreeksymbolforthenumber3.14159265358979323846“whichhas
somesignificanceintrigonometry).Thecosinewavethushasanoffset
of90°orπ/2.Thisoffsetiscalledthephaseofasinusoid,solooking
atourcosinewaveweseethatitisasinusoidwithaphaseoffsetof
90°orπ/2relativetothesinewave.
Sowhat’sthisphasebusinessallabout.Aswecan’trestrictoursignal
tostartoutatzerophaseor90°phaseallthetime(sincewearejust
observingasignalwhichmightbebeyondourcontrol)itisofinterest
todetermineitsfrequency,amplitudeandphasetouniquelydescribeit
atanyonetimeinstant.Withthesineorcosinetransform,we’re
restrictedtozerophaseor90°phaseandanysinusoidthathasan
arbitraryphasewillcauseadjacentfrequenciestoshowspuriouspeaks
(sincetheytryto‘help’theanalysistoforce-fitthemeasuredsignal
toasumofzeroor90°phasefunctions).It’sabitliketryingto
fitaroundstoneintoasquarehole:youneedsmallerroundstonesto
fillouttheremainingspace,andevenmoreevensmallerstonestofill
outthespacethatisstillleftempty,andsoon.Sowhatweneedis
atransformthatisgeneralinthatitcandealwithsignalsthatare
builtofsinusoidsofarbitraryphase.
Step6:TheDiscreteFouriertransform.
ThestepfromthesinetransformtotheFouriertransformissimple,
makingitinawaymore‘general’.Whilewehavebeenusingasinewave
foreachfrequencywemeasureinthesinetransform,weusebothasine
andacosinewaveintheFouriertransform.Thatis,foranyfrequency
wearelookingatwe‘compare’(or‘resonate’)ourmeasuredsignal
withbothacosineandasinewaveofthesamefrequency.Ifoursignal
looksmuchlikeasinewave,thesineportionofourtransformwillhave
alargeamplitude.Ifitlookslikeacosinewave,thecosinepartof
ourtransformwillhavealargeamplitude.Ifitlooksliketheopposite
ofasinewave,thatis,itstartsoutatzerobutdropsto-1instead
ofgoingupto1,itssineportionwillhavealargenegativeamplitude.
Itcanbeshownthatthe+and–signtogetherwiththesineandcosine
phasecanrepresentanysinusoidatthegivenfrequency.
We’restillleftwiththeproblemofhowtogetsomethingusefulout
oftheFourierTransform.IhaveclaimedthatthebenefitoftheFourier
transformovertheSineandCosinetransformisthatweareworkingwith
sinusoids.However,wedon’tseeanysinusoidsyet,therearestillonly
sinesandcosines.Well,thisrequiresanadditionalprocessingstep:
AfterrunningthecodefragmentshowninListing1.3onourDFToutput,
weendupwitharepresentationoftheinputsignalasasumofsinusoid
waves.Thek-thsinusoidisdescribedbyfrequency[k],magnitude[k]
andphase[k].UnitsareHz(Hertz,periodsperseconds),dB(Decibel)
and
°(Degree).Pleasenotethatafterthepost-processingofListing1.3
thatconvertsthesineandcosinepartsintoasinglesinusoid,wename
theamplitudeofthek-thsinusoidtheDFTbin?magnitude?,asitwill
nowalwaysbeapositivevalue.Wecouldsaythatanamplitudeof-
1.0correspondstoamagnitudeof1.0andaphaseofeither+or-180°.
Intheliterature,thearraymagnitude[]iscalledtheMagnitude
Spectrumofthemeasuredsignal,thearrayphase[]iscalledthePhase
SpectrumofthemeasuredsignalatthetimewherewetaketheFourier
transform.
Asareferenceformeasuringthebinmagnitudeindecibels,ourinput
waveisexpectedtohavesamplevaluesintherange[-1.0,1.0),which
correspondstoamagnitudeof0dBdigitalfullscale(DFS).Asan
interestingapplicationoftheDFT,listing1.3can,forexample,beused
towriteaspectrumanalyzerbasedontheDiscreteFourierTransform.
Conclusion
Aswehaveseen,theFouriertransformandits‘relatives’,thediscrete
sineandcosinetransformprovidehandytoolstodecomposeasignalinto
abunchofpartialwaves.Theseareeithersineorcosinewaves,or
sinusoids(describedbyacombinationofsineandcosinewaves).The
advantageofusingboththesineandcosinewavesimultaneouslyinthe
Fouriertransformisthatwearethusabletointroducetheconceptof
phasewhichmakesthetransformmoregeneralinthatwecanuseitto
efficientlyandclearlyanalyzesinusoidsthatareneitherapuresine
orcosinewave,andofcourseothersignalsaswell.
TheFouriertransformisindependentofthesignalunderexaminationin
thatitrequiresthesamenumberofoperationsnomatterifthesignal
weareanalyzingisonesinglesinusoidorsomethingelsemorecomplicated.
ThisisthereasonwhytheDiscreteFouriertransformiscalleda
nonparametrictransform,meaningthatitisnotdirectlyhelpfulwhen
an‘intelligent’analysisofasignalisneeded(inthecasewherewe
areexaminingasignalthatweknowisasinusoid,wewouldpreferjust
gettinginformationaboutitsphase,frequencyandmagnitudeinsteadof
abunchofsineandcosinewavesatsomepredefinedfrequencies).
Wenowalsoknowthatweareevaluatingourinputsignalatafixed
frequencygrid(ourbins)whichmayhavenothingtodowiththeactual
frequenciespresentinourinputsignal.Sincewechooseourreference
sineandcosinewaves(almost)accordingtotastewithregardtotheir
frequency,thegridweimposeonouranalysisisartificial.Havingsaid
this,itisimmediatelyclearthatonewilleasilyencounterascenario
wherethemeasuredsignal’sfrequenciesmaycometoliebetweenthe
frequenciesofourtransformbins.Consequently,asinusoidthathasa
frequencythathappenstoliebetweentwofrequency‘bins’willnot
bewellrepresentedinourtransform.Adjacentbinsthatsurroundthe
binclosestinfrequencytoourinputwavewilltryto‘correct’the
deviationinfrequencyandthustheenergyoftheinputwavewillbe
smearedoverseveralneighbouringbins.Thisisalsothemainreasonwhy
theFouriertransformwillnotreadilyanalyzeasoundtoreturnwith
itsfundamentalandharmonics(andthisisalsowhywecallthesineand
cosinewavespartials,andnotharmonics,orovertones).
Simplyspeaking,withoutfurtherpost-processing,theDFTislittlemore
thanabankofnarrow,slightlyoverlappingbandpassfilters
(’channels’)withadditionalphaseinformationforeachchannel.It
isusefulforanalyzingsignals,doingfilteringandapplyingsomeother
neattricks(changingthepitchofasignalwithoutchangingitsspeed
isoneofthemexplainedinadifferentarticleonDSPdimension.com),
butitrequiresadditionalpostprocessingforlessgenerictasks.Also,
itcanbeseenasaspecialcaseofafamilyoftransformsthatusebasis
functionsotherthanthesineandcosinewaves.Expandingtheconcept
inthisdirectionisbeyondthescopeofthisarticle.
Finally,itisimportanttomentionthatthereisamoreefficient
implementationoftheDFT,namelyanalgorithmcalledthe?FastFourier
Transform?(FFT)whichwasoriginallyconceivedbyCooleyandTukeyin
1969(itsrootshowevergobacktotheworkofGaussandothers).The
FFTisjustanefficientalgorithmthatcalculatestheDFTinlesstime
thanourstraightforwardapproachgivenabove,itisotherwiseidentical
withregardtoitsresults.However,duetothewaytheFFTisimplemented
intheCooley/Tukeyalgorithmitrequiresthatthetransformlengthbe
apowerof2.Inpractice,thisisanacceptableconstraintformost
applications.TheavailableliteratureondifferentFFTimplementations
isvast,sosufficeittosaythattherearemanydifferentFFT
implementations,someofwhichdonothavethepower-of-tworestriction
oftheclassicalFFT.AnimplementationoftheFFTisgivenbytheroutine
smbFft()inListing1.4below.
)simplyspeaking,trigonometricfunctionsarefunctionsthatareused
tocalculatetheanglesinatriangle(?tri-gonos?=Greekfor?three
corners?)fromthelengthofitssides,namelysinus,cosinus,tangent
andthearcustangent.Thesinusandcosinusfunctionsarethemost
importantones,asthetangentandarcustangentcanbeobtainedfrom
sinusandcosinusrelationshipsalone
)Notethatintheliterature,duetoageneralizationthatismade
fortheFouriertransformtoworkwithanothertypeofinputsignalcalled
a‘complexsignal’(complexinthiscontextreferstoacertaintype
ofnumbersratherthantoaninputsignalthathasacomplex
harmonicstructure),youwillencounterthesineandcosinepartunder
thename
‘real’(forthecosinepart)and‘imaginary’part(forthesinepart).
)Ifyou’realreadyacquaintedwiththeDFTyoumayhavenotedthat
thisisactuallyanimplementationofthe?realDiscreteFourier
Transform?,asitusesonlyrealnumbersasinputanddoesnotdealwith
negativefrequencies:intherealDFTpositiveandnegativefrequencies
aresymmetricandthusredundant.Thisiswhywe’recalculatingonly
almosthalfasmanybinsthaninthesinetransform(wecalculateone
additionalbinforthehighestfrequency,forsymmetryreasons).
|
|