精华内容
下载资源
问答
  • 傅里叶光谱原理
    千次阅读
    2018-06-21 12:48:18

    傅里叶变换光谱仪(FTS)中的相干性和能量关系


    【1】

    在量子力学层面上,最小光谱间隔取决于原子光谱自然线宽(不考虑Doppler展宽、压力展宽、仪器展宽等效应),更本质上又取决于能量-时间不确定性原理,没有一个仪器能够分辨比自然线宽更窄的谱线。因此一般自然辐射光谱对仪器而言可视为一连续信号,其光谱间隔极其狭窄近似趋近于零,从而使信号Fourier变换谱的Nyquist本征频带趋于无穷大,这个量名义上称为频率是从Fourier频谱分析角度而言的,其实际对应物理量是时间长度[s]或空间长度[m]。在FTS将入射光谱辐射调制为干涉图过程中,由于自然线宽趋于零,对应时间域/空间域的Nyquist本征频带趋于无穷大。设光谱自然线宽为Δf [Hz],则对应Nyquist本征带宽为τ=1/Δf [s],这个τ又称为相干时间,进一步称L=c×τ为相干长度,二者分别对应于上述光谱在时间域/空间域调制信号(干涉图)的延展长度,亦即极限光程差LPD。对于Michelson干涉仪,其测量的双边干涉图最大光程差总长不超过相干长度L,而实际仪器设计的±MPD光程差范围也只是LPD内的一段距离。例如,对某一单色激光,假设其线宽为Δf=10kHz,对应相干时间为τ=0.1ms,相干长度L≈30000m,这是实际仪器干涉臂所远远达不到的,所以实际仪器设计也无需担心光谱线宽极限对光程差设计的限制。

    然而,在计算机FTS模拟器中,入射光谱不可能为连续信号,而是以dwavnum极密间隔取样的离散光谱信号(近似视为连续光谱),于是便引入了一个有限大小的Nyquist本征频带SPCFs=1/dwavnum,它规定了干涉图的极限光程差LPD=SPCFs。因此,在模拟干涉图采样过程中,设计的仪器(单边)最大光程差MPD不得超过LPD/2。实际上,在“输入光谱-调制干涉图”模拟过程中,由于Fourier变换的周期性,调制干涉图会以LPD为周期进行长程周期延拓,并且在±LPD/2附近产生较强的混叠效应,因此±MPD设计值应尽可能地远离±LPD/2。例如,在计算机模拟中,假设细密光谱间隔dwavnum=1cm-1,则LPD=1cm,最大光程差±MPD设计极限为±0.5cm,实际设计为±0.4cm为宜(对应标称光谱分辨率为1.25 cm-1);若随意增加MPD(比如2cm),干涉图会以LPD=1cm为周期重复出现,而光谱分辨率并不会同步增加(至0.25 cm-1),最高只能为1cm-1,实验显示确实如此。


    【2】

    在Michelson干涉仪将入射光谱辐射调制成干涉条纹过程中,时变均匀的入射光束转变为时变非均匀的相干光束。其能量分配关系为:对入射辐亮度光谱Lin(v) [W·m-2·sr-1·(cm-1)-1]按光谱积分得到单位面积单位立体角内的辐射光功率Pin [W·m-2],即干涉仪入瞳接收的光功率,在每个动镜移动位置(或每个光程差观测时间)对应这样一份入射光功率。然而,由于干涉仪调制的干涉光是时间分布非均匀(对应相干时间)或空间分布非均匀的(对应相干长度或光程差),亦即由分束器分出的两路子光束在两干涉臂的行程不同,使得随后在分束器处再次叠加的两路子光束并非来自同一份入射光,而是分别来自时间相继的两份入射光,这一时间差(光程差)反映Michelson干涉仪这类仪器具有时间调制性质(Michelson干涉仪属于时间调制干涉仪),于是在某一时刻或某一位置测量的相干光功率并不一定等于入射光功率。而零光程差是一个特殊位置,此时两路子光束在各自干涉臂的行程相同,在分束器处再次叠加时仍属于同一份入射光,进一步测量的相干光功率也就等于入射光功率(实际上,零光程差反映了双光束干涉的同时性)。为了通过干涉图间接测量入射光功率,可以对(相干时间内一段时间生成的)干涉图Iigm[W·m-2]沿光程差积分,得到的即是这段时间内相干光束的总能量Engigm [(W·m-2)·cm],亦是均匀入射光束对应的总能量Engin,将总能量在总光程差上平均,得到的即是这段时间内相干光束的平均功率 ,亦即均匀入射光功率Pin。由这段分析不难发现,Michelson干涉仪光程差虽然是一长度量[cm],但却具有时间意义,沿光程差积分是一种时间积分。

    此外,由于Michelson干涉仪分束器的半透半反性质,只有一半入射光功率参与相干叠加并被探测器探测,另一半虽然也参与相干叠加,但却沿原光路反向传播而不能被探测器探测,所以相干光平均功率等于入射光功率的一半。

    延伸思考:所谓测量功率实质是一统计量,假若能精细控制入射光源,使光子一个接一个地进入干涉仪,且每一测量过程中只有一个光子在干涉仪内,那么随着动镜移动会有干涉条纹出现吗?对于Michelson干涉仪这种所谓“分振幅”干涉装置,单个光子如何相干?


    【3】

    FTS调制的干涉光束入射到光电探测器上,由于探测器平方律响应特性(光电效应),一份光能量(功率单位)转化为电子量(电流/电压单位)。因此,对电学干涉图通过Fourier变换计算的光谱是为幅度谱S(v),虽然它具有电压/电流单位[Elc.V/cm-1],但却与光学功率谱等价( ),于是对幅度谱按光谱的积分量∑[S(v)]等于相干光平均功率 ,等于电学干涉图零光程差幅度值Iigm-peak,等于干涉仪入射均匀光功率的一半Pin/2,这即是仪器输出电信号与待测光谱辐射之间的线性关系,也即是光谱辐射定标的物理基础。另一方面,单纯就电信号而言,电学干涉图总能量或平均功率(电压信号平方和∑[Iigm2(cm)])与电学能谱总能量或平均功率(幅度谱平方和∑[S2(v)])相同(Parseval定理),反映的是电信号在Fourier变换前后的能量守恒。

    更多相关内容
  • #资源达人分享计划#
  • 叙述了常用的Mertz法的基本原理并进一步讨论了其中的计算效率问题, 通过利用实序列离散傅里叶变换的性质与相位校正的具体的处理内容相结合, 优化了Mertz法的计算处理, 提高了计算效率。
  • 傅里叶变换光谱仪实验实验目的1、加深对傅里叶变换光谱测量原理的掌握2、傅里叶变换光谱仪的特点及应用3、学会用傅里叶变换光谱仪测量光谱实验原理傅里叶变换光谱仪的核心是迈克尔逊干涉装置,如图1,P是光电倍
  • 阐述了傅立叶变换红外光谱仪的原理,并对红外光谱 的特点进行了归纳和总结。希望对读者有用。
  • 傅里叶变换光谱的基本原理出发,得出傅里叶变换光谱仪内部的双光束干涉的光程差是影响其分辨率的最主要因素。介绍了几种通过改变仪器的结构设计,来增加双光束干涉的光程差,从而提高仪器分辨率的方法。在分析聚焦...
  • 在干涉成像光谱仪的光谱复原中,由于系统存在各种误差,若直接对所得干涉图进行傅里叶变换重构,得到的光谱图会产生较大误差甚至错误。介绍了 Sagnac型干涉成像光谱仪基本原理,针对上述问题得到一套对采集得到的干涉图...
  • 为了明确像场调制傅里叶变换成像光谱仪的工作机理,通过分析多级微反射镜对成像光场的相位调制特性,建立了像场调制干涉成像的理论模型。数值计算结果表明,通过对获得的干涉图像数据立方体进行图像剪切与图像拼接,可以...
  • 首先对IRFPA非均匀性校正的原理与算法进行了简要阐述,在此基础上,通过对时间调制型、空间调制型及时间空间联合调制型傅里叶变换红外成像光谱原理结构的深入分析,针对这三种制式的光谱仪,分别研究出了相应的非均匀...
  • 为了获得更高精度的目标光谱偏振特性,结合现有的标准探测器技术与偏振探测技术原理,研制了一种光谱偏振分析仪。采用连续等角度间隔旋转检偏器对目标信号进行强度调制,通过傅里叶级数算法能够显著降低检偏器定位...
  • 提出了一种基于多级微反射镜的静态化新型红外傅里叶变换成像光谱仪结构。系统不含狭缝和可动部件,因此光通量大、结构稳定。介绍了该成像光谱仪的工作原理和光程差的产生方式。根据系统原理对后置成像光学系统进行了...
  • 以某微生物孢子气溶胶为对象开展探测原理试验,研究了生物气溶胶红外光谱特性。给出了生物气溶胶室内吸收池试验结果,用辐射传输模型解释了试验现象。通过试验验证了探测模型的合理性、有效性;结合探测模型和室内试验...
  • 定标原理如下:将测量目标的温度区间分为n个子区间, 测量并记录目标温度区间内n+1个不同温度黑体对应的红外光谱辐射计输出数据, 并分别计算各个子区间的定标系数; 进行红外光谱辐射测量时, 比对红外光谱辐射计输出...
  • 在以多级微反射镜为核心器件的静态傅里叶变换红外光谱仪中,基于传统单孔径透镜组结构的缩束系统由于口径较大、透镜数量较多和整体尺寸较长等原因,其体积和重量在整个光谱仪中占有很大的比重,严重制约了光谱仪的微...
  • 给出了傅里叶变换光谱测量方法的设计原理、硬件系统构成以及LabWindows/CVI软件数据采集系统,包括实现干涉图数据采集、实时显示、数据分析处理和光谱分辨率的精确计算等功能。实验结果证明该系统满足设计需求。
  • 该方法通过双强度调制干涉信号的和差处理,即保证单强度调制偏振光谱测量的优点,又能有效消除混叠现象,提高有效光程差和光谱分辨率,再对和差处理后的干涉信号进行傅里叶变换即可得到被测偏振光谱。文章对双强度调制...
  • 该方法不仅保留了现有PEM的优点,而且调制光电流的频率降低了2个,类似于3个数量级(10个类似于500 Hz),并且可以通过普通阵列检测器进行检测,并且入射光谱可以通过对调制信号中的低频分量进行傅里叶-贝塞尔变换...
  • 该方法通过双强度调制干涉信号的和差处理,即保证单强度调制偏振光谱测量的优点,又能有效消除混叠现象,提高有效光程差和光谱分辨率,再对和差处理后的干涉信号进行傅里叶变换即可得到被测偏振光谱。文章对双强度...
  • 阐述了偏振光谱成像技术的基本原理,同时利用光学元件构建了静态偏振光谱成像仪系统模型,推导了入射光斯托克斯矢量光谱调制干涉光强的数据公式和该矢量的解调公式,并对整个过程进行了数值仿真分析,理论分析及模拟...
  • 基于光纤色散方法,对比分析了基于色散展宽和基于时间透镜的光学傅里叶变换的实现原理。对两类方案中5种不同系统的输出结果、色散量、光谱分辨率、时间带宽积等影响因素进行分析,分别阐述了其应用及改进方案。
  • 干涉型计算层析成像光谱仪是一种将空间调制傅里叶变换成像光谱仪(FTIS)的原理与计算层析成像光谱仪(CTIS)的原理相结合的一种新型成像光谱仪,具有高通量、高光谱分辨力以及高空间分辨力的特点。分析和讨论了干涉型...
  • 介绍空间调制干涉成像光谱技术,提出基于实体迈克尔逊干涉仪的空间调制干涉成像光谱仪,推导了这种仪器的基本几何参数公式与主要误差容限计算公式,总结了其主要特点。文中还归纳出高通量大视场干涉仪原理模型
  • 由于其原理特殊,定标方法尚不成熟。介绍了一种用于空间调制干涉光谱成像仪光谱辐射度定标的方法,即干涉光谱成像仪和光谱辐射度计同时采集目标辐射强度,复原采集干涉图得到的光谱图,与光谱辐射度计采集绝对光谱进行...
  • 根据光谱光学相干层析(OCT)技术的原理,建立了光谱光学相干层析系统,并利用该系统测量了膜系的反射光谱。该方法通过傅里叶变换实现了光学相干层析信号光谱信息的提取,克服了传统反射测量方法测量结果不直观、受非测量...
  • 高分辨率反射式转镜干涉光谱仪是一种新型傅里叶变换光谱仪,它利用倾斜旋转镜的转动产生非线性的光程差(OPD)。光程差的研究对于干涉光谱仪的性能指标分析、设计和研制具有重要意义。在介绍反射式转镜干涉光谱仪的...
  • 阐述了傅立叶变换红外光谱仪的原理,并对红外光谱 的特点进行了归纳和总结。
  • 基于傅里叶变换合成法的基本原理, 合成了一个K9基底上的负滤光片, 合成的渐变折射率薄膜具有期望光学特性, 但实际制备难度很大, 因此将其细分为足够多层离散折射率的均匀薄膜, 由于实际薄膜材料种类有限, 不能获得...
  • 从双光束干涉的基本原理出发,分析了高斯分布光源条件下利用傅里叶变换解调光纤法布里珀罗传感器的原理。针对高斯分布光源特点,提出直接对波长均匀采样得到的光谱数据进行快速傅里叶变换,然后利用传感器腔长与频域...
  • 就我而言,我其实不需要深入了解傅里叶变换的原理,但又不得不必须学会如何使用FFT对采样的数据进行处理并得到最终的图像,比如频谱图,比如色散图。所以,对于傅里叶变换的作用的直观理解,应该是从它的结果图像中...



    青山依旧在,几度夕阳红。----杨慎《临江仙·滚滚长江东逝水》


    前言

    本文在参考网上众多有关傅里叶变换资料的基础上,以如何使用程序对采样得到的数值进行傅里叶变换为角度,总结了在程序中对采样得到的离散数值进行快速傅里叶变换(FFT)的相关流程。本文不仅总结了众所周知的对时间的傅里叶变换,也总结了对空间的傅里叶变换,还总结了如何对在时间和空间上采样得到的数值进行二维傅里叶变换。本文以平面波y=A*sin(w*t - k*x)为例子,通过三个简单的程序来演示以上三种情况,分别是:①固定空间为一个点,采样得到这个空间点随时间变化的一维数值序列,然后对这个一维序列进行傅里叶变换,即可得到频谱图(A'(f)图)。②固定时间为一个时刻,在这个时刻,采样得到随空间变化的一维数值序列,然后对这个一维序列进行傅里叶变换,即可得到波矢&振幅图(A'(k)图)。③对时间范围内和空间范围内采样得到的数值矩阵进行二维傅里叶变换,以得到色散图(A'(k,f)图)。这些内容都是适用于各个领域的通用基础技巧,本文最后还以微磁学中的自旋波为例,通过解析一个现有的频谱和色散分析程序包的源码来演示如何从微磁模拟得到的数据中提取自旋波的频谱和色散信息。本文不涉及傅里叶变换的原理和公式,只是将之视为一个信号分析工具,来演示普通人该怎么操作这个工具。

    一、傅里叶变换的作用

    1.傅里叶变换的直观理解

    如果你在网上看过以下关于傅里叶变换的资料:
    傅里叶分析之掐死教程(完整版)更新于2014.06.06
    Python 中 FFT 快速傅里叶分析
    形象理解二维傅里叶变换
    【numpy】几种fft函数的使用
    那么相信你已经对傅里叶变换有了基本的认识,如果没有看过也没关系!就我而言,我其实不需要深入了解傅里叶变换的原理,但又不得不必须学会如何使用FFT对采样的数据进行处理并得到最终的图像,比如频谱图,比如色散图。所以,对于傅里叶变换的作用的直观理解,应该是从它的结果图像中来理解。比如在频谱图中,可以看到原始信号中每个频率分量的占比情况,根据频率带宽的分布可以得到系统允许频率和禁止频率。
    维基百科上的这个图像很形象表达了傅里叶变换的作用:时域上复杂的原始信号可以分解成在时域上不同频率和振幅的简单的正弦波,而把分解得到每一个正弦波的频率的和对应的振幅绘制出来,即可得到傅里叶变换的结果即频谱图 A ~ ( f ) \tilde{A}(f) A~(f)
    在这里插入图片描述
    对于获取原始信号中的频率分量的过程,其实我觉得还有一个非常直观例子,即三棱镜对白光的色散:
    在这里插入图片描述

    当一束白光(输入的原始信号)经过三棱镜后,三棱镜会将白光展开为不同的颜色(输出的各个频率的波)。那么傅里叶变换的作用也就好比三棱镜,当原始信号进行傅里叶变换后,就可以获取组成原始信号的频率贡献情况。

    2.傅里叶变换家族

    针对原始信号在时域上是否连续,是否周期的情况,傅里叶变换家族出现了四种获取相应频谱的方法:
    在这里插入图片描述

    (a)连续傅里叶变换(FT):在时域上连续的,非周期的波形,经过FT之后,得到的频谱是连续的,非周期的。
    (b)连续傅里叶级数(FS):在时域上连续的,周期的波形,经过FS之后,得到的频谱是离散的,非周期的。
    (c)离散傅里叶变换(DTFT):在时域上离散的,非周期的波形,经过DTFT之后,得到的频谱是连续的,周期的。
    (d)离散傅里叶级数(DFS):在时域上离散的,周期的波形,经过DFS之后,得到的频谱是离散的,周期的。

    就实际情况而言,我们对原始信号经过等时间间隔采样并保存得到的都是在时域上离散的数值,而且计算机对这些输入的离散数值进行傅里叶变换之后的得到的结果也是且只能是离散的数值。因此,看起来以上四种方法中适合计算机使用的只有DFS,只有该方法的输入输出都是离散数值,但由于它对输入的离散数值具有周期性的要求,所以又有了它的改进方法即DFT,这种方法对输入的离散数值没有特别的限制,大多数情况所说的傅里叶变换就是指DFT。DFT会对输入的离散数值进行如下运算:
    假设对一时变信号进行采样后得到N个数值: x ( 0 ) , x ( 1 ) , x ( 2 ) . . . x ( N − 1 ) x(0),x(1),x(2)...x(N-1) x(0)x(1)x(2)...x(N1),那把这N个数值作为DFT的输入,对于N个输入的数值来说:
    f ( k ) = ∑ n = 0 N − 1 x ( n ) ⋅ e − j ⋅ 2 π ⋅ n N ⋅ k , k = 0 , 1 , 2... N − 1 f(k)= \displaystyle\sum_{n=0}^{N-1}{x(n)} \sdot e^{-j \sdot \frac{2\pi \sdot n}{N}\sdot k},k=0,1,2...N-1 f(k)=n=0N1x(n)ejN2πnk,k=012...N1
    这个公式表示:①假设DFT有N个输入数值的话,那么经过傅里叶变换后也得到了N个输出点,因此频谱图中就有N个点。

    ②式中的 e − j ⋅ 2 π ⋅ n N ⋅ k e^{-j \sdot \frac{2\pi \sdot n}{N}\sdot k} ejN2πnk是一个复数,因此DFT的结果也是一个复数。但是对于第一个输出点 f ( k = 0 ) f(k=0) f(k=0)来说,它的DFT的结果值为 f ( 0 ) = ∑ n = 0 N − 1 x ( n ) ⋅ e − j ⋅ 2 π ⋅ n N ⋅ 0 = x ( 0 ) + x ( 1 ) + x ( 2 ) + . . . + x ( N − 1 ) f(0)= \displaystyle\sum_{n=0}^{N-1}{x(n)} \sdot e^{-j \sdot \frac{2\pi \sdot n}{N} \sdot 0}=x(0)+x(1)+x(2)+...+x(N-1) f(0)=n=0N1x(n)ejN2πn0=x(0)+x(1)+x(2)+...+x(N1) ,即第一个点的DFT结果是一个实数值。

    ③对于输出的每一个点来说,DFT都需要进行N次运算,对于总共有N个输入的采样值来说,那么DFT的算法运算复杂度为 O ( N 2 ) O(N^2) O(N2)。可见输入的点数越多,运算次数也越多。

    囿于原始DFT算法的运算速度,于是一种用于加速DFT计算的算法横空出世,那就是快速傅里叶变换(FFT),相比DFT,它的运算复杂度为 N l o g N NlogN NlogN,当输入的采样点数越多,FFT算法的速度相比原始的DFT快得多。所以在程序的数学库中用于一维傅里叶变换的函数为fft(),该函数便实现了上式的运算流程。

    二、时间、空间、时空的傅里叶变换

    1.时间的傅里叶变换

    首先需要明白一点,用于当做时间傅里叶变换的输入的采样值是怎么获得的?在现实世界中,或者通过示波器的探头接触某个器件的输出位置点进行等时间间隔采样,以此获得N个随时间变化的采样值;或者通过放在空间中某一位置点的传感器等时间间隔采样,也能获取N个随时间变化的采样值。。。它们都有一个特点:把空间固定为一点,只等时间间隔采样在该空间位置点上的原始输出信号

    假设原始信号的最高频率分量为f(非角频率w),时间采样频率为Fs_time,采样点数为N。
    以上三个条件之间有如下关系:①根据奈奎斯特采样定理,采样频率至少是原始信号的最高频率分量的两倍,即 F s _ t i m e ≥ 2 ∗ f Fs\_time≥2*f Fs_time2f。②采样频率的倒数即采样周期Ts,采样周期乘以采样点数就是采样时间。
    接下来,将从时间的FFT得到的频谱图为切入点,来讲解如何获取频谱图中每一个点,并绘制出频谱图的基本流程。
    首先得到频谱图中每个点的横坐标:
    对N个采样数值进行FFT后,也有N个复数结果值: f ( k ) , k = 0 , 1 , 2.. N − 1 f(k),k=0,1,2..N-1 f(k)k=012..N1,它们对应N个频率点,这些频率点对应的频率值为 F s _ t i m e N ⋅ k \frac{Fs\_time}{N} \sdot k NFs_timek ,其中 F s _ t i m e N \frac{Fs\_time}{N} NFs_time称为频率分辨率,即频谱图( A ~ ( f ) \tilde{A}(f) A~(f))中的两个点横坐标之间的最小间隔,因此要得到更小的频率分辨率(也就是让频谱曲线看起来更加连续),就只能增加采样点数(但不能降低采样频率),也就是延长采样时间。
    比如当频率点k=0时,其代表的频率值即为 f = F s _ t i m e N ⋅ 0 f=\frac{Fs\_time}{N} \sdot 0 f=NFs_time0 = 0Hz,当频率点k=1时,其代表的频率值即为 f = F s _ t i m e N ⋅ 1 f=\frac{Fs\_time}{N} \sdot 1 f=NFs_time1 = F s _ t i m e N \frac{Fs\_time}{N} NFs_timeHz。。。于是便得到了频谱图中所有点的横坐标。显然,如此计算的话,频谱图中的点的最大横坐标为Fs_time,也就是说,频谱图的横轴范围是受时间采样频率Fs_time影响的。
    接着是得到频谱图中每个点的纵坐标:
    对于FFT得到的除了第一个点(0Hz)以外的所有频率点,其傅里叶振幅 A ~ \tilde{A} A~为:首先对该点的复数结果取模,接着再除以采样点数N,最后再乘以2即可。对于第一个点(即直流分量0Hz),它的傅里叶振幅为取模之后,只除以采样点数N即可。如此,便得到了频谱图( A ~ ( f ) \tilde{A}(f) A~(f))的N个点的纵坐标。
    看起来有了频谱图中每个点的横坐标和纵坐标,按理说就可以绘制出一个完整的频谱图了。但是,实际上,对于FFT之后的前面一半的频率点k=0,,, N 2 − 1 \frac{N}{2}-1 2N1,它们对应的频率确实为0Hz,,, F s _ t i m e N ⋅ ( N 2 − 1 ) \frac{Fs\_time}{N} \sdot (\frac{N}{2}-1) NFs_time(2N1)Hz,但是对于后一半部分的频率点k= N 2 \frac{N}{2} 2N,,,N-1,它们对应的频率为 − N 2 ⋅ F s _ t i m e N -\frac{N}{2} \sdot \frac{Fs\_time}{N} 2NNFs_time,,, − F s _ t i m e N -\frac{Fs\_time}{N} NFs_time。显然对于时间的傅里叶变换得到的频谱图来说,负频率没有物理意义,所以需要舍弃负频率的点,使频谱图的宽度减少一半。

    什么,你居然问我怎么晓得如此操作的,有什么依据吗?
    在这里插入图片描述

    其实,网上关于傅里叶变换的资料众多,质量也是参差不齐,让人找不到北,而恰好MATLAB官方的fft()函数的文档中有关于时间的傅里叶变换的例程,见MATLAB官方参考链接,这个示例演示了如何获取频谱。而更巧的是Python的Numpy库的fft()例程居然没有获取信号频谱的例子,那我就搬一下砖,直接将MATLAB官方的傅里叶变换的例程移植到Python程序中,如下:

    #################
    #author:YQYUN
    #date:22/6/18
    #desc:从MATLAB官方例程移植而来,演示时间的傅里叶变换
    #desc:固定空间为一个点,得到在这个位置点上采样得到时间范围内的信号
    #################
    
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.pylab import mpl
    
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
    mpl.rcParams['axes.unicode_minus'] = False  # 显示负号
    
    #时间采样频率
    Fs_time =1000
    #采样周期
    T = 1/Fs_time
    #采样点数
    M = 1500
    #采样时间序列
    t = np.arange(0, M-1) * T
    
    # 设置需要采样的信号
    y = 0.7 * np.sin(2 * np.pi * 50 * t)+np.sin(2 * np.pi * 120 * t)
    #加入噪声信号
    y = y + 2 * np.random.randn(np.size(t))
    
    #画原始波形
    plt.subplot(121)
    plt.plot(1000*t[0:50], y[0:50])
    plt.xlabel('t(s)')
    plt.ylabel('A')
    plt.title('原始波形')
    
    #快速傅里叶变换得到一个复数序列,取复数的绝对值,即复数的模(双边频谱)
    abs_y = np.abs(np.fft.fft(y))
    #除了第一个点(0Hz)外,对模值要先除以采样点数量,再乘以2
    fft_y = abs_y / M * 2
    #对于第一个点(0Hz),只需对模值除以采样点数量即可 
    fft_y[0] /= 2
    #幅度取一半区间(单边频谱)
    fft_y = fft_y[0:int(M / 2)]
    #频率取一半区间
    F = np.arange(0,int(M / 2)) * Fs_time / M
    
    #将FFT结果保存到文件
    FFTResult_time = np.zeros((int(M / 2),2))
    FFTResult_time[:,0] = F[:]
    FFTResult_time[:,1] = fft_y[:]
    np.savetxt('FFTResult_time.txt', FFTResult_time,delimiter=' ')
    print('时间FFT结果已经写入FFTResult_time.txt,第一列为频率,第二列为振幅')
    
    #画图
    plt.subplot(122)
    plt.plot(F, fft_y)
    plt.xlabel('f(Hz)')
    plt.ylabel('A')
    plt.title('单边频谱')
    plt.show()
    

    首先看一下运行结果,对比MATLAB(上)和Python(下)的结果,如下:
    在这里插入图片描述

    由于加入了随机的噪声信号,所以两者的时域波形是不一样的,比如再次运行一下这个程序:
    在这里插入图片描述
    以上三个图像的时域波形看起来都是乱七八糟的,没有什么规律可言,但是从它们的频谱图中可以很明显的看出50和120Hz频率分量的对整个频谱图的贡献。傅里叶变换为我们提供了另一种视角来观测世界。

    首先是构造对原始信号进行等时间间隔采样获得的数值,此处以振幅A=0.7,频率f=50Hz,和A=1,f=120Hz的正弦波为基础波形,然后加入一个随机噪声,从而生成一个随着采样时间变化的函数,在规定的采样频率和采样点数M=1500个的条件下即可通过这个函数生成采样值。总共需要的采样时间为采样周期T乘以采样点数M。此处要注意进行傅里叶变换的前提条件:采样频率Fs_time设置为信号需要识别的最高频率分量的两倍以上,此处为Fs_time=1000Hz(原始波形中除噪声外的最高频率f=120Hz)。

    #时间采样频率
    Fs_time =1000
    #采样周期
    T = 1/Fs_time
    #采样点数
    M = 1500
    #采样时间序列
    t = np.arange(0, M-1) * T
    # 设置需要采样的信号
    y = 0.7 * np.sin(2 * np.pi * 50 * t)+np.sin(2 * np.pi * 120 * t)
    #加入噪声信号
    y = y + 2 * np.random.randn(np.size(t))
    

    将M个采样值作为函数fft()的输入参数,并将输出的复数结果取模得到abs_y ,然后按照上文所述的那样,对输出结果的第一个值(0Hz)和其他值分别处理,得到每个频率点的傅里叶振幅fft_y,也就是频谱图中每一个点的纵坐标。接着是将频率分辨率分别乘以频率点以获取每个频率点对应的频率值,也就是频谱图中每一个点的横坐标。最后舍弃掉频谱图的后一半部分的数据,即负频率的部分,只保留fft_y和F序列的前一半区间 。

    #快速傅里叶变换得到一个复数序列,取复数的绝对值,即复数的模(双边频谱)
    abs_y = np.abs(np.fft.fft(y))
    #除了第一个点(0Hz)外,对模值要先除以采样点数量,再乘以2
    fft_y = abs_y / M * 2
    #对于第一个点(0Hz),只需对模值除以采样点数量即可 
    fft_y[0] /= 2
    #幅度取一半区间(单边频谱)
    fft_y = fft_y[0:int(M / 2)]
    #频率取一半区间
    F = np.arange(0,int(M / 2)) * Fs_time / M
    
    

    如此便得到了频谱图中每个点的横坐标和纵坐标,那么可以直接使用绘图库的函数绘制这些点。若要保存频谱图的话,则可以将每一个点的横坐标和纵坐标装进一个 M 2 \frac{M}{2} 2M行,2列的矩阵中,其中第一列放点的横坐标,第二列放点的纵坐标。最后调用savetxt()函数保存该矩阵即可。

    #将FFT结果保存到文件
    FFTResult_time = np.zeros((int(M / 2),2))
    FFTResult_time[:,0] = F[:]
    FFTResult_time[:,1] = fft_y[:]
    np.savetxt('FFTResult_time.txt', FFTResult_time,delimiter=' ')
    print('时间FFT结果已经写入FFTResult_time.txt,第一列为频率,第二列为振幅')
    
    

    对于不了解MATLAB的人来说,可能对MATLAB官方程序的这行代码感到奇怪:

    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    
    

    这是由于在MATLAB中,数组的下标是从1开始的,而其他语言都是从0开始的,在移植MATLAB代码时需要注意。

    2.空间的傅里叶变换

    网上关于时间的傅里叶变换的资料是最多的,耐心查阅的话,不难搞明白它的使用技巧。接下来将介绍空间的傅里叶变换,以平面波为例,平面波的形式为 y = A ⋅ s i n ( w ⋅ t − k ⋅ x ) y=A \sdot sin(w \sdot t-k \sdot x) y=Asin(wtkx),其中k是波矢,稍微展开一下:
    y = A ⋅ s i n ( 2 π ⋅ f ⋅ t + 2 π ⋅ ( − 1 λ ) ⋅ x ) y=A \sdot sin(2\pi \sdot f \sdot t+2\pi \sdot (-\frac{1}{\lambda}) \sdot x) y=Asin(2πft+2π(λ1)x)
    其中的时间项为 2 π ⋅ f ⋅ t 2\pi \sdot f \sdot t 2πft,空间项为 2 π ⋅ ( − 1 λ ) ⋅ x 2\pi \sdot (-\frac{1}{\lambda} )\sdot x 2π(λ1)x,表示函数值y是随着时间和空间变化的。

    在上文中,把空间固定为一点,等时间间隔采样得到该空间位置点上的原始信号在某一时间范围内的数值,该一维序列经过FFT后,得到了( A ~ ( f ) \tilde{A}(f) A~(f))图。时间和空间是等同地位的,通过类比,固定时间为一个时刻,在这个时刻,等空间间隔采样得到随空间变化的一维数值序列,然后对这个一维序列进行傅里叶变换,即可得到( A ~ ( − 1 λ ) \tilde{A}(-\frac{1}{\lambda}) A~(λ1))的图像,不过通常还需要将横坐标 ( − 1 λ ) (-\frac{1}{\lambda}) (λ1)转化为波矢k,即对横坐标进行运算: ( − 1 λ ) ⋅ ( − 1 ) ⋅ ( 2 π ) = k (-\frac{1}{\lambda}) \sdot (-1) \sdot (2\pi)=k (λ1)(1)(2π)=k ,于是便得到了波矢&振幅图( A ~ ( k ) \tilde{A}(k) A~(k)
    在时间的傅里叶变换后,舍弃掉了负频率的数据(负频率无意义),但是波矢呢?波矢k是一个矢量,在此处举的例子中,k是在限制在了x轴方向上变化,但它的值有正有负,表示波的传播方向,所以就不能舍弃任何空间傅里叶变换后得到的点。
    类比上文时间的傅里叶变换,空间的傅里叶变换示例如下:

    #################
    #author:YQYUN
    #date:22/6/19
    #desc:以平面波为例,演示空间的傅里叶变换
    #desc:固定时间为一个点,在这个时刻上采样得到空间范围内的信号
    #################
    
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.pylab import mpl
    
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
    mpl.rcParams['axes.unicode_minus'] = False  # 显示负号
    
    #平面波方程:A*sin(w*t - k*x)
    
    #####构造信号#####
    
    #信号振幅:24
    A = 24
    #信号频率:8Hz
    f = 8
    #信号角频率:2*pi*f
    w = 2*np.pi*f
    #信号波长:2*pi*(1/3)
    lamb = 2 * np.pi *(1 / 3)
    #信号波矢:2*pi/lamb=3m^-1
    k = 2 * np.pi / lamb
    
    #对于固定时间点t=0s,对指定x范围内进行FFT可以得到A'(k)图
    t_fixed = 0
    
    #x空间范围:0到98m
    x_start,x_end = 0,98
    
    #空间采样率设为信号波长倒数的3倍
    Fs_volume = 3 * (1 / lamb)
    #采样空间序列
    x = np.arange(x_start,x_end,1/Fs_volume)
    #空间采样点数
    N = np.size(x)
    
    #举例子时注意信号的波矢k'不能大于k
    y1 = (A * np.sin(w * t_fixed - k * x) +
        2 * A * np.sin(w * t_fixed - 0.9*k * x) +
        3 * A * np.sin(w * t_fixed - 0.8*k * x) +
        4 * A * np.sin(w * t_fixed - 0.7*k * x))
    #画图
    plt.subplot(1,2,1)
    plt.xlabel('x (m)')
    plt.ylabel('A')
    plt.plot(x, y1)
    plt.title('t=0,原始波形')
    
    #空间的傅里叶变换
    fft_y1 = np.fft.fft(y1)
    fft_y1 = abs(fft_y1) / N * 2
    fft_y1[0] /= 2
    
    #将采样点等间隔分为序列号:0,1,2...N-1:np.arange(-int(N / 2),int(N / 2))
    #最终波矢点为k_min...0...k_max共N个
    
    #这里调用一种现成的函数,fftfreq的第一个参数为采样点数,第二个参数为采样频率
    #波长的倒数:波数,对比时间项wt=2*pi*f和空间项-kx=2*pi*(-1/lamb),可知得到的波数要乘以-1
    wave_num = -1 * np.fft.fftfreq(N,1/Fs_volume)
    #波数*2*pi等于波矢
    kx = wave_num * (2 * np.pi)
    
    
    #保存数据
    #重新排序数据,让横轴按照从小到大的顺序存储
    kx_sorted = np.sort(kx)
    temp_y = np.zeros((N,1))
    for i in range(0,N):
        temp_kx = kx_sorted[i]
        #kx里面没有重复元素,故可以这么查找
        kx_index = np.where(kx== temp_kx)
        temp_y[i,0] = fft_y1[kx_index]
    kx[0:N] = kx_sorted[0:N]
    fft_y1[0:N] = temp_y[0:N,0]
    
    #将FFT结果保存到文件
    FFTResult_volume = np.zeros((N,2))
    FFTResult_volume[:,0] = kx[:]
    FFTResult_volume[:,1] = fft_y1[:]
    np.savetxt('FFTResult_volume.txt', FFTResult_volume,delimiter=' ')
    print('空间FFT结果已经写入FFTResult_volume.txt,第一列为波矢(2*pi/lambda),第二列为振幅')
    
    #画图
    plt.subplot(1,2,2)
    plt.xlabel('kx(m^-1)')
    plt.ylabel('A')
    plt.plot(kx[0:N], fft_y1[0:N])
    plt.title('t=0,K空间')
    plt.show()
    
    

    示例中构造了(波矢,振幅)=(1k,1A),(0.9k,2A),(0.8k,3A),(0.7k,4A)共四组平面波的叠加,但是显然由于固定了时间t为一个数,所以平面波就变成了简单的sin形式。举例子时要注意加入的平面波的波矢k’不能大于待分析波矢的最大值k。

    在此处,把时间t设置为0,将波矢的最大值k设置为3/m。将采样的空间范围x_start,x_end 设为0,98m,类比时间傅里叶变换中的时间采样频率Fs_time,空间的采样频率Fs_volume也设为至少是信号波长倒数的2倍,示例中为3倍,如此便可以得到空间的采样点数N。

    #信号波长:2*pi*(1/3)
    lamb = 2 * np.pi *(1 / 3)
    #信号波矢:2*pi/lamb=3m^-1
    k = 2 * np.pi / lamb
    #对于固定时间点t=0s,对指定x范围内进行FFT可以得到A'(k)图
    t_fixed = 0
    #x空间范围:0到98m
    x_start,x_end = 0,98
    #空间采样率设为信号波长倒数的3倍
    Fs_volume = 3 * (1 / lamb)
    #采样空间序列
    x = np.arange(x_start,x_end,1/Fs_volume)
    #空间采样点数
    N = np.size(x)
    

    序列y1中就保存着N个等空间间隔采样得到的数值。类比平面波方程中的时间项 2 π ⋅ f ⋅ t 2\pi \sdot f \sdot t 2πft,时间的傅里叶变换后得到傅里叶振幅 A ~ \tilde{A} A~与频率 f f f 的关系,那么空间项 2 π ⋅ ( − 1 λ ) ⋅ x 2\pi \sdot (-\frac{1}{\lambda} )\sdot x 2π(λ1)x,相应的,空间的傅里叶变换后会得到傅里叶振幅 A ~ \tilde{A} A~ − 1 λ -\frac{1}{\lambda} λ1 的关系。只需稍加处理即可得到傅里叶振幅 A ~ \tilde{A} A~与波矢 k k k 的关系。

    #将采样点等间隔分为序列号:0,1,2...N-1:np.arange(-int(N / 2),int(N / 2))
    #最终波矢点为k_min...0...k_max共N个
    
    #这里调用一种现成的函数,fftfreq的第一个参数为采样点数,第二个参数为采样频率
    #波长的倒数:波数,对比时间项wt=2*pi*f和空间项-kx=2*pi*(-1/lamb),可知得到的波数要乘以-1
    wave_num = -1 * np.fft.fftfreq(N,1/Fs_volume)
    #波数*2*pi等于波矢
    kx = wave_num * (2 * np.pi)
    

    注意:在计算波矢点对应的波矢值时,不同于上文根据频率点自己手动计算频率值,这里使用了一个专门用于根据频率(波矢)点计算频率(波矢)值的函数fftfreq(),只需简单的为该函数指定采样点数(第一个参数)和采样频率的倒数(第二个参数)两个参数即可,比如:
    在这里插入图片描述
    可以看到fftfreq()函数的输出序列为从0到正的最大值,然后是从负的最小值到负的最大值。那么该函数的输出序列就可以直接当做傅里叶变换之后每一个点的横坐标

    有了( A ~ ( k ) \tilde{A}(k) A~(k))图中的每个点的横坐标之后,对于每个点的纵坐标则和时间的傅里叶变换中的获取纵坐标的方式相同,即把空间采样得到的N个数值作为函数fft()的输入参数,并对复数输出结果进行处理就得到了傅里叶振幅 A ~ \tilde{A} A~

    #空间的傅里叶变换
    fft_y1 = np.fft.fft(y1)
    fft_y1 = abs(fft_y1) / N * 2
    fft_y1[0] /= 2
    

    如此,一个( A ~ ( k ) \tilde{A}(k) A~(k))图中所有点的横纵坐标都得到了,就可以开始绘图步骤了。不过为了方便绘图,还需要将图中的每个点按照横坐标从负到0,再到正的顺序排序,否则绘制的曲线会出现闭合(参见fftfreq()函数的输出序列)。
    若要保存图像数据,则需要一个N行,2列的矩阵按照第一列保存每个点的横坐标,而第二列保存每个点的纵坐标。

    #保存数据
    #重新排序数据,让横轴按照从小到大的顺序存储
    kx_sorted = np.sort(kx)
    temp_y = np.zeros((N,1))
    for i in range(0,N):
        temp_kx = kx_sorted[i]
        #kx里面没有重复元素,故可以这么查找
        kx_index = np.where(kx== temp_kx)
        temp_y[i,0] = fft_y1[kx_index]
    kx[0:N] = kx_sorted[0:N]
    fft_y1[0:N] = temp_y[0:N,0]
    #将FFT结果保存到文件
    FFTResult_volume = np.zeros((N,2))
    FFTResult_volume[:,0] = kx[:]
    FFTResult_volume[:,1] = fft_y1[:]
    np.savetxt('FFTResult_volume.txt', FFTResult_volume,delimiter=' ')
    print('空间FFT结果已经写入FFTResult_volume.txt,第一列为波矢(2*pi/lambda),第二列为振幅')
    

    那么根据以上程序得到的( A ~ ( k ) \tilde{A}(k) A~(k))图如下:
    在这里插入图片描述

    3.时空的傅里叶变换—色散关系

    大家都对函数的一维图像,比如 y = A ⋅ s i n ( w ⋅ t ) y=A \sdot sin(w \sdot t) y=Asin(wt)非常熟悉,因为它只有一个自变量t和因变量y。但对于函数的二维图像,比如 y = A ⋅ s i n ( w ⋅ t − k ⋅ x ) y=A \sdot sin(w \sdot t-k \sdot x) y=Asin(wtkx)的图像可能不太熟悉,它有两个自变量:时间t和空间x,一个因变量y。要在一个平面坐标系中表示这个二维图像,可以把两个自变量分别映射到两个坐标轴上,于是时间t和空间x的不同组合值,在图像中就表现为点的位置坐标。因变量则可以用颜色来映射它的值,比如使用颜色映射"blue(#0000FF)-white(#FFFFFF)-red(#FF0000)"来表示因变量,那么blue就表示负值,white表示0,red表示正值,而具体的数值则可以按照颜色的深浅来等比例映射表示。如下图表示的是平面波 y = 24 ⋅ s i n ( 2 ⋅ π ⋅ 8 ⋅ t − 3 ⋅ x ) y=24 \sdot sin(2 \sdot \pi \sdot 8 \sdot t-3 \sdot x) y=24sin(2π8t3x) 的图像:
    在这里插入图片描述

    前文介绍的都是一维的傅里叶变换,输入的采样值都是一个1行多列的序列,输出的图像也只是一条曲线。但是色散图,表示的是傅里叶振幅 A ~ \tilde{A} A~ 与频率 f f f 和波矢 k k k 的关系(即( A ~ ( k , f ) \tilde{A}(k,f) A~(k,f))),它是一个二维的图像,而不是一维的曲线。图像中的每个点的横坐标即该点的波矢值,纵坐标即该点的频率值,而该点的颜色则表示傅里叶振幅 A ~ \tilde{A} A~

    想当然的,肯定是要同时对时间和空间进行傅里叶变换,即二维FFT,才能得到色散图。正如前文所说,时间的傅里叶变换是将一条以时间为变量的曲线,分解成一系列不同频率,不同振幅的正弦波。那么二维傅里叶变换则是将一个以时间,以空间为变量的二维图像分解成一系列不同频率,不同振幅的平面波。二维傅里叶变换的输入数据应该是在一定时间范围内(图像的纵坐标),一定空间范围内(图像的横坐标)的二维图像,对于图像中的每一个坐标点(x,y)来说(该点的值为f(x,y)),二维傅里叶变换会进行如下运算:
    F ( u , v ) = ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) ⋅ e − j ⋅ 2 π ⋅ ( u x M + v y N ) F(u,v)= \displaystyle\sum_{x=0}^{M-1} \displaystyle\sum_{y=0}^{N-1}{f(x,y)}\sdot e^{-j \sdot 2 \pi \sdot(\frac{ux}{M}+\frac{vy}{N}) } F(u,v)=x=0M1y=0N1f(x,y)ej2π(Mux+Nvy)
    可见对于一个M行N列(输入矩阵)的原始输入图像来说,二维FFT后也会得到一个M行N列(输出矩阵)的输出图像,而且输出的每一个结果F(u,v)也是一个复数。在程序的数学库中用于二维傅里叶变换的函数为fft2(),该函数便实现了上式的运算流程。
    上式的直观表现如下图,通过该式子可以得到组成二维图像的不同的平面波贡献:
    在这里插入图片描述

    平面波是指波面为平面的波,波矢表示平面波的传播方向,它与波面垂直。由于平面波相对于一维的正弦波来说,除了振幅,频率,相位外,还有一个额外的波矢,所以仍然使用一个复数来表示其振幅和相位,而二维FFT之后得到的矩阵中的每一个点,每一个点的坐标 ( u , v ) (u,v) (u,v)的则可以表示平面波的法向量,而该向量的模值即 u 2 + v 2 \sqrt{u^2+ v^2} u2+v2 ,可用于表示该平面波的频率。比如下图中的 ( u , v ) (u,v) (u,v)距离中心越远,表示平面波的频率越大,而坐标 ( u , v ) (u,v) (u,v)则表示垂直于平面波的波面的向量。
    在这里插入图片描述
    在大概了解二维傅里叶变换的作用后,接着就要考虑如何使用这个工具了。首先是确定二维傅里叶变换的输入,它是一个二维的图像(也就是输入矩阵),输入矩阵中每一个元素的坐标,就表示该点在二维傅里叶变换后得到的输出矩阵中的坐标,而输入矩阵每一个元素的值,则会作为二维傅里叶变换的输入值。以平面波 y = A ⋅ s i n ( 2 π ⋅ f ⋅ t + 2 π ⋅ ( − 1 λ ) ⋅ x ) y=A \sdot sin(2\pi \sdot f \sdot t+2\pi \sdot (-\frac{1}{\lambda}) \sdot x) y=Asin(2πft+2π(λ1)x)为例,振幅A=24,空间范围x_start,x_end为0到98m,空间采样率Fs_volume= 3 ⋅ 1 λ 3\sdot \frac{1}{\lambda} 3λ1,时间范围t_start,t_end为0到19s,时间采样率Fs_time= 3 ⋅ f 3 \sdot f 3f,混合了4个不同频率,不同波矢的平面波:(频率,波矢)=(f,k),(0.9f,0.9k),(0.8f,0.8k),(0.7f,-0.7k)。

    
    #平面波:A*sin(w*t - k*x)
    
    #####构造信号#####
    
    #信号振幅:24
    A = 24
    #信号频率:8Hz
    f = 8
    #信号角频率:2*pi*f
    w = 2 * np.pi * f
    #信号波长:2*pi*(1/3)
    lamb = 2 * np.pi *(1 / 3)
    #信号波矢:2*pi/lamb=3m^-1
    k = 2 * np.pi / lamb
    
    #x空间范围:0到98m
    x_start,x_end = 0,98
    
    #空间采样率设为信号波长倒数(波数)的3倍
    Fs_volume = 3 * (1 / lamb)
    #采样空间序列
    x = np.arange(x_start,x_end,1/Fs_volume)
    
    #t时间范围:0到19s
    t_start,t_end = 0,19
    
    #时间采样率设为信号频率的3倍
    Fs_time = 3 * f
    #采样时间序列
    t = np.arange(t_start,t_end,1/Fs_time)
    
    #空间采样点数
    N = np.size(x)
    
    #时间采样点数
    M = np.size(t)
    
    dataMatrix = np.zeros((M,N))
    #注意举例子用的示例信号的w'和k'不要大于最大的频率w和波矢k
    for index_row in range(0,M):
        for index_column in range(0,N):
            dataMatrix[index_row,index_column] = (A * np.sin(w * t[index_row] - k * x[index_column]) +
            A * np.sin(0.9*w * t[index_row] - 0.9*k * x[index_column]) +
            A * np.sin(0.8*w * t[index_row] - 0.8*k * x[index_column]) +
            A * np.sin(0.7*w * t[index_row] + 0.7*k * x[index_column]))
            
    #绘图
    
    #plt.rcParams['figure.figsize'] = (12,8)
    #横轴为空间,纵轴为时间
    extent_x_t = [np.min(x),np.max(x),np.min(t),np.max(t)]
    plt.subplot(2,1,1)
    plt.xlabel('x (m)')
    plt.ylabel('t(s)')
    plt.imshow(dataMatrix,  cmap='bwr',extent=extent_x_t)
    #颜色代表的标量值
    plt.colorbar(label='y=A*sin(w*t - k*x)')
    
    plt.title('原始二维波形,A=24,f1=8,k1=3,f2=0.9f1,k2=0.9k1,f3=0.8f1,k3=0.8k1,f4=0.7f1,k4=-0.7k1')
    
    

    如此构造的二维图像如下,它的横轴表示空间,纵轴表示时间:
    在这里插入图片描述

    通过时间范围和时间采样率确定了在时间上采样得到M个点,通过空间范围和空间采样率确定了空间上采样得到的N个点,将采样得到的数值放在一个M行N列的矩阵中,当进行二维傅里叶变换后,那么矩阵的每一行的就表示一个频率点,而每一列就表示一个波矢点,而由该行列确定的位置点就存放该点进行二维傅里叶变换的结果值。使用二维傅里叶变换fft2()函数之后,要先将零频中心移动到输出矩阵的中心,即直接调用fftshift()函数即可,如此,输出矩阵的纵轴从上到下依次是正频率,0Hz,负频率;输出矩阵的横轴从左往右依次是负波矢,0,正波矢。之后同样的要去除M行N列输出矩阵中后一半行数的数据,即舍弃负频率的数据。
    有了频率点和波矢点,只需按照前文的方法对这些频率点和波矢点进行一些运算即可获取对应的频率值和波矢值。即( A ~ ( k , f ) \tilde{A}(k,f) A~(k,f))图中每一个点的横纵坐标就算出来了。
    接着是计算( A ~ ( k , f ) \tilde{A}(k,f) A~(k,f))图中每一个点的值,即傅里叶振幅 A ~ \tilde{A} A~。M行N列的矩阵进行二维FFT之后,输出结果也是M行N列的矩阵,经过移频并舍弃负频率的数据之后,就变成了 M 2 \frac{M}{2} 2M行N列的矩阵,而且里面的每个元素都是复数。想当然的,必须先对该复数进行取模使之变成实数值。而下一步对这个 M 2 \frac{M}{2} 2M行N列的实数矩阵的处理就比较五花八门了,有对这个实数矩阵进行归一化处理的,也有进行取对数操作的。在下面给出的示例中就只进行取模,没有进一步处理。

    #二维FFT
    fft_Matrix = np.fft.fft2(dataMatrix)  
    #将kx=0,f=0移动到输出矩阵的中心,并求复数的模
    fft_Matrix = abs(np.fft.fftshift(fft_Matrix))
    #对于频率部分,只取正频率,要舍弃矩阵一半的行数据
    fft_Matrix = fft_Matrix[0:int(M/2),:]
    
    ####计算坐标轴####
    
    #实空间转K空间,波长的倒数:波数,对比时间项wt和空间项-kx,可知得到的波数要乘以-1
    wave_num = -1 * np.fft.fftfreq(N,1/Fs_volume)
    #波数*2*pi等于波矢
    kx = wave_num * (2 * np.pi)
    
    #时间转频率
    F = np.fft.fftfreq(M,1/Fs_time)
    F = F[0:int(M/2)]
    
    #保存数据
    np.savetxt('FFTResult_2D_Matrix.txt', fft_Matrix,delimiter=' ')
    print('二维FFT结果已写入FFTResult_2D_Matrix.txt,矩阵第一个值对应图像第一个像素点(左上角为坐标原点)的值,大小为:',fft_Matrix.shape)
    np.savetxt('kx.txt', kx,delimiter=' ')
    print('图像的横坐标即波矢(2*pi/lambda)已写入kx.txt')
    np.savetxt('F.txt', F,delimiter=' ')
    print('图像的纵坐标即频率已写入F.txt')
    
    #画图
    #横轴为波矢,纵轴为频率(只取正频率)
    extent_kx_F = [np.min(kx),np.max(kx),np.min(F),np.max(F)]
    plt.subplot(2,1,2)
    plt.xlabel('kx (m^-1)')
    plt.ylabel('F(Hz)')
    #参数origin : {‘upper’, ‘lower’}将数组的[0,0]索引放在轴的左上角或左下角,‘upper’通常用于矩阵和图像
    plt.imshow(fft_Matrix,cmap='bwr',origin='upper',extent=extent_kx_F)
    #颜色代表的标量值
    plt.colorbar(label='|fft2(y)|')
    plt.title('色散关系')
    plt.show()
    

    示例中使用了savetxt()将色散图的横坐标和纵坐标分别保存,然后将最终得到的实数矩阵保存下来。最后则是用颜色绘制函数绘制出矩阵中的每个位置的值映射到色条(colorbar)中的颜色,这样就得到了 M 2 X N \frac{M}{2}XN 2MXN个像素点的色散图。比如矩阵中第一行第一列的元素即第一个元素matrix[0,0],它就应该绘制在色散图的左上角,即第一个像素。
    将以上两段代码合并后运行,得到的输出图像如下所示:
    在这里插入图片描述
    可以在色散图中找到组成原始图像的4个平面波的贡献情况:(频率,波矢)=(f,k),(0.9f,0.9k),(0.8f,0.8k),(0.7f,-0.7k)。
    以上内容便是适用于各个领域的傅里叶变换的通用基础技巧,接下来是介绍如何获取磁体系中传播的自旋波的频谱和色散关系,大家在阅读这部分的内容后,对于分析其他领域的时空波形来说可能启发一些相似的思路。

    三、自旋波的频谱,色散

    1.已有的自旋波分析程序

    微磁模拟软件OOMMF和Mumax3都是使用OVF格式的文件来保存磁体系的磁化状态的,在进行频谱和色散分析之前,务必保证每个磁化文件之间的是按照相同间隔时间(即采样周期)来保存的。 比如模拟运行5ns,按照每5ps保存一次磁体系的磁化状态,那么模拟结束后会得到1000个磁化文件。对应到前文傅里叶变换中的参数,也就是时间采样周期Ts为5ps(时间采样频率Fs_time=1/Ts),时间采样点数M为1000。而空间采样周期(空间采样频率Fs_volume的倒数)呢?由于磁体系的空间分为x,y,z三个方向,所以就有三个方向上的空间采样周期,分别等于用户在模拟时在三个方向上设置的单元格尺寸。
    用于分析自旋波频谱和色散关系的程序包有很多,比如在笔记03中提到过的MFA,又比如在笔记04中提到过的semargl3,其中MFA程序包是以Python源码的方式给出的,它要求磁化文件的数据格式是text文本格式,用户需要在代码里设置需要分析的参数,然后直接运行即可,由于有非常清晰的说明文档,所以MFA使用起来是非常简单的。而semargl3则是一个可执行程序,它要求磁化文件的数据格式是b8或者b4,而且要求每个磁化文件必须以.ovf作为后缀名,用户需要在semargl3的图形界面设置参数。它对自旋波的分析类型是我见过的程序包中最多的,不过它的说明文档不太详细,所以使用起来有一定的门槛。
    此外,硕士论文《在磁畴壁中传播的自旋波束缚态的探测原理》的附录中给出了使用MATLAB获取自旋波的色散的代码,在OOMMF官方教程的最后几页PPT里面也有使用Python获取自旋波的色散的代码。
    那么此处以MFA的源码为例,将解析它的内容(实际上就是加一点注释方便理解而已)来看一下怎么处理磁化文件来获取自旋波的相关信息。

    2.频谱

    (1)首先需要从所有磁化文件中导入三个磁化分量mx,my,mz中的一个磁化分量,该磁化分量是用于表征自旋波的振幅的。MFA采用多进程同时读取所有的磁化文件的某一磁化分量,返回得到了一个长度为文件数量的列表list(array[…],…),列表的每一个元素是一个一维数组,数组中保存着该磁体系所有单元格的某一磁化分量:

    #读取一个磁化文件中的某一磁化分量,name:文件名称;colum:读取的矢量分量
    def fileload(name,colum):
        colum = int(colum)
        #磁化文件中'#'是注释
        m_oneFile = np.loadtxt(name, dtype = float,
                          comments = '#', usecols = (colum,))
        return m_oneFile
    
    #该函数的返回值为[array([...]), array([...])...],列表的一个元素表示保存一个磁化文件的某个分量的浮点型数组
    #列表长度为当前文件夹中所有子文件数量
    #dirc:存放磁化文件的目录
    #start_stage,end_stage:需要的磁化文件的范围
    #exportX,exportY,exportZ:选择磁化矢量的那个分量作为自旋波的振幅
    #num:占用计算机核心的数量
    def readUseMultiCore(dirc,start_stage,end_stage,exportX,exportY,exportZ,num):
        #磁化文件中的第三列表示Z分量
        if exportZ == 1:
            column = '2'
        #磁化文件中的第二列表示Y分量
        elif exportY == 1:
            column = '1'
        #磁化文件中的第一列表示X分量
        else:
            column = '0'
    
        #定义一个保存磁化文件名称的空列表
        files = []
        #walk函数会以文件夹为单位自顶而下遍历指定文件目录,
        #返回的三个参数分别表示当前文件夹, 当前文件夹中的子文件夹, 当前文件夹中的子文件
        for i,j,k in os.walk(dirc):
            #获取当前文件夹中所有磁化文件,并保存到列表中
            #因此,这里的files是二维列表:files[k[...]]
            files.append(k)
        #将files转化为一维列表
        files = files[0]
        #对列表中的磁化文件进行按名称升序排序
        files.sort()
        files_sort = files[:]
        #创建一个保存磁化文件全路径名称的列表
        files_new = []
        for i in range(start_stage,end_stage):
            files_new.append(dirc+ '/' + files_sort[i])   
        if num == 0:
            #没有指定核心数量时,表示使用所有cpu核心
            pool = Pool()
            #将需要扩展的的原始函数作为partial函数的第一个参数,原始函数的一系列参数作为partial的后续参数,最后生成一个带默认参数的新函数
            #接着pool.map会把新函数(第一个参数)和新函数的一系列参数(第二个参数即列表)映射到多进程池中运行
            #此处表示,返回读取所有磁化文件中的指定磁化分量的值
            m_allFiles = pool.map(partial(fileload,colum=column),files_new)
            #关闭进程池
            pool.close()
            #主进程阻塞等待子进程的退出,join方法要在close或terminate之后使用
            pool.join()
        else:
            pool = Pool(processes = num)
            m_allFiles = pool.map(partial(fileload,colum=column),files_new)
            pool.close()
            pool.join()
        return m_allFiles
    
    

    (2)接着是对每一个磁化文件中的指定区域进行处理,从磁体系所有单元格中找到用户指定区域范围内的单元格,并对这些单元格重新排序使之成为一个一维的单元格链。
    在这里插入图片描述

    比如第一个磁化文件中的指定空间范围:x_start到x_end,y_start到y_end,z_start到z_end,那么根据三个方向上的单元格尺寸就可以得到指定区域左下角的第一个单元格的索引index,然后按照先x方向,后y方向,最后z方向找到指定区域的单元格,并把这些单元格中重新存储到一个一维链中得到cell_1,cell_2…。如此操作结束后,将这对第一个磁化文件处理后得到的单元格链放入新矩阵的第一行,然后对剩余的磁化文件也是如此处理,将得到的每一个单元格链依照保存的时间顺序放在新矩阵的每一行。最终得到的新矩阵的形状是M行N列的,行表示时间采样点,列表示空间采样点。

    #将正在遍历的磁化文件中的磁化矢量的某个分量装载到特定格式的矩阵中
    #x_num,y_num,z_num:指定区域包含的三个坐标方向上的单元格数量
    #x_cellSize,y_cellSize,z_cellSize:xyz三个坐标方向上的单元格的尺寸
    #x_start,y_start,z_start:用户指定区域的开始坐标
    #x_totalNum,y_totalNum:整个磁体系的x和y方向上的单元格数量
    #M_matrix:行数为文件数量,列数为指定区域包含的总单元格数量,元素的值为零
    #file:保存着一个矢量场文件的某个分量的数组
    #i:正在遍历的矢量场文件的索引
    def loadDataToMatrix(x_num,y_num,z_num,x_cellSize,y_cellSize,z_cellSize,
               x_start,y_start,z_start,x_totalNum,y_totalNum,M_matrix,file,i):
        for z_index in range(0,z_num):
            for y_index in range(0,y_num):
                for x_index in range(0,x_num):
                    #待填充的矩阵的索引,按照一维展开的方式来排序
                    indexOfMatrix = int(z_index * (x_num * y_num) + y_index * x_num + x_index)
                    #指定区域的单元格的索引,从左下角开始
                    indexOfSelectedRegion = (int(x_start / x_cellSize) + x_index) + ((int(y_start / y_cellSize) + y_index) * x_totalNum) + (
                    (int(z_start / z_cellSize) + z_index) * (x_totalNum * y_totalNum))
                    #填充时间点(行)单元格(列)矩阵
                    M_matrix[i,indexOfMatrix] = file[indexOfSelectedRegion]
        return M_matrix
    
    def GetDocs(m_allFiles,start_stage,end_stage,x_totalSize,y_totalSize,z_totalSize,x_start,y_start,z_start,
                x_end,y_end,z_end,x_cellSize,y_cellSize,z_cellSize,
                exportX,exportY,exportZ):
        #整个磁体系的x和y方向上的单元格数量
        x_totalNum,y_totalNum = int(np.round(x_totalSize / x_cellSize)), int(np.round(y_totalSize / y_cellSize))
        #指定区域包含的三个坐标方向上的单元格数量
        x_num = int(np.round((x_end - x_start) / x_cellSize))
        y_num = int(np.round((y_end - y_start) / y_cellSize))
        z_num = int(np.round((z_end - z_start) / z_cellSize))
        #指定的文件范围和指定区域包含的总单元格数量
        selectedFlieNum = end_stage - start_stage + 1
        selectedCellNum = x_num * y_num * z_num
        if exportX == 1:
            #np.zeros(m,n)返回一个m行,n列,用0填充的矩阵
            Mx_matrix = np.zeros((selectedFlieNum,selectedCellNum))
        else:
            Mx_matrix = 0
        
        if exportY == 1:
            My_matrix = np.zeros((selectedFlieNum,selectedCellNum))
        else:
            My_matrix = 0
        
        if exportZ == 1:
            Mz_matrix = np.zeros((selectedFlieNum,selectedCellNum))
        else:
            Mz_matrix = 0
        
        #读取到的磁化文件的数量
        flieNum = len(m_allFiles)
        # 读取指定区域单元格的磁化分量
        if exportZ == 1:
            for i in range(0,flieNum):
                file = m_allFiles[i]
                Mz_matrix = loadDataToMatrix(x_num,y_num,z_num,x_cellSize,y_cellSize,z_cellSize,
                x_start,y_start,z_start,x_totalNum,y_totalNum,Mz_matrix,file,i)
        else:
            Mz_matrix = Mz_matrix
        
        if exportY == 1:
            for i in range(0,flieNum):
                file = m_allFiles[i]
                My_matrix = loadDataToMatrix(x_num,y_num,z_num,x_cellSize,y_cellSize,z_cellSize,
                x_start,y_start,z_start,x_totalNum,y_totalNum,My_matrix,file,i)
        else:
            My_matrix = My_matrix
        
        if exportX == 1:
            for i in range(0,flieNum):
                file = m_allFiles[i]
                Mx_matrix = loadDataToMatrix(x_num,y_num,z_num,x_cellSize,y_cellSize,z_cellSize,
                x_start,y_start,z_start,x_totalNum,y_totalNum,Mx_matrix,file,i)
        else:
            Mx_matrix = Mx_matrix
        
        print('所有磁化文件的磁化分量值已读取到行(时间)列(单元格)的矩阵中')
        return Mx_matrix, My_matrix, Mz_matrix
    

    这一部分代码的难点在于用户要对OVF格式的矢量文件要有一个深入的了解,默认情况下空间坐标原点在磁体系的左下角,指定区域的范围是:x_start,y_start,z_start,x_end,y_end,z_end,那么该区域在三个方向上的单元格数量分别为:x_num,y_num,z_num,指定区域的左下角的第一个单元格索引为x_start除以x_cellSize,对于指定区域每一行的单元格,索引值就用一个范围为x_num的for循环进行遍历即可得到该行内的单元格,而对于指定区域的y方向的单元格索引,相应的用一个范围为y_num的for循环进行遍历即可,且在y方向每增加一次,则应该跳过磁体系在x方向的所有单元格,z方向同理,z方向每增加一次,则应该跳过磁体系在xy平面的所有单元格。

    (3)计算频谱:通过前文的描述,我们需要固定空间为一个点(此处最小的空间点就是一个单元格),得到该点关于时间变化的磁化分量。但是此处用户的指定区域有N个单元格,即N个空间采样点,该如何压缩使之成为一个点呢?
    可能有人马上想到了,在同一个时间采样点上,直接把N个点的磁化分量加起来,然后除以N,即取用户指定区域所有单元格的平均磁化分量,如此即可把M行N列矩阵压缩为M行1列的形状了。然后直接对这个一维序列进行一维傅里叶变换,得到的频谱即为最终结果。
    可能还有人想到另一种方法,即先对指定区域的每一个单元格进行一维傅里叶变换(即矩阵的列方向为傅里叶变换的方向),结果仍是一个M行N列的矩阵,之后在对指定区域的所有单元格取傅里叶振幅的平均值,将这个平均频谱作为最终结果。
    以上方法都是可行的!参考文献《Micromagnetic Simulations in Magnonics》(DOI:10.1007/978-3-642-30247-3_8),文中描述了三种处理多个单元格频谱的方法:
    在这里插入图片描述
    (8.1)表示将从所有磁化文件中减去基态磁化状态,并除以饱和磁化强度来实现将要分析的磁化分量进行归一化。
    (8.3a)对应在每一个时间点上,先把空间上所有单元格的磁化分量平均为一个单元格的磁化分量,然后对这个单元格的所有时间点上的磁化分量进行一维傅里叶变换,然后就可以得到频谱图了。
    (8.3b)表示先对每一个单元格在所有时间点上进行一维傅里叶变换,接着对所有单元格的傅里叶变换结果取平均值,如此也得到了一个频谱图。
    (8.3c)则表示对在所有时间点上对所有单元格进行二维傅里叶变换,然后可以选择特定波矢区域的单元格的频谱图。

    在这篇参考文献中,比较了这三种方法的结果差异如下:

    在这里插入图片描述
    对于(8.3a),由于是先在空间上取平均,它会降低空间中非均匀模式的自旋波的频率贡献。从图中可以看到,相比(8.3b)和(8.3c),它的频谱曲线没有太多的峰,也没有丰富的频谱细节。
    (8.3b)和(8.3c)则是先对空间中的所有点求出频谱,如此便不会降低空间中非均匀模式的自旋波的频率贡献,因此可以得到任何模式的自旋波的的频谱,但是正如参考文献中所说的,由于数值模拟的误差和采样原因,可能会出现一些人为的假谱线。

    从MFA的源码中可以看到,MFA采用的方法是(8.3b),即先对指定区域中的所有单元格求出频谱,然后再取该区域中的平均频谱。

    #对矩阵(行代表时间,列代表空间单元格)进行一维傅里叶变换
    #Mx,My,Mz:三者之一保存着待分析的矩阵
    #export_x,export_y,export_z:待分析的磁化分量
    #T:采样周期
    #f_min,f_max:设置最终频谱图的频率上下限
    def matrixFFT(Mx,My,Mz,export_x,export_y,export_z,T,f_min,f_max):
        if export_z == 1:
            M_matrix = Mz
        elif export_y == 1:
            M_matrix = My
        else:
            M_matrix = Mx
        
        #得到M_matrix的行(M)和列(N)数,
        #行是采样时间间隔:t0,t1,t2...tM-1,列是指定区域单元格数量cell0,cell1...cellN-1
        M,N = M_matrix.shape
        #时间采样频率
        Fs_time = 1 / T
        #将采样点等间隔分为序列号:0,1,2...M-1
        n = np.linspace(0,M-1,M)
        #根据每个序列号获取对应的频率值(即横轴坐标),频率分辨率为采样频率/采样点数
        fn = Fs_time / M * n
        #N个实数点进行FFT之后也有N个复数点,这里需要对每个单元格进行FFT,所以结果点数为N*M
        fftMatrix = np.zeros((N,M))
        #在时间(M行)-单元格(N列)矩阵中,按照列的方向来求每个单元格的FFT
        for i in range(0,N):
            #将FFT得到的M*N个数据点按照N行(表示单元格),M列(表示时间)来存放在fftMatrix中
            fftMatrix[i,:] = np.abs(np.fft.fft(M_matrix[:,i]))
            #对FFT后的复数结果的模值进行处理,得到正确的幅值:
            #除了第一个点(0Hz)外,对模值要先除以采样点数量,再乘以2
            fftMatrix[i,:] = fftMatrix[i,:] / M * 2
            #对于第一个点(0Hz),只需对模值除以采样点数量即可
            fftMatrix[i,0] = fftMatrix[i,0] / 2
        
        #对于每个频率点来说,对应的幅值是对所有单元格在该频率点下的幅值取平均值
        aveMatrix = np.zeros((1,M))
        #axis为0(1)表示对矩阵的每一列(行)取均值
        #对每一列的所有单元格求平均值,返回 1*M列的矩阵
        aveMatrix[0,:] = np.mean(fftMatrix[:,0:M], axis = 0,keepdims = True)
        
        #只取前一半部分正的频率点,舍弃后一半部分负的频率点
        F = fn[range(int(M / 2))]
        aveMatrix = aveMatrix[0,range(int(M / 2))]
        
        #归一化幅度
        aveMatrix = aveMatrix / np.max(aveMatrix)
        #将结果保存到文件
        FFTResult_time = np.zeros((len(aveMatrix),2))
        #第一列为频率
        FFTResult_time[:,0] = F
        #第二列为振幅
        FFTResult_time[:,1] = aveMatrix
        #列之间用空格分隔
        np.savetxt('FFTResult_time.txt', FFTResult_time,delimiter=' ')
        print('FFT结果已经写入FFTResult_time.txt,第一列为频率,第二列为振幅')
        
        #####绘图#####
        #变换横坐标(HZ->GHz)
        F = F / 1e9
        #矩阵.T返回矩阵的转置矩阵
        plt.plot(F,aveMatrix.T,linewidth = 2)
        plt.xticks(fontsize=18)
        plt.yticks(fontsize=18)
        plt.xlabel(r'$\mathrm{Frequency\enspace (GHz)}$',fontsize=18)
        plt.ylabel(r'$\mathrm{Normalized\enspace FFT\enspace power\enspace (arb.units)}$',fontsize=18)
        #频谱图中显示的频率上下限
        plt.xlim((f_min,f_max))
        plt.show()
        return;
    

    该部分代码没有任何难点,有了前文的铺垫,是很容易理解的。

    3.色散

    在计算色散时,需要的是x,y,z三个空间方向中的一个。在计算自旋波的频谱时,我们得到了用于一维FFT的M行N列的输入矩阵,其中N是三个方向总的单元格数量,所以需要将它压缩到一个方向。MFA源码中为x方向,具体处理方式为:将三个方向总的单元格数量(x_num X y_num X z_num)压缩到待分析的波矢的方向的单元格数量(x_num),即对y方向和z方向的单元格中的磁化分量取平均值。如此,M行N列矩阵变成了M行x_num列的矩阵。
    对这个M行x_num列的矩阵进行二维傅里叶变换,之后也得到了一个M行x_num列的输出矩阵,而且输出矩阵中的每个元素都是复数。接着要将零频(f=0,k=0)移动到输出矩阵的中心,即直接调用fftshift()函数,然后和前文一样,舍弃输出矩阵中一半的行数据(负频率),于是输出矩阵变成了 M 2 \frac{M}{2} 2M行x_num列了。之后还要对这个 M 2 \frac{M}{2} 2M行x_num列的矩阵取模,归一化,取对数等操作。

    #对矩阵(行代表时间,列代表空间单元格)进行二维傅里叶变换
    #Mx,My,Mz:三者之一保存着待分析的矩阵
    #export_x,export_y,export_z:选择作为自旋波的振幅的磁化分量
    #Ts,Vs:时间采样周期和空间采样周期(即波矢方向的单元格尺寸)
    #T_resample,X_resample:是否对时间和空间重采样
    #win,Resample_Switch:是否加窗,重采样
    #f_cutoff,klim_cutoff:截止频率和截止波矢
    #x_start,x_end:x方向的开始和结束坐标,确定了波矢为x方向
    #x_cellSize:单元格x方向的尺寸
    
    def matrixFFT2(Mx,My,Mz,export_x,export_y,export_z,Ts,Vs,
              T_resample,X_resample,win,Resample_Switch,
              f_cutoff,klim_cutoff,x_start,x_end,x_cellSize):
        if export_z == 1:
            M0_matrix = Mz
        elif export_y == 1:
            M0_matrix = My
        else:
            M0_matrix = Mx
        #M为时间,n为指定区域的总单元格(x*y*z方向)
        M, n = M0_matrix.shape
        #指定区域的x方向包含多少个单元格
        N = int((x_end - x_start) / x_cellSize)
        M_matrix = np.zeros((M,N))
        #只需要x方向(波矢方向)的单元格,即对波矢的垂直方向的单元格求平均
        start = 0
        end = N
        count = 0
        while end <= n:
            M_matrix = M_matrix + M0_matrix[:,start:end]
            start += N
            end += N
            count += 1
        M_matrix = M_matrix / count
    
        ###重采样和加窗###
        #对原始矩阵进行重采样后,矩阵的形状会改变
        M_fft, N_fft = M, N
        #定义一个新矩阵来保存对原始矩阵进行重采样后的结果
        Matrix = np.zeros((M_fft,N_fft))
        
        #重采样
        if Resample_Switch == 1:
            #对原始矩阵进行重采样后,矩阵的形状会改变
            M_fft, N_fft = round(M / T_resample), round(N / X_resample)
            j = 0
            for i in range(0,N-1):
                if np.mod(i,X_resample) == 0:
                    Matrix[:,j] = M_matrix[:,i]
                    j += 1
        else:
            Matrix = M_matrix
        
        #加窗
        if win == 1:
            from scipy.signal import chebwin as chebwin
            attenuation = 50.0
            win1 = chebwin(M_fft,attenuation)
            win2 = chebwin(N_fft,attenuation)
            win1 = win1[np.newaxis,:]
            win2 = win2[np.newaxis,:]
            win1 = win1.T
            win = np.dot(win1, win2)
            Matrix = Matrix * win
        else:
            win = win
        
        #二维傅里叶变换
        fftMatrix = np.zeros((M_fft,N_fft))
        fftMatrix = np.fft.fft2(Matrix)
        #移频之后取模值
        fftMatrix = np.abs(np.fft.fftshift(fftMatrix))
    
        #(1 / Ts)为时间采样频率,频率分为两半,并把单位转化为GHz
        F = (1 / Ts) / 2 / 1e9
        #(1 / Vs)为空间采样频率,波矢分为两半,并把单位转化为nm^-1(计算之后为负值)
        kx = (1 / Vs) / 2 / 1e9 * (-2 * np.pi)
        
        #fftMatrix的行是时间->频率(只选取正频率),列是空间->波矢(选取所有点)
        fftMatrix = fftMatrix[0:round(M_fft/2), 0:N_fft-1]
        #对FFT结果模值的归一化取分贝值
        fftMatrix = 10 * np.log10(fftMatrix / np.max(fftMatrix))
    
        #根据矩阵中所有元素获得一个平均值
        M_mean = np.mean(fftMatrix)
        #获取矩阵中最大值
        M_max = np.max(fftMatrix)
        
        #clip(a, a_min, a_max, out=None, **kwargs):将矩阵a中的值限制在最小值和最大值之内
        fftMatrix = np.clip(fftMatrix,M_mean,M_max)
        
        #画图
        #图像横向方向的放大因子
        multi_horizontal = 2e2
        X_neglim, X_poslim = kx * multi_horizontal, -kx * multi_horizontal    #在横轴方向放大图像
        Y_neglim, Y_poslim = 0, F
        extent = [X_neglim,X_poslim,Y_neglim,Y_poslim]
        plt.figure()
        plt.rcParams['figure.figsize'] = (9.0,8.5)
        #参数:origin : {‘upper’, ‘lower’}将数组的[0,0]索引放在轴的左上角或左下角。‘upper’通常用于矩阵和图像。
        #cmap = plt.cm.jet将标量数据映射到色彩图:蓝-青-黄-红
        #extent:(left, right, bottom, top)数据坐标中左下角和右上角的位置。 如果为“无”,则定位图像使得像素中心落在基于零的(行,列)索引上。
        plt.imshow(fftMatrix, cmap = plt.cm.jet, origin='upper',extent=extent)
        
        #在横轴方向放大图像
        klim_cutoff = klim_cutoff * multi_horizontal
        if klim_cutoff in [20, 40, 60, 80, 100, 120, 140, 160, 180, 200]:
            plt.xticks([-200,-180, -160,-140,-120,-100,-80,-60,-40,-20,
                0,20,40,60,80,100,120,140,160,180,200],
                ['-1.0','-0.9','-0.8','-0.7','-0.6','-0.5','-0.4','-0.3','-0.2','-0.1',
                '0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1.0'],fontsize = 18)
        elif klim_cutoff in [2, 6, 10, 14, 18]:
            plt.xticks([-18,-14,-10,-6,-2,0,2,6,10,14,18],
                ['-0.09','-0.07', '-0.05','0.03','0.01', 
                '0','0.01', '0.03', '0.05','0.07','0.09'],fontsize = 18)
        else:
            plt.xticks([-16,-12,-8,-4,0,4,8,12,16],
               ['-0.08', '-0.06', '-0.04','-0.02','0','0.02', '0.04', '0.06','0.08'],fontsize = 18)
        plt.yticks(fontsize = 18)
        plt.xlim(-klim_cutoff,klim_cutoff)
        plt.ylim(0,f_cutoff)
        #将色条减少为图形高度的一半
        plt.colorbar(shrink = 0.5)
        plt.xlabel(r"$\mathrm{Wave\enspace vector}\enspace k_x\enspace \mathrm{(nm^{-1})}$",
                   fontsize = 17)
        plt.ylabel(r'$\mathrm{Frequency\enspace (GHz)}$',fontsize = 17)
        plt.savefig('Dispersion curve.eps', dpi=500)
        plt.show()
        return fftMatrix
    

    最后绘制这个 M 2 \frac{M}{2} 2M行x_num列的矩阵的图像时进行了额外的处理,即关于坐标轴的刻度标注问题:
    ①纵轴直接按照从小到大的顺序标注即可,不过需要注意,纵轴的起始点为0Hz。
    ②横轴则是进行放大处理,此处就是将色散图在横向方向上放大了2e2倍。而且相比于MFA原来只显示正波矢的的限制,此处则取消了这个限制,运行MFA自带的示例,对比取消波矢限制的前(左图)后(右图)如下:
    在这里插入图片描述


    后记

    终于对傅里叶变换有了一个较为深入的认识了。还记得本科时,《信号与系统》和《数字通信原理》两个课程简直是工科生质检员,而我一直到本科毕业也没有搞懂傅里叶变换,现在想来也是一阵唏嘘啊。文中举例的用的数字来自一位朋友的生日:九八年,八月初三。最后希望当大家在生活中遇到看起来无比复杂的困难时,都能像傅里叶变换一样,找到困难的本质并解决它。
    展开全文
  • 分析了双通道剪切干涉以及基于微偏振阵列调制的傅里叶变换光谱偏振成像原理, 论述了光谱信息反演方法以及偏振信息提取方法。搭建了实验装置, 对实际场景目标进行光谱偏振成像实验, 获得了目标的高空间分辨率光谱图像...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 670
精华内容 268
关键字:

傅里叶光谱原理