联系我们

OASP 2 - 傅里叶变换

发表时间: 2023-09-23 14:15:59 作者: 6LNZ系列组合米机

  傅里叶变换是很有用的工具;不仅在时序上能确定周期/频率,之后还会看到一个光谱仪也相当于一个傅里叶变换器。在这里有一点需要明确的是我们以后会经常讨论光谱,所以傅里叶变换的

  对于傅里叶变换的直观理解知乎上已经有了很多不错的答案,我觉得只要了解可以将任何函数拆成各种各样的三角函数相加就可以了。F,f一般都是复函数,只有当F(x)为实偶函数的时候对应的f(σ)才是实函数。

  如果讨论谱线,此时傅里叶变换的意义就是谱线和xx轴围成的面积,也就是总吸收量。

  这里给出了常见的傅里叶变换以及图像。fft的结果因为是取模后画的图所以相当于解析的结果取了绝对值。

  单个\delta函数的傅里叶变换的振幅是常数,相位随着\sigma线性增加;两个正而且对称的\delta函数的傅里叶变换是余弦函数,一正一负(对称)的傅里叶变换的振幅是正弦函数。

  实际上观测到的数据都是离散的点,所以我们要考虑这样的离散化会带来什么样的结果。

  仍然从连续函数出发,离散化意味着将连续函数与Shah函数相乘(以波长为横坐标):

  实际测量的函数等于Shah函数与理想的连续函数相乘。同时一般来说,无论是时变数据还是光谱数据,我们都不可能无限地测量下去(望远镜时间不是无限的,光谱仪有截止波长),所以我们还需要加上一个响应度函数B(\lambda),来限制波长范围以及修改每个波长上的响应度(时间域上的数据应该不需要响应度):

  这个时候如果想求D(\lambda)的傅里叶变换的话就必须了解到几个函数相乘的傅里叶变换是什么,就涉及到卷积了。

  将其中一个函数以0为中心翻转、从左往右平移,两函数每点相乘在求和的结果就是这两个函数在第一个函数对称中心处的值。

  实际上两函数相乘的傅里叶变换就是这两函数各自傅里叶变换的卷积(倒过来也对),即

  上有证明;这里为了理解方便将证明过程倒过来写一遍。简单起见我们证明两函数各自傅里叶变换的卷积的逆傅里叶变换为两函数的乘积,可用富比尼定理,则:

  同理,洛仑兹与洛仑兹函数的卷积也还是个洛仑兹函数,只不过现在\beta_c = \beta_a + \beta_b。

  高斯函数和洛仑兹函数的卷积叫做福格特(Voigt)函数,它的傅里叶变换为

  我们现在可以回到D(\lambda)中。除了截断和采样之外,观测仪器当然也会在最终的结果中留下一些特征,叫做仪器轮廓(instrumental profile)。一般来说这些特征是一个\delta函数;也就是说一个很窄的\delta函数在经过观测仪器之后会变成另一个宽一点的\delta函数,或者说两个\delta函数卷积了。仪器轮廓可以被测出来(通过观测很细的谱线章),我们在这里记作I(\lambda)。所以现在的D(\lambda)变成了

  可以看到最终傅里叶域的函数被乘上了一个i(\sigma)。一般仪器轮廓是一个比较窄的高斯或者洛仑兹轮廓,所以它的傅里叶变换是一个相对宽的轮廓,对\sigma域中的“低频”部分没有特别大的影响;但是轮廓会在高\sigma部分下降,所以会抹掉高\sigma部分的信息。在\lambda域上这在某种程度上预示着每个数据点不再是一个孤立的点,而因为与仪器轮廓的卷积会和旁边的数据点耦合,将高\lambda的信息“平均”掉。

  采样对高\sigma信号的也会产生一定的影响(可能比仪器轮廓更大)。信号与Shah函数的卷积III(\sigma) * f(\sigma)会使得信号在Shah函数每个不为0处附近重复,使得一样的图样发生重叠,同样抹掉了高x的信息。因为Shah函数傅里叶变换后成为间隔为\frac{1}{\Delta x}的Shah函数,所以我们大家都认为在这个间隔的一半区间内,信号是没有重叠的,叫做奈奎斯特频率:

  在这里因我们对信号进行了采样,所以得到的其实就是离散的信号点;我们并不知道点和点之间究竟是如何的(虽然经常认为是平滑过渡的),而点和点之间丢失的信息就对应着频域中大于奈奎斯特频率的部分,它们重叠了。与仪器轮廓的不同在于采样是离散化,而仪器轮廓是模糊。

  我们想做的是经过测量D(\lambda)求得d(\sigma),再转换成f(\sigma)最后得到F(\lambda)。从d(\sigma)到f(\sigma)这一步有III(\sigma)的作用,它们都抹掉了高\sigma的信息。如果f(\sigma)在高\sigma处趋于0,那么转换这一步绝大多数都是正确的。但是如果重叠的部分很宽,这部分就不会有唯一的解,处理方法是减少采样的间隔。实际上f(\sigma)一般的确在高\sigma处趋于0,而且i(\sigma)的存在也会强制使得f(\sigma)在高\sigma的值减小,所以转换是成立的。

  最后一个需要提及的是b(\sigma)。只要B(x)够宽,问题就不大;但是不够宽的时候b(\sigma)本身的宽度和在\sigma \frac{1}{W}处的起伏都可能会导致数据的模糊,所以有必要修改B(x)的边缘以消除起伏。b(\sigma)本身的宽度会使和f(\sigma)的卷积很大,所以有必要稍微减小f(\sigma)。12章会有更多介绍。

  用定义计算有些时候会带来非常大的计算量;而我们大家常常遇到的数据是均匀采样的,这样一个时间段用快速傅里叶变换(FFT)可以很方便地计算。

  FFT中的输入和结果都是离散的,也就是说D(x)变成了D(j\Delta x),j是整数,从0到N。同理,将结果也离散化成\sigma = k \Delta \sigma,从傅里叶变换的定义出发我们大家可以得到

  从上一节我们大家都知道,采样后数据的傅里叶变换只在奈奎斯特频率以下包含有用数据,所以结果的计算也只到\sigma_N = \frac{1}{2\Delta x}为止。同时我们要求结果的长度和输入的数组长度一致,而且中间为0,所以频域的步长

  因为采样的缘故,fft只计算到奈奎斯特频率的位置。而如果原函数不是实数函数,那么傅里叶变换关于0点对称的复数值将不会共轭,所以numpy中的fft先计算0到\sigma_N的值,再计算0到-\sigma_N的值,并使傅里叶变换的长度与源函数长度一致(一半分给正频率,一半给负频率)。这给我们画图带来了一些麻烦,所以numpy提供了numpy.fft.fftfreq(n[, d])来给出傅里叶变换的x值,只需输入原函数长度和步长即可。上面的画图中因为结果已知而且基本上特征都在0附近,所以人为地截断到了和原函数同样的x值上;也能看出fft的结果分辨率比解析的差一点,是采样带来的后果。

  有关模,fft结果乘上步长之后取模才是正确的大小,是因为fft只计算了前两项,没有乘上最后的\Delta x;当然这是因为程序不知道步长的缘故了。

  代进定理8导出来的式子里面,再加上\Delta x和\Delta \sigma的关系,我们有

  也就是说想减少傅里叶变换的误差,能够大大减少测量误差或者缩短采样间隔,或者增大采样范围。实际测量的时候只把关心的数据点纳入(比如说某条谱线),将其他的点改成平均值或者1这种误差为0的数值能够大大减少测量误差。

  当误差与\lambda,\sigma无关,并且采样间隔小到除了\sigma域的主体之外,还有一部分高\sigma的变换也在\sigma_N之内时,我们通过可以高\sigma部分直接测量S_\sigma从而算出S_x。

  对于有周期的信号来说,长时间观测之后可以求出周期;但是观测时间不长或者周期性不明显的时候,可优先考虑使用傅里叶变换。但是要注意由于采样引起的假信号(Shah函数),不一定每个峰都是真的。