精华内容
下载资源
问答
  • 我们利用正交偏振方法将傅里叶谱的实部和虚部分开,从而测定了一维和二维傅里叶谱的相位。文中报道了位移狭缝和位移方孔傅里叶谱相位的测量结果,它们与用位移量计算出的相位值符合很好,因此证实了该方法的可靠性。在...
  • 傅里叶谱和地震动中的傅里叶振幅是有些许差别的。主要表现在幅值的调整上。 参考资料:《地震动的分析入门》大崎顺彦 傅里叶谱 一般的傅里叶谱的幅值利用fft变换后乘以2再除以N进行调整。 import numpy as...

    引言

    傅里叶谱和地震动中的傅里叶振幅谱是有些许差别的。主要表现在幅值的调整上。

    参考资料:《地震动的谱分析入门》大崎顺彦

    傅里叶谱

    一般的傅里叶谱的幅值利用fft变换后乘以2再除以N进行调整。

    import numpy as np
    import math
    import matplotlib.pyplot as plt
    
    sampRat = 100
    T = 6
    
    t = np.linspace(0, T, T*sampRat, endpoint=False)
    y = 5*np.sin(2*math.pi*t*5)+10*np.sin(2*math.pi*t*25)
    
    plt.figure()
    plt.title('时程曲线')
    plt.plot(t, y)
    plt.xlabel('时间')
    plt.ylabel('幅度')
    plt.show()
    
    f = np.linspace(0, sampRat, T*sampRat, endpoint=False)
    ff = np.fft.fft(y)
    ff = np.abs(ff)
    ff = ff*2/sampRat/T
    
    plt.figure()
    plt.title('傅里叶谱')
    plt.plot(f, ff)
    plt.show()


    地震动傅里叶振幅谱

    地震动傅里叶振幅谱利用fft变换后除以采样频率即可。

    # -*- coding: utf-8 -*-
    
    import numpy as np
    import math
    import matplotlib.pyplot as plt
    
    sampRat = 100
    T = 6
    
    t = np.linspace(0, T, T*sampRat, endpoint=False)
    y = 5*np.sin(2*math.pi*t*5)+10*np.sin(2*math.pi*t*25)
    
    plt.figure()
    plt.title('时程曲线')
    plt.plot(t, y)
    plt.xlabel('时间')
    plt.ylabel('幅度')
    plt.show()
    
    f = np.linspace(0, sampRat, T*sampRat, endpoint=False)
    ff = np.fft.fft(y)
    ff = np.abs(ff)
    ff = ff/sampRat
    
    plt.figure()
    plt.title('傅里叶振幅谱')
    plt.plot(f, ff)
    plt.show()










    展开全文
  • 介绍一种利用运动目标图象傅里叶谱位相信息实现对目标定位的新方法,该方法具有较高定位精度,且不受运动目标姿态改变的影响。文中介绍原理、给出计算机仿真和实验测量结果,并讨论实际应用的可行性。
  • Matlab求地震波傅里叶谱
  • 傅里叶变换及切比雪夫方法(配点)法求解PDE 以下的内容介绍的是傅里叶谱方法求解PDE、切比雪夫方法求解PDE(这里指的是配点法)以及一种先进行差分离散,再对离散系统的每个变量使用离散傅里叶级数展开的求解...

    傅里叶变换及切比雪夫谱方法(配点)法求解PDE

    以下的内容介绍的是傅里叶谱方法求解PDE、切比雪夫谱方法求解PDE(这里指的是配点法)以及一种先进行差分离散,再对离散系统的每个变量使用离散傅里叶级数展开的求解PDE的方法。

    因为时间仓促,文中的很多公式,我没能自己手打出来,而用简单地截图来替代,因此也导致了本文在排版上比较混乱,图片大小不一。另外,也因为截图的原因,导致符号不一致(因为截自不同的书籍),符号系统略显混乱,希望读者在阅读时心中自明。

    本文内容如下:第一,我想谈谈和傅里叶变换有关的一些东西,以及一种方法使用离散傅里叶变换去求解pde的方法。这和傅里叶谱方法和切比雪夫谱方法是不同的,它是一种离散优先的方法。第二,就是所谓的傅里叶谱方法求解PDE,这和那个“高阶有限元方法”的那个谱方法不尽相同。第三,我想简单地谈一下chebyshev谱方法。

    差分离散的傅里叶方法求解PDE

    一般的傅里叶变换的定义是:

    在这里插入图片描述

    我们倾向于实用角速度,而不是频率,也就是说:

    在这里插入图片描述

    那么,傅里叶变换就变为了:

    在这里插入图片描述

    傅里叶变换有一些性质,比如说:

    在这里插入图片描述

    如果我用双向箭头来表示傅里叶变换,即:

    在这里插入图片描述

    那么,我们有:

    在这里插入图片描述

    更重要的是,如果我定义卷积(convolution)如下:

    在这里插入图片描述

    那么,有一个卷积定理:

    在这里插入图片描述

    并且,还有Parserval定理:

    在这里插入图片描述

    Parserval定理告诉我们,傅里叶变换在时域(频域)上是平方可积的,那么在(频域)时域上也是平方可积的。当然,傅里叶变换还有很多其他的性质,我不想把它们全都列出来。


    下面,我想介绍一个很重要的东西,即香农采样定理。有了它,离散傅里叶变换也就成了顺理成章的事情了。

    我们用Δ\Delta来表示时域上的采样时间间隔,我们能都达到一系列离散的采样点:

    在这里插入图片描述

    我们定义一个特别的频率fcf_c,被称为Nyquist critical frequency(奈奎斯特),

    在这里插入图片描述

    它只是和时间间隔Δ\Delta是有关系的。

    什么是采样定理呢?它是说:

    在这里插入图片描述

    我们用ω\omega来替代2πf2\pi f,那么写成中文的形式,有:

    在这里插入图片描述

    这个告诉我们什么呢?我们能让采样间隔足够小,以至于f^\hat f的支集包含在Nyquist 频率之间,那么采样就没有信息损失。

    这个定理非常非常地重要,就像泰勒展开在数值计算中的作用一样。

    那么,如果这个条件不被满足,那么会发生什么呢?alising,中文中有的人称为“混叠”,也有叫“截频”之类的,我不知道他们是怎么翻译的。

    如下图所示:

    在这里插入图片描述

    那个aliased的那条线其实啊,就是用原来的那个(12.1.3)得到的,只不过不是真正的傅里叶变换之后频域上的取值。其实,误差会以某种方式被推到并且累积到Nyquist区间中(其实有个两者关系的一个表达式)。在我的观点中,Nyquist采样定理是连接连续和离散的桥梁。


    下面我们来谈一谈离散傅里叶变换。

    我们取:

    在这里插入图片描述

    为了方便,我们假设N是偶数。那么,让我来寻找在傅里叶变换之后,函数值在下述这些点上的估计:

    在这里插入图片描述

    那么,如果我们的采样满足前面的采样条件,且h(t)h(t)的非零值都集中在一个有限的区间内的话,我们做一个傅里叶变换的离散,能得到以下的公式。可以看到,这个约等号可以使用采样定理h(t)h(t)的表达式得到,因为那个类似于sinc函数一样的东西,做个傅里叶变换得到的就是指数函数。从另一个角度看,我们通过梯形求积分公式,也能得到类似的表达。

    在这里插入图片描述

    那么,我们再令:

    在这里插入图片描述

    有:

    在这里插入图片描述

    (12.1.7)就是所谓的离散傅里叶变换。容易看到,离散傅里叶变换后的值是周期性的,有周期为N。为了保持一致,我们喜欢变换后的点的下标n从0到N-1。但是从(12.1.5),我们能知道,频率的值是从-N/2 到 N/2。所以呢,有:

    在这里插入图片描述

    DFT有一些性质,和连续傅里叶变换是类似的,没有必要去提它了。

    最后,我们定义反(inverse)傅里叶变换如下:

    在这里插入图片描述


    傅里叶变换的用途是什么呢?它可以帮我们求解微分方程,我等等会提到这个。通过以上的讨论,我们知道,对于周期函数,我们可以通过采样,做离散的傅里叶展开,那么对于非周期边界条件呢?

    下面,我想介绍一下,快速正弦变换,和快速余弦变换。在如下所示的情况下,我们可以用sine变换,和consine变换。

    在这里插入图片描述

    首先,我需要介绍一下什么是sine变换。

    在这里插入图片描述

    对于fjf_j,首先可以关于j=N点,做一个twice的奇延拓,以至于,我们有:

    在这里插入图片描述

    然后,我们再用傅里叶变换:

    在这里插入图片描述

    我们可以把这个式子分成两部分,后面一部分,做一个替换:

    在这里插入图片描述

    把两部分合成一部分,我们有:

    在这里插入图片描述

    当然,我们也能够通过fft算法来进行离散sine变换。

    那么离散余弦变换是什么呢?和离散正弦变换的推导相同,此时我们只是在N处做偶延拓,使得:

    在这里插入图片描述

    相同的推导过程,我们能得到:

    在这里插入图片描述

    这就是离散余弦变换。它也有一些快速的算法通过fft。离散余弦变换的另一种形式是:

    在这里插入图片描述

    这里的一撇表示当k=0的时候,那一项前面有个系数0.5。它是原有序列关于点N-1/2做偶延拓得到的。

    二维的傅里叶变换定义为:

    在这里插入图片描述

    容易看到,它只是在不同的方向上做一维的傅里叶变换的结果。


    到这里,我们已经讨论了离散傅里叶变换,离散正弦变换和离散余弦变换。下面我想说一种使用FFT求解PDE的一种方法,它不同于傅里叶谱方法,和chebyshev谱方法,它是离散优先的一种方法。快速傅里叶方法直接可用,当方程空间表达项前面的系数是常数的时候。

    我们以椭圆方程的例子举例。对于问题:

    在这里插入图片描述

    由有限差分方法得到:

    在这里插入图片描述

    重新唤起二维的离散傅里叶变换:

    在这里插入图片描述

    在这里插入图片描述

    带入到差分问题。我们能得到:

    在这里插入图片描述

    从这里,我们可以看到,我们只需要做一次傅里叶变换,做一次除法,接着傅里叶反变换,我们就完成了pde的求解。

    以上的过程对周期性边界是有效的。也就是说:

    在这里插入图片描述

    在这里插入图片描述对于齐次的dirichlet边界变换问题,我们使用sine展开:

    在这里插入图片描述

    和前面相似,过程是:

    非齐次边界条件怎么办?书上有一些技巧去处理它。这里不再提了。有兴趣的读者,可以翻翻书。

    对于齐次的纽曼边界条件来说,我们使用cos展开:

    在这里插入图片描述


    对于以下的一个问题,我们可以设计一个傅里叶变换解法。我们使用cos变换展开。

    在这里插入图片描述

    时间上,我使用隐式欧拉格式。

    离散我们的问题,可以得到:
    uj,ls+1=uj.ls+Δth2(uj1,ls+1+uj+1,ls+1+uj,l1s+1+uj,l+1s+14uj,ls+1)u_{j,l}^{s+1}=u_{j.l}^s+\frac{\Delta t}{h^2}(u_{j-1,l}^{s+1}+u_{j+1,l}^{s+1}+u_{j,l-1}^{s+1}+u_{j,l+1}^{s+1}-4u_{j,l}^{s+1})

    替换cos变换进去,cos变换的定义在这。容易知道:
    uj1,ls+1+uj+1,ls+1=2J2Lm=0Jn=0Lu^m,ns+1cos(πlnL)[cos(πjmJ)cos(πmJ)+sin(πjmJ)sin(πmJ)+cos(πjmJ)cos(πmJ)sin(πjmJ)sin(πmJ)]=2J2Lm=0Jn=0Lu^m,ns+1cos(πjmJ)cos(πlnL)2cos(πmj)u_{j-1,l}^{s+1}+u_{j+1,l}^{s+1} = \frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi ln}{L}) \\ [cos(\frac{\pi j m}{J})cos(\frac{\pi m}{J})+sin(\frac{\pi j m}{J})sin(\frac{\pi m}{J})\\ +cos(\frac{\pi j m}{J})cos(\frac{\pi m}{J})-sin(\frac{\pi j m}{J})sin(\frac{\pi m}{J})]\\ =\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})2cos(\frac{\pi m}{j})

    可以得到:
    uj1,ls+1+uj+1,ls+1+uj,l1s+1+uj,l+1s+14uj,ls+1=2J2Lm=0Jn=0Lu^m,ns+1cos(πjmJ)cos(πlnL)2[cos(πmj)+cos(πnL)2]u_{j-1,l}^{s+1}+u_{j+1,l}^{s+1}+u_{j,l-1}^{s+1}+u_{j,l+1}^{s+1}-4u_{j,l}^{s+1}\\ =\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})2[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]

    所以,我们有:

    2J2Lm=0Jn=0Lu^m,ns+1cos(πjmJ)cos(πlnL)=2J2Lm=0Jn=0Lu^m,nscos(πjmJ)cos(πlnL)+Δth2[2J2Lm=0Jn=0Lu^m,ns+1cos(πjmJ)cos(πlnL)2[cos(πmj)+cos(πnL)2]]\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})\\ =\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L}) \\ +\frac{\Delta t}{h^2}[\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})2[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]]

    最后,我们得到(为什么能去掉sum符号?因为e指数表达式一组基)
    u^m,ns+1=u^m,ns+Δth2[u^m,ns+12[cos(πmj)+cos(πnL)2]]\hat u^{s+1}_{m,n} = \hat u^{s}_{m,n} +\frac{\Delta t}{h^2}[\hat u^{s+1}_{m,n}2[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]]

    也就是说,

    u^m,ns+1=u^m,ns/{12Δth2[cos(πmj)+cos(πnL)2]}\hat u^{s+1}_{m,n} = \hat u^{s}_{m,n}/\{1 - \frac{2\Delta t}{h^2}[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]\}

    思考: 这里看起来和sin变换的出来的结果是一样的,这是一个巧合吗?是的。不过,我们这里对uu做的是cos变换,故而整个结果还是不一样的。


    下面,我想通过一个一维的例子说说所谓的傅里叶谱方法求解PDE。

    在这里插入图片描述

    过程是:

    在这里插入图片描述

    在这里插入图片描述

    在这里插入图片描述

    那么,我们可以用龙格库塔格式求解了。也就是:

    在这里插入图片描述


    接下来,我想介绍得是chebyshev谱方法。它有点像差分方法,只不过它用的是切比雪夫点,而不是均匀的网格点。

    什么是切比雪夫点呢?

    在这里插入图片描述

    在这里插入图片描述

    在这里插入图片描述

    我们可以将单位圆上半圆周平均分为N份,等分节点在x轴上的投影就是切比雪夫点。

    下面我想介绍的是使用切比雪夫求导矩阵去求一些函数在切比雪夫点上的导数。

    想要做的就是:给你一个函数在切比雪夫点上的值,你需要寻找这些点上的导数值。

    我们得到离散导数值可以通过两步,即先插值,得到插值多项式,在求导,去这些点上的值:

    在这里插入图片描述

    这个算子是线性的,所以可以写为:

    在这里插入图片描述

    一个矩阵乘上一个向量。先考虑简单的情况。

    如果N=1N=1,那么插值点为x0=1,x1=1x_0=1,x_1=-1。拉格朗日插值为:

    那么求导,得到:

    在这里插入图片描述

    所以呢,求导矩阵为:
    在这里插入图片描述

    考虑N=2N=2的情况,插值节点为1 、0、 -1,拉格朗日插值为:

    在这里插入图片描述

    导数(derivative)为:

    在这里插入图片描述

    这一次,求导矩阵为:

    在这里插入图片描述

    更一般,我们有一个定理。

    在这里插入图片描述

    在这里插入图片描述

    有了谱求导矩阵,我们就可以使用如下近似(approximation):

    在这里插入图片描述


    如何使用切比雪夫求导矩阵去求解PDE呢?让我们看一些简单的例子。

    考虑一维泊松问题:

    在这里插入图片描述

    通过离散,它可以写为:

    在这里插入图片描述

    如何处理狄利克雷边界条件?这里有N-1个未知数(unknown number),但是有 N+1 个方程。条件太多啦。怎么处理?如下图:

    在这里插入图片描述

    我们先去掉求导矩阵的第一行和最后一行,以及右端项f的第一个值和最后一个值。再由狄利克雷边界条件,删掉求导矩阵的第一列和和最后一列,以及u的第一个值和最后一个值。我们用D~,u~\tilde D,\tilde u来表示处理后的结果。最后,求导结果就是:

    在这里插入图片描述

    如果区间不是从-1到1,我们只要在求导矩阵前面乘以一个缩放因子(Scale Factor)。

    另一个例子是:

    在这里插入图片描述

    u~\tilde u表示除边界外所有函数的值,即:

    在这里插入图片描述

    D~\tilde D也是一样的意思。那么我们有这样的近似:

    在这里插入图片描述

    那么纽曼边界条件怎么办?这也有一个小例子。

    在这里插入图片描述

    边界条件可以写为:

    在这里插入图片描述

    我们可以交换求导顺序(为什么可交换?):

    在这里插入图片描述

    最后,我们能得到这样一个问题:

    在这里插入图片描述

    我们可以使用Dn2D_n^2去表示2/x2\partial ^2/\partial x^2。在时间方向上,我们可以用龙格库塔格式。问题是,边界条件怎么办?

    not so easy!

    如果一个问题有纽曼边界条件,我们有:

    在这里插入图片描述

    u(±)=0u'(\pm)=0的条件下:

    在这里插入图片描述

    到这,我们将边界上的狄利克雷约束条件转化为了简单的边界值表达。

    回到原问题,我们将u/t\partial u / \partial t视作uu,那么我们得到u/t\partial u / \partial t在边界上的一个表达,直接替换这个表达到切比雪夫离散的右端项上,这样我们时间步(time-stepping)上的格式能够进行。对于切比雪夫离散矩阵,我们将边界上的值替换为其内部(interior)值得一个线性表出。

    现在,让我们把这种方法应用到我们的问题(之前写过)上,即:

    utΔu=0 in Ω=[0,1]2uηΩ=gηΩu(0,x,y)=g(x,y)u_t - \Delta u = 0 \ in \ \Omega=[0,1]^2\\ \frac{\partial u}{\partial \eta}\mid _{\partial \Omega}=\frac{g}{\eta}\mid _{\partial \Omega}\\ u(0,x,y) = g(x,y)

    容易convert这个问题为:

    utΔu=0 in Ω=[0,1]2ηut=0u(0,x,y)=g(x,y)u_t - \Delta u = 0 \ in \ \Omega=[0,1]^2\\ \frac{\partial }{\partial \eta}\frac{\partial u}{\partial t}=0\\ u(0,x,y) = g(x,y)

    然后,我们在不同的方向上,对如上的问题如法炮制,我们能够很快地求解这个问题。


    下一文章,我想介绍的是,傅里叶级数(三角形式和指数形式)、连续傅里叶变换、离散时间傅里叶变换、傅里叶变换的所需要的条件(对函数的要求),以及他们之间的一个关系。敬请期待。

    在这里插入图片描述

    展开全文
  • 傅里叶谱方法及其工程实现

    千次阅读 多人点赞 2018-11-19 20:28:11
    傅里叶谱方法求解PDE 快速傅里叶变换 用一维举例,在(−∞,+∞)(-\infty,+\infty)(−∞,+∞)有定义且绝对可积,它的傅里叶变换及其逆变换为: u^(k)=∫−∞+∞u(x)e−ikxdxu(k)=12π∫−∞+∞u^(x)eikxdx \hat u(k)=...

    傅里叶谱方法求解PDE

    快速傅里叶变换

    用一维举例,在(,+)(-\infty,+\infty)有定义且绝对可积,它的傅里叶变换及其逆变换为:
    u^(k)=+u(x)eikxdxu(k)=12π+u^(x)eikxdx \hat u(k)=\int_{-\infty}^{+\infty}u(x)e^{-ikx}dx\\ u(k)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\hat u(x)e^{ikx}dx
    傅里叶变换及其逆变换称为傅里叶变换对,傅里叶变换对不是唯一的。两个定义式子中的系数可以随意修改,只要他们的乘积是12π\frac{1}{2\pi}就行了。当然定义式中的eikxe^{-ikx}eikxe^{ikx}也是可以互换的,它们叫做记分变换的核。
    傅里叶变换,大家喜欢说是从时域(空域)到频域上的变换。看图稍微理解一下,不再细述。因为不是关键。

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    频域的基本单元也可以理解为一个始终在旋转的圆。简单地说就是,如果傅里叶变换表示的是从时域到频域上的变换,那么kk表示的就是频率(角速度,有人用ω\omega来表示),u^(x)\hat u(x)表示的就是频谱。类似,如果xx表示的是空间坐标,那么傅里叶变换就是从空域到频域上的一个变换。那么kk表示的就是波数(2π2\pi 长度上的波长的个数,代表空间上的变化速度,暗含频率的意思),此时u^(k)\hat u(k)称为波数谱。

    简单地说,傅里叶变换及其逆变换就是从时域(空域)到频域之间相互转换的工具。

    DFT和FFT

    DFT,离散傅里叶变换,就是傅里叶变换的离散化。简单地说,就是把一个N长的向量,通过乘以一个MxN的矩阵,变成M长的向量(一般可以去M=N)。重要的是,这个变换矩阵是什么呢?其实它是一个范德蒙德矩阵,所以DFT其实就可以写为:

    在这里插入图片描述

    那么这个ω0,ω1...ωM1\omega_0,\omega_1...\omega_{M-1}是个什么玩意儿呢?它其实就是复单位圆周上的M个等分点位置。如下:

    在这里插入图片描述

    一言以蔽之,要将一个N长的向量变成M长的向量,将复平面的上单元圆分成M分,取第1个等分点,乘方出来N个数,和原来的N长向量做內积,得到傅里叶变换出来的第1个数,依次类推,每个等分点都能出来一个数,那么M个等分点就能出来M个数。

    我想到这,我应该把离散傅里叶变换将清楚了。非要用公式来表达的话,可以写成这样。

    u^k=j=0N1ei2πkMjuj,j=1,2,...,N\hat u_k = \sum\limits_{j=0}^{N-1}{e^{i\frac{2\pi k}{M}\cdot j}\cdot u_j},j=1,2,...,N

    上面提到的都是一维的傅里叶变换,二维傅里叶变换其实就是在两个维度上分别做傅里叶变换。

    在实际计算中,人们不直接采用上述离散傅里叶变换的定义式进行计算,而是使用快速傅里叶变换(Fast Fourier Transform),我们可以简单地先把快速傅里叶变换认为是一种做离散傅里叶变换的快速算法。

    快速傅里叶变换被誉为二十世纪最伟大的算法之一。DFT算法的运算量为O(N2)O(N^2),而快速傅里叶变换将运算量降为了O(NlogN)O(NlogN)。但是,为了获得较高的运算速度,待变换序列的元素数量必须是2n2^n形式个。

    傅里叶谱方法

    傅里叶谱方法的思路:对PDE方程进行FFT变换,得到只对时间微分的常微分方程组。然后,使用一个基于时间步的数值格式来解ODE方程组。

    在matlab中,采用的离散傅里叶变换的定义为(下面只考虑N=M):

    在这里插入图片描述

    因为matlab的下标是从1开始,这样定义式比较好的。注意到上述定义中的归一化系数也可以其他选择,但它们的成绩必须为1N\frac{1}{N}。并且,由此定义,我们可以看到u^k=u^k+N,uj=uj+N\hat u_k=\hat u_{k+N},u_j=u_{j+N},也就是说,离散傅里叶变换,隐含着周期性边界条件。

    在matlab中,调用傅里叶变换和逆变换的函数为fft(u)\text{fft(u)}ifft(u)\text{ifft(u)}。若u为以矩阵,则fft返回矩阵每一列的离散傅里叶变换。若要计算矩阵u的每一行的离散傅里叶变换,只需调用fft(u,[],2)即可,它比先转置做傅里叶变换,再转置回来要快一些。

    需要极其注意的一点是,fft作用于离散函数向量值uu,它所对应的u^\hat u,对应的频域上的定义域(频率kk)不是从小到大排列的,而是在0位置砍断之后的一个左右平移交换。什么意思呢?下面需要仔细说说这个东西,这个对傅里叶谱方法的实现尤为重要。

    空域(时域)上的序列u1,u2,,uNu_1,u_2,\ldots ,u_N在x轴上的间隔是LN\frac{L}{N},相应的横坐标为:

    在这里插入图片描述

    为了方便,这里的区间我们取为[L/2,L/2][-L/2,L/2],因为我们分成N分,产生了N+1个点,由于周期性边界的假设,我们抛弃最右边那个点,得到了N个点。那么uu序列做傅里叶变换的输出序列u^1,u^2,,u^N\hat u_1,\hat u_2,\ldots ,\hat u_N,其对应的频率的顺序排列为:

    在这里插入图片描述

    注意到,这正是在2πN/L2\pi N/L的区间长度上的等分点,在中间割开之后,左右进行了互换。

    那么在变换过去的频域上的区间长度和两点间隔是多少呢?看下面这个表格。

    在这里插入图片描述

    怎么记?有这样一个规律,上表格的对角线相乘是横等于2π2\pi的,所以只要知道了空域坐标的分割情况,就知道频域坐标的分割。为什么这样?现在我们在乎的工程实现,其中的道理先不用深究。

    在matlab当中啊,有fftshift函数和ifftshift函数来调整输出的u^\hat u的顺序,使其与顺序的频域坐标对应上。比如fftshift(fft(u)),就使得,输出的u^\hat u对应的频域坐标是按从小到大(递增)的顺序排列的。当然,在程序中我们考虑变频域坐标是更为方便的,否则在做逆变换的时候,就需要将u^\hat u变过去,再变回来。什么意思呢?后面在程序中就能看到这一点。ifftshift函数用于进行与fftshift函数相反的操作。

    二维的傅里叶变换matlab调用语法为fft2(u),由于二维傅里叶变换是做两次一维傅里叶变换,所以其等价于fft(fft(u),[],2),即先对列做傅里叶变换,再对行做傅里叶变换。同样,二维傅里叶变换,变换过去的值,对应的频率(空间坐标)也不是对应的。那么,它也有一个移位函数fftshift(u)。易知,这个函数是将u^(x,y)\hat u(x,y)关于x轴左右互换,再关于y轴上下互换得到的,其实就是将其一、三象限交换,二、四象限交换并返回。其实就是先左右倒一下,再上下倒一下。

    也就是说,做傅里叶变换之后,低频部分到了两端,高频部分到了中间。在实际中,我们常常使用abs函数求绝对值,来得到频谱的振幅。

    下面的一个问题是,求解PDE的时候为什么要做傅里叶变换?它到底有什么用呢?

    若用F表示傅里叶变换,若满足:在x|x| \rightarrow \infty时,有u(x)0u(x)\rightarrow 0,那么由分部积分公式容易得到:

    F[u(n)(x)]=(ik)nF[u(x)]{F[u^{(n)}(x)]}=(ik)^nF[u(x)]

    这里的k是频域坐标。导数的傅里叶变换,相当于傅里叶变换穿进去,出来一个ikik的n方。也就是说,函数的求导运算在傅里叶变换的作用下,可以转化为相对简单的代数运算。也就是要求n阶到,只需要做个傅里叶变换,乘(ik)n(ik)^n,在变回去就完事了。这是其在求上的意义。实际操作当中,我们要注意到这个kkF[u(x)]F[u(x)]不是直接对应的,差了一个移位关系,所以在做乘法的时候要注意。对于高维,比如说二维,这个表达式也是类似的,我们可以直接根据定义,对两个方向分别做傅里叶变换和分布积分得到结果,也可以直接把xx视为一个向量,利用格林公式来处理,得到的是相同的结果。

    同样地,对上式进行稍加处理,得到:

    F1{F[u(n)(x)]/(ik)n}=u(x)F^{-1}\{F[u^{(n)}(x)]/(ik)^n\}=u(x)

    那么傅里叶变换就具有了求不定积分的意义。用这个来计算积分,在实际操作中,可能会出现k=0k=0,而导致式子出现无穷大,通常可以通过将kk修改为一个很小的数,或者分母加上很小的数来解决。

    傅里叶变换还可以用来求解PDE,什么意思呢?我们知道非定常的发展PDE,方程中既带有对时间的导数,又带有对空间的导数,从上面可以看到,对于空间的导数,我们关于空间变量做一个傅里叶变换,就不含对空间(这里的空间指的是频率空间)的导数了,而对时间的导数项,做傅里叶变换,不影响空间坐标上的变化。所以,通过傅里叶变换,我们可以将u,u/x,2/x2u,\partial u/\partial x ,\partial ^2 / \partial x^2 \ldots变为u^,iku^,(ik)2u^\hat u,ik\hat u,(ik)^2\hat u。容易知道,通过这样的变化,PDE就变成了一个只关于时间变量t的ODE。然后,我们再用一般的ODE的数值方法求解即可,比如用matlab内置ode45,一般使用时ode相关函数时,优先选择这个。最后,把结果,再做一个傅里叶逆变换,变回来,就完事了。是不是很简单呢?

    需要注意的一点是,傅里叶谱方法求解PDE一般默认边界条件为周期性的边界添加,也就是说一端边界处的函数值将对另一端边界处的函数值产生影响。对于一些非周期性的文艺,我们如何确保边界处的值对整体影响不大或者能确保边界处的值一直为一个特定的常数,是一个值得思考的问题。

    一个简单的数值算例

    一个简单的数值算例如下:

    在这里插入图片描述

    这里的R0=1/4R_0=1/4,我们计算t=3/256t=3/256时刻的uu值。假设它是周期性边界条件,将Ω\Omega扩充到R2R^2上,那么有:

    (Δu)^(k)=ΩΔu(x)eikxdx=Ωuηeikxds+ikΩu(x)eikxdx=ikΩu(x)eikxdx=ikΩu(x)eikxds+(ik)2Ωu(x)eikxdx=k2Ωu(x)eikxdx=k2u^(x)\hat{(\Delta u)}(k)=\int _{\Omega}\Delta u(x)e^{-ikx}dx\\ =\int_{\partial \Omega}\frac{\partial u}{\partial \eta}e^{-ikx}ds+i|k|\int_{\Omega }\nabla u(x)e^{-ikx}dx\\ =i|k|\int_{\Omega }\nabla u(x)e^{-ikx}dx\\ =i|k|\int_{\partial \Omega}u(x)e^{-ikx}ds+(i|k|)^2\int_{\Omega }u(x)e^{-ikx}dx\\ =-|k|^2\int_{\Omega }u(x)e^{-ikx}dx= -|k|^2\hat u(x)

    注意到,这里有一个问题,就是要使上式成立,要保证在xx \rightarrow -\infty时,有u(x)0u(x)\rightarrow 0,并且趋于0的速度要大于指数速率,才能保证边界项为0。

    如果对二维格林公式的操作不熟,我们也可以将uu写为u(t,x,y)u(t,x,y),然后对于Δu(t,x,y)\Delta u(t,x,y)分别对xxyy做两次一维的傅里叶变换,结果是一样的。即:
    Δu^(k)=(k12+k22)u^(x)\hat{\Delta u}(k) = -(k_1^2+k_2^2)\hat u(x)

    那么,记u(0,x,y)=g(x)u(0,x,y) = g(x)对上述问题,做傅里叶变换,就得到了:

    u^t+(k12+k22)u^=0u^(0)=g^\frac{\partial \hat u}{\partial t}+(k_1^2+k_2^2)\hat u = 0\\ \hat u(0) = \hat g
    对于这样一个ODE问题,我们可以设计欧拉格式、龙哥库塔等方法求解。写出来一个matlab脚本如下:

    function T = Runge_Kutta_fft_spectral(D,h,tau)
    %谱方法求解二维热传导方程
    L=1; N=round(1/h);%划分
    x=L/N*[-N/2:N/2-1]; y=x;%拆成两个维度来
    kx=(2*pi/L)*[0:N/2-1 -N/2:-1]; ky=kx;%对应的频域
    [X,Y]=meshgrid(x,y);
    [kX,kY]=meshgrid(kx,ky);
    K2=kX.^2+kY.^2;%非顺序下的K方
    %初始条件
    D0 = D(1:end-1,1:end-1);
    u0 = D0;
    ut0 = fft2(u0);
    ut= ut0; ut_flatten = ut(:);%flatten
    K2_flatten = K2(:);
    %求解
    t=[0 tau];
    %options = odeset('RelTol',1e-5,'AbsTol',1e-8);
    [t,utsol]=ode45(@(t,ut_flatten) heat_equation_2D(t,ut_flatten,K2_flatten),t,ut_flatten);
    T_temp = ifft2(reshape(utsol(end,1:N^2),N,N));
    
    T_temp = [T_temp T_temp(:,1)];
    T = [T_temp;T_temp(1,:)];%周期性边界扩充
    end
    function dut=heat_equation_2D(t,ut_flatten,K2_flatten)
    dut=[-K2_flatten.*ut_flatten];
    end
    

    程序比较好懂,注意看到的是我们的kx,kyk_x,k_y的取值为(2*pi/L)*[0:N/2-1 -N/2:-1],是为了计算(k12+k22)u^(k_1^2+k_2^2)\hat u的时候,频域坐标对应的频谱值是正确的。

    以上就是求解PDE的傅里叶谱方法,还有一个需要注意的地方是,但时间上的导数不只一阶的时候,我们往往通过引入一些新的变量(正如处理高阶ODE那样),将高阶的ODE变成一阶的ODE方程组。在做傅里叶变换,从而使用傅里叶谱方法。举个简单的例子。

    对于无界区域上的二维波动方程:

    在这里插入图片描述

    a=1a=1,初始条件为:

    在这里插入图片描述

    引入函数v=utv=\frac{\partial u}{\partial t},进行降阶,得到:

    在这里插入图片描述

    对上式等号两边做傅里叶变换,得到:

    在这里插入图片描述

    再用数值方法求解此ODE方程组即可。

    其他一些东西

    • 滤波法
      对于刚性比较大的微分方程,所以的刚性大,就是扰动太大,变化太快,我们可以使用滤波法来降低刚性。

    • 谱求导矩阵
      谱方法还可以用来做插值。使用谱求导矩阵,可以确定插值过程中的插值函数,并将其用于计算离散数据的各阶导数、偏微分方程(组)的数值求解。

    • 切比雪夫谱方法
      不管是傅里叶谱方法还是谱求导周期插值方法,在进行数值求导时都暗含了周期性边界条件,很有局限性。对于一般的第一类、第二类和第三类边界条件问题就不好处理。所以就有了切比雪夫谱方法。也就是说,切比雪夫谱方法可以处理非周期区间内的问题。

    谱方法

    谱方法和有限元法的思想很类似,差别在于:谱方法以一系列全局连续的函数(可以是三角函数,也可以是多项式)的叠加来近似真实解,而有限元法则是使用分片的简单函数的叠加来近似近似真实解。即,有限元的插值函数,只在该单元内作用。

    谱方法方法是解偏微分方程的一种数值方法。其要点是把解近似地展开成光滑函数(一般是正交多项式)的有限级数展开式,即所谓解的近似谱展开式,再根据此展开式和原方程,求出展开式系数的方程组。对于非定常问题,方程组还同时间有关。

    谱方法实质上是标准的分离变量技术的一种推广。一般多取切比雪夫多项式和勒让德多项式作为近似展开式的基函数。对于周期性边界条件,用傅里叶级数和调和级数比较方便。谱方法的精度,直接取决于级数展开式的项数。

    一般说,谱方法远比普通一、二阶差分法准确。由于快速傅里叶变换之类的技术不断发展,谱方法的运算量越来越少,一般是很合算的。特别是对于二维以上的问题,用差分法计算必须设置足够多的网格点,造成计算量的增加,而用谱方法一般不需取太多的项就可得到较高精度的解。因此谱方法在计算流体力学复杂流场的问题中有广泛应用。

    展开全文
  • 为了准确预测和分析多重网格法用于大地电磁法正演计算的收敛效果,对多重网格法的收敛性进行了广义傅里叶谱分析。通常情况下,系数矩阵的傅里叶谱是复数,为了直观地判断、提取收敛信息,将转换到实数域。在实数域...
  • 傅里叶谱 def myfft(inp,delta): #输入的x为竖向的array from scipy.fftpack import fft,ifft inp=np.array(inp).reshape(-1,1) inpf=np.fft.fft(inp,axis=0) #默认的是对行向量fourier变...

    主要用于地震动处理种的代码:

    1. 傅里叶谱
          def myfft(inp,delta): #输入的x为竖向的array
              from scipy.fftpack import fft,ifft
              inp=np.array(inp).reshape(-1,1)
              inpf=np.fft.fft(inp,axis=0) #默认的是对行向量fourier变换,所以此处要定义下轴
              inpf=inpf*2/inp.shape[0]
              inpf2=inpf[0:int((inp.shape[0]+1)/2)]#取一半
              n_points=inpf2.shape[0]#地震波点数
              frequency=1/delta
              fseries=np.linspace(0,frequency/2,n_points)
              #     plt.plot(fseries,inpf2)
              return inpf2,fseries

       

    2. 功率谱
          def mypowverspectra(inp,delta):
              from scipy.fftpack import fft,ifft
              inp=np.array(inp).reshape(-1,1)
              inpf=np.fft.fft(inp,axis=0) #默认的是对行向量fourier变换,所以此处要定义下轴
              inpf=inpf*2/inp.shape[0]
              inpf2=inpf[0:int((inp.shape[0]+1)/2)]#取一半
              inpf2=np.abs(inpf2)**2
              n_points=inpf2.shape[0]#地震波点数
              frequency=1/delta
              fseries=np.linspace(0,frequency/2,n_points)
              #     plt.plot(fseries,inpf2)
              return inpf2,fseries

       

    3. 反应谱
      def myresponse(inpacc,delta,damp,Tmax):
          #Tmax是反应谱周期最大值
          inpacc=np.array(inpacc).reshape(-1,1)
          count=inpacc.shape[0]
          displace=np.zeros([count])
          velocity=np.zeros([count])
          absacce=np.zeros([count])
          ta=np.arange(delta,Tmax,delta)  # 频段范围-1/deltaHZ
          mdis=np.zeros([len(ta)])#只记录一种阻尼比的响应幅值
          mvel=np.zeros([len(ta)])
          macc=np.zeros([len(ta)])
          frcy=2*math.pi/ta
          damfrcy=frcy*np.sqrt(1-damp**2)
          e_t=np.exp(-damp*frcy*delta)
          s=np.sin(damfrcy*delta)
          c=np.cos(damfrcy*delta)
          d_f=(2*damp**2-1)/(frcy**2*delta)
          d_3t=damp/(frcy**3*delta)
          for i in range(len(ta)):
              A=np.zeros([2,2])
              A[0,0]=e_t[i]*(s[i]*damp/np.sqrt(1-damp**2)+c[i])
              A[0,1]=e_t[i]*s[i]/damfrcy[i]
              A[1,0]=-frcy[i]*e_t[i]*s[i]/np.sqrt(1-damp**2)
              A[1,1]=e_t[i]*(-s[i]*damp/np.sqrt(1-damp**2)+c[i])
              B=np.zeros([2,2])
              B[0,0]=e_t[i]*((d_f[i]+damp/frcy[i])*s[i]/damfrcy[i]+(2*d_3t[i]+1/frcy[i]**2)*c[i])-2*d_3t[i]
              B[0,1]=-e_t[i]*(d_f[i]*s[i]/damfrcy[i]+2*d_3t[i]*c[i])-1/frcy[i]**2+2*d_3t[i]
              B[1,0]=e_t[i]*((d_f[i]+damp/frcy[i])*(c[i]-damp/np.sqrt(1-damp**2)*s[i])-(2*d_3t[i]+1/frcy[i]**2)*(damfrcy[i]*s[i]+damp*frcy[i]*c[i]))+1/(frcy[i]**2*delta)
              B[1,1]=e_t[i]*(1/(frcy[i]**2*delta)*c[i]+s[i]*damp/(frcy[i]*damfrcy[i]*delta))-1/(frcy[i]**2*delta)
              for k in range(0,count-1):
                  displace[k+1]=A[0,0]*displace[k]+A[0,1]*velocity[k]+B[0,0]*inpacc[k]+B[0,1]*inpacc[k+1]
                  velocity[k+1]=A[1,0]*displace[k]+A[1,1]*velocity[k]+B[1,0]*inpacc[k]+B[1,1]*inpacc[k+1]
                  absacce[k+1]=-2*damp*frcy[i]*velocity[k+1]-frcy[i]**2*displace[k+1]
              mdis[i]=np.max(np.abs(displace))
              mvel[i]=np.max(np.abs(velocity))
              if i==0:
                  macc[i]=np.max(np.abs(inpacc))
              else:
                  macc[i]=np.max(np.abs(absacce))
      #     plt.plot(ta,macc)
          return macc,ta

      反应谱是matlab代码修改的算法,详细见http://blog.sciencenet.cn/blog-419879-446739.html

    展开全文
  • 这里将结合opencv中的sample文档(\source\sample目录下的dft.cpp),对一副灰度图像进行傅里叶变换,并显示其傅里叶谱图像。 2. 傅里叶谱显示 #include #include int main(void) { //加载一副灰度图像 cv::...
  • 学习傅里叶图像变换的过程中,看到知乎上一个问题“为什么用图像二维傅里叶变换的相位进行反变换,能够大致得到原图的形状,而幅度则不行呢?”,然后思考了一下,理清楚思路,记录下来。 在傅里叶变换...
  • 傅里叶谱与矩不变结合的图像阈值分割,先对图像进行基于矩不变的阈值分割,在对分割后的图像进行傅里叶变换,对图像进行阈值调整
  • 关于傅里叶谱和相角的进一步说明: 在频域下,原始图像经过dft变换后,可得到两个矩阵,分别是复数的实部和虚部;Z=A+Bi 在复数坐标系下,某一点像素可由值和角度表示; 如果只保留相角信息,则只需将r=1,即...
  • 使用差分法和傅里叶谱方法求解拉普拉斯算子%Test laplacian operator with difference method and fourier spectral method clear;close all; Nx = 128;Ny=Nx;dx=1;dy=dx; %% Generage Gaussian Distribution X = 0 ...
  • 基于MATLAB的地震反应傅里叶谱计算分析XXX摘要:地震反应特征参数包括平台值和特征周期,其分别表征了地震动的强度特性和频谱特性;傅里叶谱表示地震波的重要意义有两个,一是检出时间过程中所含的频率分量,...
  • 傅里叶变换  二维DFT的极坐标表示  幅度或频率为 R(u,v)和I(u,v)分别是F(u,v)的实部和虚部  相角或相位为 ...
  • 转载自百度百科 52xenos |十级采纳率37% 冈萨雷斯版里面的解释非常形象:一个恰当的比喻是将傅里叶...当我们考虑光时,讨论它的光谱或频率。同样, 傅立叶变换使我们能通过频率成分来分析一个函数。 玻璃棱镜将自然
  • 频域滤波(一) 傅里叶谱与相位

    千次阅读 2014-05-09 15:25:02
    % 矩形的傅里叶变换2D phase_rec=angle(fft_rec);% 矩形的相位 rec_fre=abs(fft_rec);% 矩形的振幅 subplot(2,4,1),imshow(phase_rec,[]);title('ori rec phase 矩形的相位'); % figure,imshow(abs(fft_rec),[]),...
  • 提出傅里叶功率法评价分辨率标板极限分辨率值的理论判据。 这是一种用光学信息处理原理对标板图案傅里叶谱进行分析处理的客观评价方法。 文中建立了与分辨率标板图案相适应的傅里叶谱数学模型, 并通过对傅里叶谱的...
  • 方波的离散傅里叶变换分析,利用DTFT对方波进行傅里叶谱分析,并且设计恰当的数字滤波器,过滤一级谐波和二级谐波,很适合新人学习!!
  • 短时傅里叶变换估计
  • 傅里叶变换 相位 幅度

    万次阅读 2017-11-29 10:28:46
    周期信号的频谱  为了能既方便又明白地表示一个信号在不同频率下的幅值和相位,可以采用成为频谱图的表示方法。... 三角形式的傅里叶级数频率为非负的,对应的频谱一般称为单边;指数形式的傅里
  • //快速傅里叶变换 void fft2Image(InputArray _src, OutputArray _dst) { //得到Mat类型 Mat src = _src.getMat(); //判断位深 CV_Assert(src.type() == CV_32FC1 || src.type() == CV_64FC1); CV_Assert(src....
  • 整理的图像处理离散傅里叶变换频谱相位幅度关系ppt,及matlab代码
  • 本文基于二阶傅里叶微分矩阵来近似二阶导数,得到一个泊松方程的全离散傅里叶格式. 运用FFT理论分析了该数值格式,推导了快速方法,最后进行了数值试验. 数值试验显示数值方法求解速度快、方便实施,且高精度,...
  • 傅里叶级数(离散) [!] 以下笔记内容可能会出现部分错误的地方,恳请各位师生批评指正,谢谢! 主要内容 三角函数及正交性 傅里叶展开 频率与相位 傅里叶级数的复数形式 写在开头 由于目前我还是一名大一...
  • 在实验中用频谱分析仪器记录了法布里珀罗型光纤水听器透射光强的傅里叶变换,与软件仿真的理论谱线较好地吻合。通过实验进一步指出反馈临界点频率的选择方法,检验了反馈控制法在法布里珀罗型光纤水听器相位工作点...
  • 简单的求取下灰度图像的幅度和相位并进行双重构: 直接上代码: clear all Picture = imread('E:\others\Picture\Library.jpg'); Picture_Gray = rgb2gray(Picture);%灰度处理 Picture_FFT = fft2(Picture_...
  • 针对低信噪比条件下雷达信号脉内调制方式识别算法识别率低的问题,提出了一种基于分数阶傅里叶变换(FRFT)和循环的雷达信号识别方法。通过分数阶傅里叶变换搜索出最大峰值对应的分数阶,把信号粗分为非调频信号和...
  • 以语音信号的语图作为处理对象,提出了基于语图二次傅里叶变换对特定人二 字词汇识别的方法.首先对语图二次傅里叶变换频域图的图像意义以及相应的语音特性表征 进行了详细剖析;然后对语图频域图像进行二进...
  • 图像傅里叶变换,幅度,相位

    千次阅读 2018-09-09 10:36:20
    <span style="font-size:18px;">cl; img=imread('lena.jpg');... %傅里叶变换 f=fftshift(f); %使图像对称 r=real(f); %图像频域实部 i=imag(f); %图像频域虚部 margin=...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 48,113
精华内容 19,245
关键字:

傅里叶谱