手机 网站 开发,淘宝店铺网站策划书,wordpress如何清缓存,wordpress 内容引言 目标 傅里叶变化#xff08;Fourier transform#xff09;是一种信号处理技术#xff0c;它可以将时间信号转换为频率信号#xff0c;即将一组具有相同数量频率的正弦波叠加在一起#xff0c;形成一组新的正弦波。如果我们把时间信号从频域转换到时域#xff0c;那么…引言 目标 傅里叶变化Fourier transform是一种信号处理技术它可以将时间信号转换为频率信号即将一组具有相同数量频率的正弦波叠加在一起形成一组新的正弦波。如果我们把时间信号从频域转换到时域那么就会得到一个“频谱”。所以傅里叶变化就是通过分析不同频率正弦波的相位差来确定它们的频率。 傅里叶变化在信号处理中有广泛应用。比如可以用它来检测和定位信号中的突变点可以用它来确定信号中是否存在非周期信号可以用它来分析复杂系统的行为还可以用它来计算复杂系统的频谱等。 那么我们为什么要使用傅里叶变化呢 意义 如果你是一个对音乐敏感的人那么你一定会发现我们今天讨论的这个主题似乎非常熟悉。在任何技术领域几乎都能找到类似的事情从音频和视频的分析中发现信号中的非周期性成分。当然这种分析是复杂的有时甚至是不可能的因为在现实世界中这些信号可能是由许多不同的原因产生的。但是它们通常都会以特定的方式出现。这就是为什么我们要使用傅里叶变化。 傅里叶变化在技术分析中有很多用途。它可以帮助我们确定信号是否包含一些频率之外的频率成分。或者如果我们想找出特定信号或噪声的频率范围我们可以使用傅里叶变化来确定它是否会出现在我们要找的范围内。 如果你还不知道如何使用傅里叶变化那么这篇博文可以帮您非常深入的了解。
背景知识 声音 我们知道声音是一种机械波是由振动产生的。我们生活的世界就是由各种不同频率、不同振幅、不同方向的振动组成的。通常把频率范围为20Hz-20kHz的振动称为声音频率范围为1 kHz-20 kHz的振动称为乐音频率范围为2 kHz-20 kHz的振动称为声音频率范围为3 kHz-20 kHz的振动称为高音。另外还存在有某些频率范围在10-6~10-7 Hz之间如人耳可听到的次声波的声音。可以说每个人都能发出声音但发出什么样的声音取决于振动器官所处环境和声波频率。但是声音究竟是什么样子的呢我们使用一种特殊的拍照可以查看声波穿过空气分子时看到空气分子的运动如图1.1。这是Michael Hargather 博士和 Gary Settles 博士进行的一项实验他们在纹影装置内的摄像头前放置了一个发出单一音调的扬声器。并拍摄了周围空气分子发生的情况。中间红色的椭圆是扬声器通过空气远离它的是声波当空气分子聚在一起时会得到一个波峰显示出一个高压区 当空气分子再次被拉开时会看到一个低气压区域的形成。
图1.1 当在扬声器面前放置一个假球时可假想球经过声波的震动如何上下摆动的如图1.2。一直持续下去我们可以绘制出声音波形的表示。而我们在数学上有一个函数就可以表示此波形SIN函数。 图1.2 如果在三角形上取①点的正弦值并将其乘以斜边的长度就可以计算出三角形定点的高度如图1.3。 图1.3 随着我们三角形的角度增加它会绕着旋转三角形的顶点成为一个圆圈如果把它对应在二维坐标系中三角形旋转的不同角度的形状看起来与之前的声波信号非常相似图1.4。这就是为什么sin信号是声音一个很好的近似值。但如果只是单纯的一个正弦波那么听起来就是一种很刺耳的声音如果想让正弦波听起来想现实生活中听到的声音呢 图1.4 稍微以不同的方式看待声波信号的波浪假设正弦波的基本形状永久保持不变它总是随着时间的推移上升和下降的。那么我们可以改变正弦波的那个属性频率、幅度、相位先从频率改变缩短正弦波或者拉长正弦波幅度可以拉高正弦波或拉低正弦波。最终所呈现的声音也会随之发生改变。如图1.5表示的是在时域下查看幅度与频率的关系。如图1.6所表示的是在频域下查看幅度与频率的关系。 图1.5 图1.6 而多频率信号叠加在一起时候的时域下如图1.7。频域下查看如图1.8. 图1.7 图1.8 显示世界中发出的声音实际是不同频率和振幅的正弦波这一系列的正弦波加在一起就是所谓的傅里叶系列。 相位 取一个频率为440Hz的正弦波与一个频率为660Hz的正弦波将二者叠加在一起然后开始改变660Hz的相位相加后的波形会随着第二个波形的相位改变而发生变化。因此如果想要使用傅里叶级数准确的表示该信号为一系列正弦波那么每个波的相位信息也很重要因为它会印象信号的整体形状那么如果在数学上表示这种相移呢如图1.9。傅里叶曾经说过一个变量的任何函数都可以在该变量的倍数的一系列正弦中展开。当我们谈论声音信号时我们的变量在X轴上发生是时间而在当时傅里叶提出他的理论时他正在测量热流所以他的变量很可能是距离因此需要一些可以应用于所有类型信号的通用变量。在谈论正弦波时角度是常见的变量然后可以缩短此变量以用于我们正在处理的任何信号。如果我们有一个频率为660Hz的声波这意味着我们耳朵周围的空气正在以每秒660Hz的速度被压缩和膨胀。意味着1个压缩/扩展周期所需的时间是1/660秒。 图1.9 完成一个正弦波的一个周期需要360°也就是在园内旋转一周如图2.0。因此我们可以将声音信号的一个周期所花费的时间映射到360°是圆上因此如果我们想将信号相移30°则30°是360°的12度正弦波的一个周期因此对于需要1/660秒完成的一个周期的660Hz信号1/660秒的12秒0.000126秒如图2.1. 图2.0 图2.1 因此要将正弦波的角度变量转换为信号的时间变量我们将角度除以360°并乘以周期信号在本例中为1/660秒。需记住乘以信号周期等于除以频率。频率是每秒的信号周期数660Hz表示每秒660个完整的正弦波周期。如果我们将信号的时间变量转换为正弦波的角度变量而在这种情况下我们除以信号周期并乘以360°或者因为这次我们除以周期我们可以简单的乘以频率。因此要将时间转换为角度我们乘以频率然后乘以360°。 所以我们在数学上表示30°的相移注意看第二个正弦波比第一个正弦波滞后30°如图2.2。如果角度为30°蓝色正弦波的幅度为0.5红色正弦波不会达到0.5。当角度为60°再达到30°。因此红色波的方程SINX-30°。 图2.2 如果我们将红波移向另一个方向使红波领先与蓝波30°。而红波现在的方程为SINX30°如图2.3。 图2.3 如果我们将正弦波移动90°会得到另一个三角函数余弦函数描述的波。所以SINx90° COSX。 如图2.4所示.如果我们取这个角的正弦值乘以斜边的长度可以计算出三角形的高度。取角度的余弦值并乘以斜边长度可以计算出三角的宽度。 图2.4 如果我们将三角形旋转90°三角的宽度现在变成了它的高度。所以任何角的正弦称为sin(X90°) 等于 该角的余弦值。所以到目前为止有两种方法可以查看正弦和余弦函数的办法。要么相位可以改变的波要么是角度可以改变的三角形。 当我们查看旋转三角形时看到这个角度对于改变正弦波相位的量。我们可以通过两个不同的波余弦波和正弦波相加来描述任何相移相同的频率只需某种方法改变它们的振幅。如图2.5 图2.5 我们将从振幅为32的余弦开始如图2.6可以看到这只是一个没有相移的普通的余弦波。还有一个振幅位29的正弦波同样是一个没有相移的普通正弦波。 图2.6 将两个波形进行相加如图2.7。这是一个大约偏移42°的余弦波如果想知道它的振幅是多少首先回到三角形建模。如图2.8。余弦波为32幅度由三角形底边长表示正弦波为29幅度由三角形高长度表示。那么斜边的长度需要用到勾股定理求出斜边上的平方等于其他两条边上的平方和。而这就是相移后的波形幅度值。 图2.7 图2.8 根据余弦波与正弦波的幅度计算角度在一次把波形看作成三角形。如果增加三角形的宽度而高度保持不变角度会缩小。如果增加三角形的高度但宽度保持不变角度会增加。如果增加三角的宽度那么为了保持角度不变就必须增加它的高度。这意味着三角的角度、高度和宽度之间必须存在关系。保持三角形的宽度相同并增加高度会使角度一起增加。随着除法运算的分子增加除法结果也会增加如图2.9就像三角形的角度一样所以高度必须构成除法的分子。保持三角形的高度相同并增加宽度会使角度缩小。随着除法运算分母的增长除法结果会缩小如图3.0就像三角形的角度一样因此宽度必须构成除法的分母。
图2.9 图3.0 在图形的X轴上绘制三角形的高度除以宽度的结果要测量三角形的角度并将其绘制在Y轴上如图3.1。可以清晰的看到三角形的宽度高度和角度之间存在的关系但它不是直接甚至线性关系。也就是说高度除以宽度还不足以给出角度如果只知道三角形的的宽度和高度还需要进行其他操作才能恢复三角形的角度。 图3.1 为了计算三角形的高度我们将斜边乘以角度正弦。计算三角形的宽度将斜边乘以角的余弦。通过这俩个表达式来恢复角度。斜边乘以角的正弦除以斜边乘以角的余弦。经过化简又得到一个正切函数。所以恢复角度的公式为高度/宽度的反正切如图3.2所示。 图3.2 通过上述例子中正弦波幅度为29余弦波幅度32角度计算后得出结果42°。因此通过将频率相同但是幅度不同的正弦波和余弦波相加可以改变合成波的相位。 众所周知一个圆圈是360°然而还有一种角度测量方法使用弧度。一个弧度的定义是这三个长度都相等的角度如图3.3所示。因此如果圆的半径为1那么角度在一弧度时该弧度的长度也将为1。如果将角度增加到2弧度那么求得弧长需要知道圆弧的长度总是等于半径乘以弧度角。所以在单位圆上就是一个半径为1的圆弧长为2就像弧度角一样。我们知道圆的周长等于π3.1415926。如果弧度角等于圆弧的长度那么在一个圆上一定有2π弧度。 图3.3 欧拉公式 著名的美国物理学家理查德·费曼研发出一项特别的公式他称之为“我们的宝石”和“数学中最杰出的公式”。他所指就是如图3.4所示。是菜昂哈德·欧拉于1748年首次表达的欧拉身份。昂哈德·欧拉 于1707年4月15日出生瑞士巴塞尔。英国数学家数学家、天文学家、物理学家、力学和天体力学的奠基人之一。欧拉对数学物理做出了许多重大的贡献被公认为是十七世纪最伟大的数学家。欧拉在数学、几何学、代数学、力学和天文学等方面都有很深的造诣为近代科学奠定了坚实的基础。他提出了“连续统假设”证明了欧拉数创立了“复变函数论”。欧拉对数学物理的贡献有许多是后人难以企及的他被称为“数学王子”。 图3.4 雅各布·佰努利发现如果他在年初存入1美元并在当年给自己共100%的利息那么无论他全年分摊多少利息他都无法超过总计到年底为2.72美元。确切的数字实际上是2.71828182......是欧拉以他自己命名了这个数字“e”这个就是那个卡住无限数字的符号。 普通的计算器只有4中运算加减乘除。数字e实际上可以通过无限系列的加法、乘法和除法来计算如图3.5所示。约等于2.7182787698 图3.5 那么有乘法表示就像把e写成1的幂。更一般化点就把e提到x的幂那么无穷级数变成如图3.6所示。 图3.6 现在碰巧e到x并不是唯一可以通过这种无限级数计算的东西。三角函数正弦和余弦也可以用无穷级数来计算x的正弦等于x除以1阶乘 减去x平方3的阶乘加上x到5的阶乘减去x7的阶乘等等如图3.7所示。很容易预测这个系列的延续当分母每次增加2时我们增加的x的幂也的如此。每个项的正弦在加号和减号之间保持交替x的余弦也会发生类似的事情。与正弦不同是余弦是从2开始。正弦与余弦的无穷函数看起来与e非常相似。如果我们把余弦函数和正弦函数加在一起也没有办法得到欧拉函数。因为在正弦和余弦中项直接的符号不断变化而在计算欧拉数的系列中不会出现这样的问题。所以为了使方程相等我们需要乘以x平方的值必须等于负1。当x立方时也是如此当x再次被提升到6和7的幂时也是如此依此类推。当我们对一个数字进行平方时结果总是正数不是负数。所以这个时候出现了一个数学符号虚数i就这样诞生了i是唯一一个平方后得出负数结果得-1。 图3.7 如果取欧拉数并将其增加不是x的幂而是-1如图3.8所示。虚数乘以i的项用红色表示实数项用蓝色表示将它们从新排序后得到图3.9所示。其中所有虚数都乘以i因此可以将i进行分解。而实数项用来计算x的余弦值的项虚数项用来计算正弦值的项并将其全部乘以i。因此欧拉数与余弦和正弦之间只要将e增加到ix方并将正弦项相乘以i才能有效。所以这就是欧拉公式也叫做最美公式如图4.0所示。欧拉的恒等式如图4.1所示。 图3.8 图3.9 图4.0 图4.1 复数 在傅里叶变化中i起到了非常重要的作用。它将事物分开将信号分解为其组成正弦波的算法。在此还需要用到关于相位通过将频率相同但幅值不同的非相移余弦波和非相移的正弦波加在一起制作任何相移的正弦波。用数学公式表达如图4.2所示其中A是余弦波的幅度B为正弦波的幅度。 图4.2 将图4.0用一种新型的表示法复数。复数是由实数和虚数构成的数。实数是我们日常生活中常见的数例如1、2、3等等而虚数是-1的平方根通常用i表示。所以虚数i的平方是-1即i²-1。 复数的一般形式可以写为abi其中a和b都是实数a称为实部而b称为虚部。例如32i就是一个复数其中实部是3虚部是2。同样地我们也可以将实数看作是虚部为0的复数。例如5可以写为50i。 复数的定义使我们能够更好地描述和操作一些实际问题。例如在电路分析中复数经常用于描述电阻、电感和电容等元件的阻抗。在量子力学中复数被广泛应用于描述波函数和量子态。此外在信号处理和图像处理领域复数也被用来表示信号和图像的频谱。 复数之间可以进行加法、减法、乘法和除法运算这些运算与实数的运算规则类似。例如两个复数相加时只需将它们的实部和虚部分别相加即可。而两个复数相乘时则需要使用分配律和虚数i的平方等规则进行计算。 此外复数还有一些其他的特性。例如复数的共轭是指保持实部不变而将虚部取负的操作。如果一个复数的共轭与它本身相等那么这个复数就是实数。我们也可以通过求模来计算复数的绝对值即复数与原点的距离如图4.3所示。 图4.3 将三角形的底代表余弦波的振幅A。三角形的高度表示正弦波振幅B。这样会产生一个正弦波S振幅为5由三角形的斜边表示相移为53.1°由三角形的角度表示如图4.4所示。
图4.4 还有另为一种计算三角形底边长的办法是取角的余弦乘以斜边长度乘以53.1°的余弦等于3。而高度三角形的正弦值乘以斜边的长度5乘以53.1°的正弦等于4。所以可以这一种方式来描述复数点用它的极坐标形式来写5·cos53.1° 5i·sin53.1°。极坐标就是一个点极点是指向该点的线。可以通过它的长度和角度来描述这条线。将复数的极坐标与欧拉公式联系起来可以写成5·e^(i·53.1°)。 接下来我们来看复数的四种基本运算加法、减法、乘法和除法 1. 加法和减法: 加法和减法是最直接的按照实部与实部、虚部与虚部分别进行加减即可。例如(abi) (cdi) (ac) (bd)i(abi) - (cdi) (a-c) (b-d)i。 2. 乘法: 复数的乘法稍微复杂些需要按照乘法分配律进行计算。如(abi) * (cdi) ac adi bci - bd (ac-bd) (adbc)i。 3. 除法: 复数的除法最为复杂需要通过乘以共轭复数来进行。例如(abi) / (cdi) (abi) * (c-di) / ((cdi) * (c-di)) [(acbd) (bc-ad)i] / (c^2d^2)。 指数乘积规定对于任意实数a和b指数乘积规定可以表示为 a^b e^(b * ln(a)) 其中a^b表示a的b次幂e是自然对数的底数ln(a)是a的自然对数。 性质 指数乘积规定具有多个重要的性质 1. 根据指数乘积规定a^(bc)等于a^b乘以a^c。这个规则对于简化复杂的指数计算非常有用。 2. 指数乘积规定使得两个数的乘积可以转化为它们各自的指数函数之积。例如a*b可以表示为e^(ln(a) ln(b))。 3. 指数乘积规定还可以转化为对数运算。具体来说如果a^bc那么ln(a^b)等于b * ln(a)即ln(c)。 例如欧拉函数可以以指数的形式将它写出则完全相同的方法适用于abi·cdi。以34i·92i为例为此可以使用勾股定理和反正切原则。对34i3的平方加4的平方在开根号等于5。使用反正切acrtan4/3≈53.1°。对92i进行一样的操作如图4.5所示。 图4.5 将乘法改写成5e^(53.1°·i) · 9.2e^(12.5 °·i)得出结果46e^(65.6°·i)。接下来可以使用Polar形式将其转换回笛卡尔形式并表明这俩种方式给出的答案如图4.6所示。利用复数的乘法与除法即可算出结果0.410.35i。 图4.6 傅里叶级数 变量的任何函数无论是连续的还是不连续的都可以在变量的倍数的一系列正弦中展开。傅里叶将任何信号视为许多正弦波每个正弦波都有自己的频率相位和幅度全部加在一起。这就是转换一种从更简单的构建中构建信号的方式。但傅里叶绝不是唯一一个发明变化的人。还有许多其他类型的转换例如拉普拉斯变换就完美的克服了傅里叶的一大缺点。 拉普拉斯变换的定义如下 对于一个定义在非负实数集上的函数f(t)其拉普拉斯变换F(s)定义为积分形式即 F(s) ∫[0,∞) e^(-st) f(t) dt 其中s是复变量可以是实部和虚部都大于零的复数。拉普拉斯变换的本质是将时域函数通过指数衰减的方式映射到复平面上的一条曲线这样可以方便地进行频域分析和求解微分方程。 拉普拉斯变换的一个重要特性是线性性质。即对于任意两个函数f(t)和g(t)以及对应的拉普拉斯变换F(s)和G(s)有以下等式成立 a*f(t) b*g(t) ↔ a*F(s) b*G(s) 其中a和b是任意实数。这个特性使得拉普拉斯变换在信号处理和系统分析中非常有用可以方便地进行信号加权、线性组合和系统响应的求解。 另一个重要的特性是微分和积分的关系。对于一个函数f(t)的导数f(t)和拉普拉斯变换F(s)有以下等式成立 f(t) ↔ s*F(s) - f(0) 其中f(0)是函数f(t)在t0时刻的初值。这个特性在求解微分方程和系统响应时非常有用可以将微分方程转化为代数方程从而简化求解过程。 小波变换的定义如下首先它能够提供信号在时间和频率上的局部信息使我们能够更准确地分析信号的瞬时特征。与传统的傅里叶变换相比小波变换能够在时间和频率之间建立更好的平衡从而在处理非平稳信号时更为有效。 其次小波变换具有多分辨率分析的能力。它能够将信号分解成不同尺度的子信号使我们能够在不同频率范围内捕捉到信号的细节和特征如图4.7所示。 图4.7 卷积 信号是指随时间变化的物理量例如声音、图像、电压等。在数字信号处理中信号通常以离散的形式表示即在一系列离散时间点上取样得到的数值序列。 卷积是一种数学运算用于描述两个函数之间的关系。在信号处理中我们通常将信号表示为离散时间上的函数因此信号卷积可以看作是两个离散函数之间的运算。 信号卷积的定义如下给定两个离散函数f(n)和g(n)它们的卷积记作f(n)*g(n)定义为 f(n)*g(n) ∑[f(k)g(n-k)] 其中k为求和变量取值范围为整个定义域。换言之卷积是将两个函数按照一定规则进行加权求和得到一个新的函数。 卷积的工作原理我们可以想象一下当窗口滑动到某个位置时它会与该位置及周围的数据进行相乘并将相乘结果相加得到一个新的数值。这个新的数值就是卷积运算的输出值它可以被用于后续的处理或分析。 卷积的工作原理看似简单但却非常灵活。它可以通过调整窗口的大小、形状和滑动的步长来实现不同的效果。 为了找到具有我们在信号中寻找的频率的正弦波的相位我们可以做完全相同的事。在信号上滑动或卷积相同频率的余弦波。如图4.8所示称红色的余弦波为测试波。之所以不使用正弦波是因为使用傅里叶公式先计算余弦后计算正弦如图4.9所示。 图4.8 图4.9 所以测试波为一个余弦波通过增加它的相位慢慢的将它划过原始信号。当测试波与原始信号汇聚时就已知在原始信号中找到了正弦波的相位。但是电脑无法识别则操作所以需要生成一个数字一个表示两个信号匹配程度的分数。当它们完全不匹配时则分数很低。当它们完全匹配时分数很高。 现在将测试波上面的每个点与原始信号进行相乘并将相乘后的信号绘制到新的二维图上。在0度的角测试波值为1正弦波值为0.17所以两个点相乘得0.17,。以此类推如图5.0所示。 图5.0 然后测量由“乘法图”包围的区域会得到一个单一的值这个是测试波和原始信号正弦波匹配程度的分数如图5.1所示。随着测试波的平移分数也会随之改变。随着测试波越来越接近原始波形相同时乘法图所包围的区域会增加也就意味着分数也随之增加。当两个波相同时分数达到最高值180随着测试波划出分数再次下降。当两个波最异相时分数达到最低值-180。 图5.1 因此得出结论当测试波达到180°相位时两个信号完美匹配得分数处于最高值。可以在信号中找到了正弦波的相位。 接下来这是一个不规则的信号让我们首先将频率为1、幅度为1、的测试波与原始信号进行卷积如图5.2所示。 图5.2 当我们将测试波上的所有点乘以信号上的所有点并将其绘制在相乘图上相位是什么并不重要x轴上方的图所包围的面积总是与在它下面。这意味着当我们计算图形所包围的总面积时得到的分数为0分如图5.3所示。这就表明我们的原始信号中不包含频率为1的正弦波。在此基础上逐步增加使用频率为2的正弦波对原始信号进行卷积。可以看到随着相位的增长乘法图在上下移动分数也在随之上下移动如图5.4所示。因此原始信号中确定包含频率为2的正弦波。找到分数的最高点分数90在相位30°的时候。再次进行卷积这次使用频率为3的正弦波对原始信号进行乘法图上下移动分数也在上下移动。说明在原始频率中包含频率为3的正弦波分数在60°相位处达到最高值72因此正弦波相位为60°。使用频率为4正弦波进行卷积分数为0。得原始信号中不包含频率为4的正弦波。
图5.3 图5.4 所以每当尝试信号中不存在的频率时无论那个阶段分数为0当信号中确定存在频率时会得到一个非零的分数值。当找到分数最大点时就找到了正弦波的相位。为每个频率找到的最大分数正是正弦波的幅度成正比。 通过将分数除以信号周期的一半这个信号周期是360°所一半为180°对于这个测试波频率为3最高分数为72。所以72/1800.4这是实际信号中改频率的幅度。这个原始信号已经得出结论是由两个频率的正弦波组成的第一个频率为2、相位为30°、相对振幅为90。第二个频率为3、相位60°、相对振幅为72如图5.5所示。当两个正弦波加在一起时得到的结果看起来就想原始信号。因此使用卷积已经成功的将原始信号分解为正弦波的组成这就是傅里叶级数。 图5.5 任何相移的正弦波都可以通过添加一个频率相同但幅度不同的非相移余弦波和一个非相移正弦波来描述。所以不必执行整个卷积操作只需要先将我们的信号乘以余弦波然后乘以正弦波使用Phythagoras和反正切从两个分数计算幅度和相移每一波。 上述出现的第一个频率是2如果使用非相移余弦波作为测试波计算分数得到78则使用非相移正弦波测试波计算分数得到45使用勾股定理78平方根45平方根等于90。在使用反正切最终得30°与直接卷积算出来的答案一致。 卷积二 傅里叶变换将一个函数从时间域转换到频率域。在时间域中我们可以看到函数随着时间变化的情况而在频率域中我们可以看到函数中包含的不同频率成分。 傅里叶变换的数学表达式是 F(ω) ∫[f(t) * e^(-iωt)] dt 其中F(ω)是函数在频率域的表示f(t)是函数在时间域的表示e^(-iωt)是一个复指数函数ω是角频率。 傅里叶变换的工作原理可以通过以下步骤来理解 1. 首先我们选择一个函数f(t)它可能是一个信号或图像等。 2. 然后我们将这个函数表示为一系列复指数函数的和每个复指数函数具有不同的频率和幅度。 3. 这些复指数函数与函数f(t)进行内积运算以获得在频率域中的表示。 4. 最后我们得到一个频谱它显示了函数f(t)中包含的不同频率成分的强度和相位。 提取了一段来自小提琴的声音信号与余弦波相乘时得到的结果如图5.6所示。现在需要找到乘法信号图下的区域以获得表示该余弦波对整个信号的贡献幅度分数。想要计算一个矩形的面积可以用矩形的高度乘以宽度。所以可以将信号做成很多矩形划分成很多片把每一片当作一个矩形然后把所以矩形的面积加在一起如图5.7所示。 可以看到矩形块如何近似于相乘信号的形状。再次减少每个矩形直接的宽度并增加矩形的数量会得到一个更好的近似值。矩形越薄数量越多就越接近真实信号的面积值。 图5.6 图5.7 如果每个矩形无限薄并将矩形的数量增加到无穷大可以精确计算面积。这称为积分在数学上表示∫。因此为了找到乘积余弦图下的面积取乘积信号x的t乘以cos2πft相对于t的积分即每个切片的厚度为t秒其中一个t和下一个t之间的差异无限小。在数学上将其写为∫x(t)·cos(2πft)dt的积分积分的上下限在信号的一个周期就足够了因为此信号周期完成后该信号只会重复自己。在原始信号中没3.5毫米所有内容都会重复一次因此信号周期为3.5毫秒。因此积分的时间等于0和3.5毫秒之间如图5.8所示。 图5.8 现在由于这是一个重复的信号只要咋子恰好1个周期内进行积分哪里设置限制实际上并不重要。假设这个信号一直持续下去。信号会一直一直重复下去。也可以在-1.75毫秒到1.75毫秒进行积分最终结果依然是一样的只要积分的时间保持好3.5毫秒。只限于这个特殊的信号 可以通过上述公式计算出第一个频率的余弦波对原始信号的贡献程度分数。重复刚刚的所有步骤只是现在使用相同频率的正弦波。将原始信号乘以正弦波然后积分使用数学公式表达∫x(t)·sin(2πft)dt得到的该正弦波对原始信号贡献幅度的分数。 现在已经有两个独立的积分称之为余弦积分和正弦积分如图5.9所示。如果想将两个公式合二为一那么必须要用到函数i如果将正弦积分乘以i就可以将两项积分进行相加处理甚至将一个与另一个相减也不会对结果产生任何影响因为i将他们分开。 图5.9 在进行组合之前我们要知道傅里叶变化是一种将信号分解为其组成正弦波的方法。因此如果傅里叶知道如何分解信号他也可以将这些正弦波重新组合在一起以恢复信号并且确定可以使用傅里叶逆变换如图6.0所示。 图6.0 将两个信号的积分合并为一个不会改变结算结果。 因为i将两项式子分开不会改变结果t信号x乘以正弦项和余弦项可以将其作为公因式提出,并将其它所有内容放在括号中。在通过欧拉公式将括号内的数据进行转换如图6.1所示。这项公式诠释了在我们第一个测试频率下计算余弦和正弦项的分数。现在需要对许多不同频率一次又一次的运行这个方程可以通过不断改变“f”的值做到。 图6.1 公式中n是一个索引如图6.2所示当n等于1时我们用第一个测试频率F(1)测试我们的信号得到结果C(1)。其中C(1)是一个复数实部值代表余弦波在频率1的分数而虚部代表正弦波在频率1的分数。当n等于2时用下一个频率F(2)测试原始信号显然易见得到的结果也是C(2)以此类推。这是傅里叶级数原理。 图6.2 傅里叶变换 傅里叶级数非常适合模拟重复信号但重复信号不携带任何新信息。我们在日常生活中遇到的绝大多数信号都不会重复。因此对于那些则需要用到傅里叶变换。 下述方程如图6.3所示在表面上看确实相似但是它们之间有两个重要的区别。首先是积分的限制不同在傅里叶级数方程中积分的范围从信号周期的负一半到信号周期的正一半而在傅里叶变换方程中它是从负无穷大到正无穷大。其次傅里叶级数中的频率项F有一个小下标n这个下标n也会出现在方程的结果中。但是在傅里叶变换方程中没有这样的下标n所以方程的结果大不相同。 图6.3 积分是在图表下面积的地方在傅里叶级数中使用该值作为分数来告诉我某个频率对整体信号的分数有多大。积分的极限表示正在查看多少图。在傅里叶级数中只需要查看信号的一个周期因为在该周期完成后信号会重复。而在傅里叶变换中需要查看整个信号这样做的原因是想要对永不重复的信号进行建模因此没有循环所以需要从时间开始到时间结束。时间是无穷的所以需要查看从负无穷到正无穷的信号。 指数信号就像它的名字所暗示的那样继续以指数增长在负数部分可能很小单从某一刻开始当时间为零时它会随着时间迅速增大。所得傅里叶无法对此信号进行求解。傅里叶算法的初衷是用于查找信号中存在的正弦波如果将信号诚意给定频率的余弦波其次测量乘法余弦波信号图特定频率下的分数。在将原始信号与特定频率的正弦波进行相乘计算积分得出分数。以此类推直到找到信号中的所有频率。最终结果只会等于0 傅里叶的离散型 在上述中描述到项C(n)实际上是一个复数列表每个复数代表一个具有不同频率幅度和相位的下标n是正弦波。当f下标改变时也会出现相同的频率如图6.4所示。以此类推如果想要取该列表中的每个正弦波并将它们加在一起得到结果的波看起来就像我们的原始信号。傅里叶变换方程中的任何地方都没有这样的下标n因为与傅里叶级数给出的正弦波列表不同傅里叶变换为我们提供了一个连续的函数通过测试连续的频率列表产生的sin(x)由字母F表示。 图6.4 信号的带宽意味着从信号中的最小频率到最大频率每个频率都在某种程度上存在这就是为什么没有在傅里叶变换方程中任何地方都下标n因为根本无法为傅里叶变换表示的信号中的每个频率单独编号因为从最小频率到最大频率有无数个中频。 傅里叶变化 首先介绍一下卷积公式如图6.5所示。其中F(τ)是任意函数函数就是信号的另一种叫法。τ信号如图6.6所示。g给定一个余弦信号之所以选用余弦是应为余弦波反转后还是余弦波而正弦波反转后是-正弦波将两个信号绘制在同一个二维坐标系中X轴作为时间系随着时间的变化移动余弦信号并将两个信号相乘将结果绘制一个新的二维坐标系中。随着t的变化两个信号的乘积在x轴上上下下移动如图6.7所示。 图6.5 图6.6 图6.7 当乘积结果达到最大值时此时的分数会告知匹配程度如何。这就是傅里叶卷积傅里叶变化转换信号从时域进入频域任何信号都可以构建不同频率的正弦波而傅里叶变化告诉我们信号中存在那些频率。将结果得到的分数对应每个频率绘制图表将频率防止X轴Y轴防止对应频率的最高分数如图6.8这就是傅里叶变化对此信号。 图6.8 拿到余弦波与原始信号卷积后的的分数与正弦波与原始信号卷积后的的分数。如同一频率下余弦波为0.039正弦波为-0.119使用勾股定理公式可以计算出√0.039平方-0.119平方0.125幅度。而相位公式为acrtan-0.119/0.039。 快速傅里叶 随机获取一段信号将此信号进行傅里叶变换变换后的结果如图6.9所示。每行在该列表中表示不同的正弦波。正弦波有三个属性频率幅度和相位。在傅里叶变换开始时它并不知道原始信号是由那些正弦波构建的所以从理论上讲它需要测试无限数量的不同正弦波覆盖所有不同的可能性显然无法在计算机中做到这一点。因此必须更有选择性并且测试有限数量的正选曲线。 图6.9 频率所获取的N项结果全部在列表中K为每的索引N为列表中的项目总数R为采样率。计算每项的频率公式Fk/n·R。举个例子假设我想找到此列表中第二项的频率如上图所示它的索引K为1以采样率对信号进行采样R50Khz有1024个采样点代入公式F1/1024·50 得48.83Hz。现在知道了每一项的频率如何求得 复数的实部表示余弦分量虚部表示正弦分量。当然这些数字用三角形也可以表示。想象一下在一个圆圈内旋转的三角形随着三角形角度的增加它的高度增加一担角度通过90°高度再次开始降低当角度超过180°时高度开始增加但是方向相反一担角度超过270°高度就开始了再次减少一切都重新开始。如果将其绘制在带有角度的二维坐标系中X轴上是三角形的高度在y轴上只需要看它描述了那个功能如图7.0所示。 图7.0 证明正弦波和余弦波与三角形相连选择第3项用作计算它所代表的正弦曲线的幅度和相位。复数描述的余弦波和正弦波的振幅如果要求正弦曲线需要另做运算。 实例制作三角形的宽度等于复数的实部表示余弦波的振幅三角形的高度等于复数的虚部表示正弦波的振幅。将正弦波与余弦波进行相加得到正弦曲线通过勾股定理计算得出正弦曲线的幅度如图7.1所示。 图7.1 三角形的角度与其宽带及其高度增加高度角度增加反之亦然增加宽度角度较小反之亦然。但是三角形的高度和宽度不是线性的它描述的函数是逆函数切线函数。将三角形的高度除以宽度在进行反正切即可求得角度。 对称性 当运行完FFT计算后会注意到关于中心频率的对称因为信号现在是离散的。也就是说只知道信号在哪里位于某些离散时间点。意味着拥有的信号仅存在于某些离散时间点。 实例通过采用DFT进行实验简单的2hz余弦波该信号以16hz采样信号中有16个样本因此DFT将测试16种不同的正弦波间隔为1Hz。测试的第一个正弦波为0Hz1Hz第三个2Hz以此类推一直向上至15Hz对于每个频率测试信号用余弦波和反正弦波将相乘后的结果绘制在标图上计算分数查询每个频率与原始频率的匹配程度。使用勾股定理将余弦波的分数与正弦波的分数进行计算计算后的结果再次绘制在新的一张标图上如图7.2所示。对于第一个频率0Hz余弦与正弦得到的分数都为零重复过程1Hz会产生另一个零结果当测试2Hz频率时候余弦分数不为零意味着2Hz是存在的频率在原始信号中但是正弦分数依然为零代表余弦分数是为原始信号在2Hz频率上贡献所有能量。以此类推测试剩余其他频率.......。当用14H波测试信号时候余弦分数给出完全相同的分数就像在测试2Hz频率时一样。 图7.2 在14Hz时测试波正好占据采样值与2hz测试波相同关于两个波所能获取的信号就仅存在与采样点上而巧合14Hz与2Hz的那些采样点重合。没对频率都会发生同样的事情中心频率约为8Hz对于7Hz和9Hz频率也是如此等等......这就是为什么幅度输出的原因的FFT是对称的中间的频率是一个非常特殊的频率称为奈奎斯特频率。 DFT计算方法 如图7.3所示。输入不同的K数据从0一直到N-1个进行与测试波形进行乘积最后分别求和得出复数输入不同的K得到不同的结果。为了让式子看起来更加简介对负指数项做个定义如图7.4所示。称为旋转因子$W_N^k e^{-j2\pi k/N}$其中$N$为FFT长度$k$为旋转因子的下标。旋转因子可以预先计算并存储在一个数组中以便在计算FFT时使用。在蝴蝶运算中旋转因子的下标$k$与输入数据的下标$i$和$j$有关。将两者式子进行组合如图7.5所示。 图7.3
图7.4 图7.5 举例说明取长度N4的信号x(n),求得DFT结果列式如图7.6所示。把K0一直到K3全部列出来如图7.7所示。将旋转因子提出来可以写成矩阵形式如图7.8所示其中X为复数结果x为输入的原始信号与旋转矩形相乘。 1.原时间序列变为频率序列每个点均对应一个频率分量 2.X(k)中每个点均为复数X(k)与X(N-k)为共轭状态 3.摸值和相位代表当前点对应频率分量的幅度和相位余弦函数形式 4.有效频率分量为N/2个点相邻点的频率间隔为F(s)/N。 图7.6 图7.7 图7.8 式列信号为f(t)cos(2π·3·tπ/3)以16Hz频率采样即3个信号周期采集16个点如图7.9所示。已知该信号频率3Hz初始相位π/3幅度为1。使用旋转因子矩阵乘以原始信号的点数等到复数结果如图8.0所示。因此采样了16个点所以N16则旋转因子为W下标16上标nk则等于e^(-j2π/16)nk图8.0则是将公共部分e^(-j2π/16)提出在第一行K0时n从0-15变化取值相乘此行结果为直流信号将所有的输入信号进行相加以此类推最后一行k15n从0-15变化。计算结束后得到16个复数如图8.1所示。可以看到X(13)与X(3)是共轭的状态则有效点数只有8个点。 图7.9 图8.0 图8.1 旋转因子 以8个点单位圆为例如图8.2所示一个复数的平面横坐标为实数轴纵坐标为虚数轴ω的取值位于圆上的这些点。将这些点中的数按照一定规格进行排布就变成了如图7.8那样的矩阵形式。用公式的方式表示如同图7.3其中∑e^(-i2Π/N·nk)可以用矩阵表达式x[n]为原始信号x[k]为每个相乘后得到的结果如图8.3。所描述的是一个基矩阵与一个原始向量相乘得到一个新向量。 图8.2 图8.3 降低采样复杂度 原始信号yx^52x^4x^3x^2x0.2函数图像如图8.4所示可以将多项式进行奇偶划分但是会发现一个问题多项式最高次的次数没有改变导致划分后的多项式再次进行采点需要满足最高次的次数采点数还是一样多所以需要对多项式进行将此操作将XX^2那么整个式子的次数减半用同式表示如图8.5所示。 图8.4 图8.5 复共轭 以一个简单的5Hz正弦曲线为例如图8.6所示如果将它分解成它的余弦波和正弦波它们都有5Hz的频率但是它们的振幅不相同当当它们加在一起时它们会产生5Hz的正弦信号。如果对原始信号进行傅里叶变换将其转换为频域余弦分量如图8.7所示这是5Hz时的峰值在时域中余弦波振幅有0.8但是在频域中只有0.4还有另一个-5Hz时相同振幅的峰值。频域状态下的正弦分量也是一样效果如图8.8所示。在5Hz频率幅度是一半但在-5Hz振幅为负时域振幅一半。在时域中余弦和正弦表示每个频率的分量由一个复数余弦分量为实数部分正弦分量为虚数部分。 图8.6 图8.7 图8.8 此信号中有两个非零频率5Hz和-5Hz所以这个信号可以是由两个复数表示其中5Hz表示为(0.40.3i)。其中-5Hz表示为(0.4-0.3i)为了恢复信号从频域回到时域则需要执行傅里叶逆运算执行后会最终得到一个列表需要相加的复数其中每个复数代表贡献信号的一个特定频率这些数字中大多数将为零保存正负5Hz频率因为它们是唯一的频率存在与原始信号中如图8.9所示。把两个数字加在一起时0.47i和-0.47i抵消了这意味着任何时刻的信号幅度都可以用单个实数来描述。信号是真实的这就是为什么幅度谱任何真实信号的傅里叶变换总是沿着0Hz线反射以及为什么任何真实信号的相谱始终是沿着0Hz反射并反转因为如果不这样当加在一起时虚构的部分不会抵消并且时域信号不是真正的信号。 图8.9 傅里叶逆变换 Ω对我们考虑复平面时很有用复平面组合余弦并将信号的正弦分量输入一个二位坐标系X轴表示振幅任何时刻的余弦分量Y轴表示振幅任何时刻的正弦分量。正弦波和余弦波的频率可以认为是旋转的速度圆圈中间的一条线这是这条线的角速度如果使用Z轴绘制角度那么这条线然后圆就变成了一个螺旋如图9.0所示。生产线向上移动的速度Z轴是它的角速度然而正如线条可以沿着Z轴向上移动一样它也可以向相反的方向移动下。Z轴如果从两维角度看这相当于以另一种方式旋转当线沿Z轴向上移动是其角度速度或其频率为正当线沿Z轴向下移动时其角速度或其频率为负这就是负频率。 图9.0 以图7.0举例如果三角形以相反的方向旋转那么余弦波看起来是一样的单正弦图却是相反的在正旋转期间现在为负在负旋转期间反之亦然如图9.1所示。了解了其中的负频率就解释了为什么余弦的振幅在负5Hz的时候依然是正的。以及为什么正弦波的振幅在负5Hz的时候分量为负。 图9.1 如果想恢复信号需要利用傅里叶立方体如图9.2所示可以帮助我们快速理解傅里叶逆变换我们的信号在傅里叶变换中信号向下进入其组成正弦波而傅里叶逆变换恰恰相反将所有正弦曲线重新组合在一起以恢复原始时域信号。在不同频率下相加恢复原始信号首先需要从余弦分量和正弦分量中恢复成正弦波和余弦波公式表达如图9.3所示。其中X(f)是复数的列表表示振幅各种余弦分量与正弦分量频率每个复数乘以余弦波和正弦波在该频率下。但是两个复数相乘不足以恢复成我们需要正弦波和余弦波因为乘法的结果依然包含一个虚部非常复杂看起来像三维螺旋。 图9.2 图9.3 逆向的最后阶段傅里叶变换需要加法将所有复杂的正弦曲线相加在一起复数的正弦分量负频率与正弦分量完全相反对应的正频率将它们相加在一起时虚轴的正弦分量相减低只留下了一个真正的正弦波。所以这就是负频率为什么存在于我们的信号中它的存在是为了确定虚部傅里叶变换后进行逆变换。 FFT 回顾离散傅里叶DFT公式如图9.4所示从第一个式子可以看到在进行每一个点的DFT运算的时候要将e和-j^(2π/N)nk进行相乘在相加才能得到X(k)的结果。每一点的X(k)的值需要进行N次相乘和N-1次的复数相加这样才能得到1个点的X(k)的值。要完成全部DFT运算需要进行N^2次复数乘法和N(N-1)次复数加法。 图9.4 为了提高运算速度先从公式进行改进为了方便计算将-j^(2π/N)带指数的表达为W^nk,将X(K)和x(n)之间用矩阵的形式来表达如图9.5所示。为了更加简化运算前辈了研究出W^nk两大特点周期性和对称性。 1.W^nk周期性具体来说周期性函数W^nk的周期为2π/n。也就是说在每个2π/n的间隔内函数W^nk会重复其值。这意味着如果你将W^nk的自变量通常是角度增加2π/n函数的值将再次与初始值相等。但是需要注意的是周期性函数W^nk的周期与n的取值有关。不同的n将导致不同的周期长度。例如当n1时W^nk的周期为2π当n2时W^nk的周期为π当n3时W^nk的周期为2π/3依此类推。 2.W^nk对称性对于k 0, 1, 2, …, n-1W^nk中的每个元素都是复数单位模的n次根。因此W^nk在固定的k值下具有周期性对称性即当k增加n倍时W^nk保持不变。当n为偶数时W^nk中的元素具有共轭对称性。具体地说当k取值为k1和k2时W^nk1中的每个元素与W^nk2中对应位置的元素互为共轭。 图9.5 因为再W矩阵中存在这很多不必要的重复运算所以这两大特点用在W的矩阵中。将N点DFT运算分解成两组N/2点的DFT运算然后取和。假设N是2的整数次方将公式分成两个部分一部分是偶数数列一部分是奇数序列。用(2r)表示偶数数列(2r-1)表示奇数序列如图9.6所示。 图9.6 实列以N2^24为列的DFT的运算进行说明公式如图9.7所示。根据此表达式可以列出每个K值时候的式子如图9.8所示。化作为矩阵形式如图图9.5所示首先根据W的周期序列长度为N的W的0次幂等于1序列长度为N的W的mN次幂等于1换算后如图9.9所示。在根据对称性序列长度为N的mNN/2次幂等于-1简化后的矩阵如图10.1所示。再次根据对称性N的rN/2等于-W的r次幂。
如图10.2所示。 图9.7 图9.8 图9.9 图10.1 图10.2 将矩阵式子的第二行和第三行互换第二列和第三列互换x(1)和x(2)互换矩阵等式不变如图10.3所示。使用红色线将W矩阵分成4个小矩阵左边上下完全相等的矩阵只和x(0)x(2)有关而右边的两列只和x(1)x(3)有关并成相反的关系。与x(0)x(1)x(2)x(3)对应相乘得出每一项的DFT结果如图10.4所示列出后可以发现很多相同的项比如X(0)和X(2)它们的前两项都是x(0)x(2)此项可以定义为G(0)它们的后两项都有一个共同的x(1)x(3)可以定义为H(0)。再次观察在X(1)和X(3)时将前两项相同值x(0)-x(2)定义为G(1),后两项相同的x(1)-x(3)定义为H(1)。定义后如图10.5所示。 图10.3 图10.4 图10.5 写成这种形式计算量得到了大大的提升因为其中G(0)G(1),H(0),H(1)是通过四次加法得到的X(0)X(1)X(2)X(3)通过一次加法就得到了最终的结果。实际在进行这样的计算中加法通过8次乘法通过4次就把结果得出。 总结一下N点的快速计算方法首先将N点的DFT运算分解成两组N/2点的DFT运算然后取和。通过图9.6可以继续延续分解已分解成偶数n和奇数n相加。将N点写成N/2点如图10.6所示。将偶数部分定位G(k)将奇数部分定义成H(k)则最后的公式简为如图10.7所示。其中G(k)相当于N/2点的DFTH(k)相当于N/2点的DFT。 图10.6 图10.7 如果用G(k)和H(k)来表达全部的X(k)那么应该并用G(k)和H(k)的两个重复的周期。由周期性可知G(kN/2)G(k),H(kN/2)H(k)在利用加权系数Wn的特点序列长度为N的kN/2等于Wn的N/2乘Wn的k也就等于-Wn的k。将这些表达式代入前面的X(k)的计算中可以由G(k)和H(k)来决定X(k)的全部关系式。 在图10.7中已得出公式那么根据对称性X(kN/2)等于G(k)-Wn的k乘以H(k)如图10.8所示其中X(kN/2)就是另外主值周期还有N/2点的X(k)。如果是4点的话那么只要k2 X(2)G(0)-H(0)和 k3 X(3) G(1) - W4的1乘以H(1)就可以计算出来了。 图10.8 蝶形运算 在N4为例DFT分组的实际情况如图10.9所示x(0)和x(2)分在偶数组x(1)和x(3)分在奇数组计算的结果又混在了一起了然后进行N点的DFT。 图10.9 也可以使用蝴蝶图在描述计算过程两条线的汇合点是表示两个数值的相加线旁边标注的W是与复数做乘法的运算如图11.1所示。此蝴蝶运算图包括两次乘法和两次加法以此类推将x(1)和x(3)全部描述出来。然而计算X(0)和X(2)也同样需要两次加法和两次乘法如图11.2所示。 图11.1 图11.2 此蝶形运算中有些地方重发可以进行简化处理。H(0)乘以W4的0和-W4的0可以改写成H(0)和W4的0相乘然后在分别相加分别相减这样可以使运算量减少到只有一次复数乘法和两次复数加法。那么就可以简化为如图11.3所示。因此就能得到这样的结论G(k)和H(k)获得的X(k)的过程中可以看到的蝶形单元两个也就是N/8个运算单元每个新单元进行一次复数乘法和两次复数加减法因此也就是需要进行N/2的复数乘法和N此的复数加减法。 图11.3 左半边是按N的奇偶分别组合了两个N/2点的DFT的运算也是分成了两个蝶形单元这个两个蝶形单元的运算量也还是N/2的复数乘法和N此的复数加法。对于4次的这种只在左半边进行了奇偶分解。比如但N8的时候逐级分解为两个N/2点四个N/2点的DFT运算如图11.4所示。 图11.4 FFT算法特点基本运算单元为一个蝶形第M级的蝶形如图11.5所示前一块的入端有上节点和下节点那么出端同样也有上节点和下节点。入的上节点与下节点相加减最后得到出端的上节点和下节点这样就构成了基本的蝶形单元。 图11.5 每一个蝶形单元都是独立的而每一级有N/2个蝶形在图11.4中在最左则和中则还有最右则中分别有四个蝶形。那么在N2^n的任意情况有M级的蝶形在每一级里可以分成N/2^i组。每个蝶形单元仅需一次复数乘法两次复数加法就可以得到输出端如图11.6所示这两点构成了一个蝶形单元并且这俩点不再参与别的蝶形单元的运算。 图11.6 还有一个特点在进行快速傅里叶变换时进行了一个码位的倒置输入如果是岛序的话则输出是正序因为在进行矩阵变换的时候在N4点的运算过程中将第二行和第三行进行了互换则相应的就把第二列和第三列也进行了互换。所最后得到的DFT的结果可以按二进制排序输入二进制而输出则是二进制的码位倒置如图11.7所示。 图11.7 逆变换IFFT可以用同一个公式只需某个形式参量代入不同的实际值如图11.8所示它们的差别只在于W的-nk还有运算的结果要乘以1/N。在应用FFT中还需一些注意事项 1.信号离散时采样率要攒足奈奎斯特频率。 2.N一定是2的正数词幂若不是要补若干个零凑成2的正数次幂。 3.数据长度要取得足够长▲f1/NT。其中▲f为频率分辨率。NT是数据的实际长度N是序列的点数Ts是采样周期。 图11.8 FPGA实现 实验说明 假设信号频率为F采样频率为Fs采样点数为N那么FFT之后结果就是N个复数每一个复数对应着一个频率点频率间隔(分辨率)为Fs/N第一个点表示直流分量(频率为0Hz)第N个点的频率为Fn(n-1)*Fs/N; 各频率点下复数的模值就是该频率下正弦信号的幅度特性FFT结果具有对称性通常我们只使用前半部的结果。设Xn[4 3 2 6 7 8 9 0]经过FFT后得 39.0000 -10.77826.2929i 0-5.0000i 4.7782-7.7071i 5.0000 4.77827.7071i 05.0000i -10.7782-6.2929i 假设采样频率为48KHz那么两个频点之间的间隔48/86Khz。第一个为直流分量0Hz第二个点为6Khz第二个点频率为12KHz以此类推。对于中间点两边的值实部相等虚部相反称为共轭对称但是它们的模值相等也就是幅度相等反映在频谱上的效果就是会出现两个高度相同的脉冲并对称的。 实验实现 时序图如图11.9所示。当source_sop拉高表示一帧数据开始source_eop表示一帧数据结束。sink_real 数据的实部sink_imag数据的虚部。source_real和source_imag同样为输出后结果的实部和虚部数据当source_valid拉高时候数据有效。 图11.9 IP核搭建 第一个选项是同时进行几路数据流并行。 第二个选项是变换的实际点数如果设计可以在线更改变换点数的FFT这里选择的点数是所需要的最大的FFT点数。 第三个工作时钟Target Clock Frequency。 第四个选择一种FFT结构包括流水线Streaming、基4 Burst、基2 Burst和轻量级基2 Burst计算速度和消耗的资源依次减少。 第五个需要实时更改FFT的点数则需要更改Run Time Configurable Transform Length。
Implementation标签卡下 设置FFT的数据格式为定点Fixed Point或浮点Float Point 设置输入数据的位宽和相位因子位宽相当于旋转因子。 ARESETn这个信号一般要选上 “Output Ordering”设置FFT计算结果以自然顺序Nature Order或位/数值反序Bit/Digit Reversed Order输出。 代码编写
module FFT_TOP(input wire clk,input wire rst);reg [3:0] state;reg [10:0] cnt;reg [15:0] real_data;reg [15:0] imag_data;reg s_axis_data_tvalid;reg s_axis_data_tlast;reg s_axis_config_tvalid;reg m_axis_data_tready;reg [7:0] s_axis_config_tdata;wire s_axis_config_tready;wire [31:0] s_axis_data_tdata;wire s_axis_data_tready;wire [31 : 0] m_axis_data_tdata;wire m_axis_data_tvalid;wire m_axis_data_tlast;wire event_frame_started;wire event_tlast_unexpected;wire event_tlast_missing;wire event_status_channel_halt;wire event_data_in_channel_halt;wire event_data_out_channel_halt;wire [15:0] real_dataout;wire [15:0] imag_dataout;assign s_axis_data_tdata {imag_data,real_data};assign real_dataout m_axis_data_tdata[31:16];assign imag_dataout m_axis_data_tdata[15:0];always (posedge clk ) begin if (rst)begin state 4d0;imag_data 16d0;real_data 16d0;s_axis_data_tvalid1b0;cnt 11d0;s_axis_config_tvalid1b0;m_axis_data_tready1b0;s_axis_config_tdata8d0;end else case (state)0:begin state state 1b1;imag_data 16d0;real_data 16d0;cnt 11d0;s_axis_data_tvalid1b0;s_axis_config_tvalid1b1;m_axis_data_tready1b1;s_axis_config_tdata8d1;end 1:begin if (cnt 1000-1)begin cnt 11d0;state state 1b1;end else begin cnt cnt 1b1;state state;end end 2:begin real_data 16h4; s_axis_data_tvalid 1b1; state state 1b1; end 3:begin real_data 16h3; s_axis_data_tvalid 1b1; state state 1b1; end 4:begin real_data 16h2; s_axis_data_tvalid 1b1; state state 1b1; end 5:begin real_data 16h6; s_axis_data_tvalid 1b1; state state 1b1; end 6:begin real_data 16h7; s_axis_data_tvalid 1b1; state state 1b1; end 7:begin real_data 16h8; s_axis_data_tvalid 1b1; state state 1b1; end 8:begin real_data 16h9; s_axis_data_tvalid 1b1; state state 1b1; end 9:begin real_data 16h0; s_axis_data_tvalid 1b1; state state 1b1; end 10:begin s_axis_data_tvalid1b0;real_data16d0;statestate;end default : ;endcase end always (posedge clk ) begin if (rst)s_axis_data_tlast 1b0;else if (state 9)s_axis_data_tlast 1b1;else s_axis_data_tlast 1b0;end xfft_0 your_instance_name (.aclk (clk), // input wire aclk//低位配置工作模式可同时配置多通道1为FFT0为IFFT[7 : 0].s_axis_config_tdata(s_axis_config_tdata), //配置信息有效信号无特殊要求置1即可。.s_axis_config_tvalid(s_axis_config_tvalid), //准备接收配置信息.s_axis_config_tready(s_axis_config_tready), // [31 : 0]输入数据高n位为虚部信号低n位为实部信号.s_axis_data_tdata(s_axis_data_tdata), //输入有效信号.s_axis_data_tvalid(s_axis_data_tvalid), //输出信号代表该ip核可以进行下一次运算.s_axis_data_tready(s_axis_data_tready), //输入最后一个数据时拉高信号只是为了区分数据来源属于可选配置。.s_axis_data_tlast(s_axis_data_tlast), //携带BLK_EXP以及OVFLO信息两个信息不能同时通过该输出口显示[31 : 0].m_axis_data_tdata(m_axis_data_tdata), //状态信息有效信号.m_axis_data_tvalid(m_axis_data_tvalid), // output wire m_axis_data_tvalid//非必选外部系统准备好接收状态信号的时候输入高电平。仅用于非实时状态。.m_axis_data_tready(m_axis_data_tready), // input wire m_axis_data_tready//输出数据的最后一位。.m_axis_data_tlast(m_axis_data_tlast), // output wire m_axis_data_tlast//拉高一个时钟周期指示开启新的处理过程一般用于更新配置信息。.event_frame_started(event_frame_started), // output wire event_frame_started//当ip核没有采集到足够点数就出现了.event_tlast_unexpected(event_tlast_unexpected), // output wire event_tlast_unexpected//当采集的数据足够但s_axis_data_tlast未拉高时会产生提示。.event_tlast_missing(event_tlast_missing), // output wire event_tlast_missing//当内核准备好输出状态数据但是无法输出。仅在非实时模式。.event_status_channel_halt(event_status_channel_halt), // output wire event_status_channel_halt//当内核准备好接收数据但无有效数据输入时产生提示。.event_data_in_channel_halt(event_data_in_channel_halt), // output wire event_data_in_channel_halt//当内核准备好输出数据但是无法输出。仅在非实时模式。.event_data_out_channel_halt(event_data_out_channel_halt) // output wire event_data_out_channel_halt
);endmodule仿真验证 与上述实验结果相同但是顺序不同。问题博主查询找............. 致谢/参考文献
Bracewell, R. N. (2000). The Fourier transform and its applications (3rd ed.). New York: McGraw-Hill.Oppenheim, A. V., Schafer, R. W. (2010). Discrete-time signal processing (3rd ed.). Upper Saddle River, NJ: Prentice Hall.Strang, G. (1993). Introduction to Applied Mathematics. Wellesley, MA: Wellesley-Cambridge Press.Proakis, J. G., Manolakis, D. G. (2006). Digital signal processing: Principles, algorithms, and applications (4th ed.). Upper Saddle River, NJ: Prentice Hall.Papoulis, A. (1977). Signal analysis. New York: McGraw-HillFourier, J. B. J. (1822). Théorie analytique de la chaleur.Shannon, C. E. (1949). Communication in the presence of noise.Cooley, J., Tukey, J. (1965). An algorithm for the machine calculation of complex Fourier series.Oppenheim, A. V., Schafer, R.W. (1968). Digital signal analysis.Mallat, S. G. (1989). A theory for multiresolution signal decomposition: The wavelet representation. https://blog.csdn.net/weixin_41920053/article/details/104618551 https://www.cnblogs.com/limanjihe/p/9999526.html https://blog.csdn.net/CLL_caicai/article/details/107022319 https://blog.csdn.net/weixin_42837669/article/details/118496144 https://blog.csdn.net/weixin_40640020/article/details/85006846 https://zhuanlan.zhihu.com/p/134150702 https://blog.csdn.net/beauytlife_1985/article/details/82115522 https://blog.csdn.net/a18271623106/article/details/115752709 一小时学会快速傅里叶变换_哔哩哔哩_bilibili Mark Newman - YouTube
