精华内容
下载资源
问答
  • 空间傅里叶变换
    万次阅读
    2017-11-25 22:10:14

    空间域(spatial domain)

    空间域,或称图像空间,是以图像左上为原点,横为y竖为x的二维平面。
    这里写图片描述

    变换域

    在有些情况下,通过变换输入图像来表达处理任务,在变换域执行处理任务,然后再反变换到空间域会更好,比如使用傅里叶变换将图像转换到频率域进行滤波,再转换回空间域得到滤波后的图像。

    二维线性变换的通用形式可表示为:
    T ( u , v ) = ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) r ( x , y , u , v ) T(u, v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y)r(x, y, u, v) T(u,v)=x=0M1y=0N1f(x,y)r(x,y,u,v)
    其中, f ( x , y ) f(x, y) f(x,y)是输入图像, r ( x , y , u , v ) r(x,y,u,v) r(x,y,u,v)称为正变换核 T ( u , v ) T(u, v) T(u,v)称为 f ( x , y ) f(x, y) f(x,y)正变换

    给定 T ( u , v ) T(u, v) T(u,v)后,可以用 T ( u , v ) T(u, v) T(u,v)的反变换还原 f ( x , y ) f(x, y) f(x,y):
    f ( x , y ) = ∑ u = 0 M − 1 ∑ v = 0 N − 1 T ( u , v ) s ( x , y , u , v ) f(x, y) = \sum_{u=0}^{M-1}\sum_{v=0}^{N-1}T(u, v)s(x, y, u, v) f(x,y)=u=0M1v=0N1T(u,v)s(x,y,u,v)
    s ( x , y , u , v ) s(x, y, u, v) s(x,y,u,v)称为反变换核。两个式子一起称为变换对

    线性变换域中的一般操作方法
    如果
    r ( x , y , u , v ) = r 1 ( x , u ) r 2 ( y , v ) r(x, y, u, v) = r_1(x, u)r_2(y, v) r(x,y,u,v)=r1(x,u)r2(y,v)
    那么正向变换核是可分的。如果 r 1 ( x , u ) r_1(x, u) r1(x,u)的作用等于 r 2 ( y , v ) r_2(y, v) r2(y,v),则称变换核是对称的。从而有
    r ( x , y , u , v ) = r 1 ( x , u ) r 1 ( y , v ) r(x, y, u, v) = r_1(x, u)r_1(y, v) r(x,y,u,v)=r1(x,u)r1(y,v)
    如果 s s s同样适用,则同样适用于反变换核。

    傅里叶变换

    傅里叶变换在图像处理中是一种很常用的变换方法,可以使图像从空间域转换到频率域从而进行一些图像处理操作。

    二维离散傅里叶变换的变换核为
    r ( x , y , u , v ) = e − j 2 π ( u x M + v y N ) r(x, y, u, v) = e^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})} r(x,y,u,v)=ej2π(Mux+Nvy)
    s ( x , y , u , v ) = 1 M N e j 2 π ( u x M + v y N ) s(x, y, u, v) = \frac{1}{MN}e^{j2\pi(\frac{ux}{M} + \frac{vy}{N})} s(x,y,u,v)=MN1ej2π(Mux+Nvy)

    带入通用变换式中,可得离散傅里叶变换对
    T ( u , v ) = ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) e − j 2 π ( u x M + v y N ) T(u, v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y) e^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})} T(u,v)=x=0M1y=0N1f(x,y)ej2π(Mux+Nvy)
    f ( x , y ) = 1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 T ( u , v ) e j 2 π ( u x M + v y N ) f(x, y) = \frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}T(u, v)e^{j2\pi(\frac{ux}{M} + \frac{vy}{N})} f(x,y)=MN1u=0M1v=0N1T(u,v)ej2π(Mux+Nvy)
    傅里叶核是对称且可分的(证明见附),并且可分和对称的傅里叶核允许用一维傅里叶变换计算二维傅里叶变换(证明见附)。


    证明:傅里叶核的对称性和可分性

    根据同底幂的乘法运算 a m ∗ a n = a m + n a^m*a^n = a^{m+n} aman=am+n,可得
    r ( x , y , u , v ) = e − j 2 π ( u x M + v y N ) = e − j 2 π u x M e − j 2 π v y N = r 1 ( x , u ) r 2 ( y , v ) r(x, y, u, v) = e^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})} = e^{-j2\pi\frac{ux}{M}}e^{-j2\pi\frac{vy}{N}} = r_1(x, u)r_2(y, v) r(x,y,u,v)=ej2π(Mux+Nvy)=ej2πMuxej2πNvy=r1(x,u)r2(y,v)
    可分性证毕
    由上式可得
    r 1 ( x , u ) r 2 ( y , v ) = e − j 2 π u x M e − j 2 π v y N = r 1 ( x , u ) r 1 ( y , v ) r_1(x, u)r_2(y, v) = e^{-j2\pi\frac{ux}{M}}e^{-j2\pi\frac{vy}{N}} = r_1(x, u)r_1(y, v) r1(x,u)r2(y,v)=ej2πMuxej2πNvy=r1(x,u)r1(y,v)
    r 1 r_1 r1 r 2 r_2 r2是等价运算。对称性证毕

    证明:可分和对称的傅里叶核允许用一维傅里叶变换计算二维傅里叶变换

    T ( u , v ) = ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) e − j 2 π ( u x M + v y N ) = ∑ x = 0 M − 1 e − j 2 π x u M ∑ y = 0 N − 1 f ( x , y ) e − j 2 π y v N = ∑ x = 0 M − 1 T ( x , v ) e − j 2 π x u M T ( x , v ) = ∑ y = 0 N − 1 f ( x , y ) e − j 2 π y v N \begin{align} T(u, v) & = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y) e^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})}\\ & = \sum_{x=0}^{M-1}e^{-j2{\pi}\frac{xu}{M}}\sum_{y=0}^{N-1}f(x, y)e^{-j2{\pi}\frac{yv}{N}}\\ & = \sum_{x=0}^{M-1}T(x, v)e^{-j2{\pi}\frac{xu}{M}}\\ T(x, v) & = \sum_{y=0}^{N-1}f(x, y)e^{-j2{\pi}\frac{yv}{N}} \end{align} T(u,v)T(x,v)=x=0M1y=0N1f(x,y)ej2π(Mux+Nvy)=x=0M1ej2πMxuy=0N1f(x,y)ej2πNyv=x=0M1T(x,v)ej2πMxu=y=0N1f(x,y)ej2πNyv


    • 一维连续傅里叶公式:
      F ( ω ) = ∫ − ∞ ∞ f ( t ) e − j ω t d t F(\omega) = \int_{-\infty}^{\infty} f(t)e^{-j\omega t}dt F(ω)=f(t)etdt
      根据欧拉公式
      e − j θ = c o s ( θ ) − j s i n ( θ ) e^{-j\theta} = cos(\theta) - jsin(\theta) ejθ=cos(θ)jsin(θ)
      一维傅里叶公式等价于
      F ( ω ) = ∫ − ∞ ∞ f ( t ) [ cos ⁡ ( ω t ) − j sin ⁡ ( ω t ) ] d t F(\omega) = \int_{-\infty}^{\infty} f(t)[\cos(\omega t) - j\sin(\omega t)]dt F(ω)=f(t)[cos(ωt)jsin(ωt)]dt
    • 一维离散傅里叶公式:
      F ( ω ) = ∑ n = 0 N − 1 f ( n ) e − j ω n N F(\omega) = \sum_{n=0}^{N-1}f(n)e^{-\frac{j\omega n}{N}} F(ω)=n=0N1f(n)eNjωn
    更多相关内容
  • (很全)整理的离散傅里叶变换及图像处理磁共振k空间应用ppt、word及matlab代码
  • 空间傅里叶变换的基础上,提出一种用于柱状声源外部声场的声灵敏度计算方法,在已知振速及其对设计变量导数基础上,给出了可用于计算空间任意点的声学量关于声源的结构参数?材料特性等设计变量的灵敏度分析公式。...
  • 如果你在网上看过以下关于傅里叶变换的资料: 傅里叶分析之掐死教程(完整版)更新于2014.06.06 Python 中 FFT 快速傅里叶分析 形象理解二维傅里叶变换 【numpy】几种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. 观察窗口内的信号如果不是平稳的怎么办?

    1、傅里叶变换:信号处理的里程碑

    在19世纪初叶,法国数学家J. Fourier发现,任何周期函数都可以表示为复指数函数的无穷和。后来又经过许多数学家许多年的工作,慢慢形成了现在所为人熟知的傅里叶变换。

    一般情况下,若傅里叶变换不加任何限定语,则指的是连续傅里叶变换(连续函数的傅里叶变换)。

    定义傅里叶变换有许多种不同的形式,这里采用如下的定义:傅里叶变换将可积函数 x x x表示成复指数函数的积分形式。公式表示如下:

    X ( f ) = ∫ − ∞ ∞ x ( t ) ⋅ e − 2 π i f t d t X(f) = \int_{-\infty}^{\infty} x(t) \cdot e^{-2 \pi i f t} dt X(f)=x(t)e2πiftdt

    f f f为任意实数,在这里表示频率(单位是Hz), t t t表示时间(单位是s)。在适当的条件下,傅里叶变换是可逆的:

    x ( t ) = ∫ − ∞ ∞ X ( f ) ⋅ e 2 π i f t d f x(t) = \int_{-\infty}^{\infty} X(f) \cdot e^{2 \pi i f t} df x(t)=X(f)e2πiftdf

    注意, x x x X X X分别代表信号在时域和频域中的表示。

    我们首先看一下第一个公式:时域信号 x ( t ) x(t) x(t)乘上一个由指定频率组成的指数项( e − 2 π i f t e^{-2 \pi i f t} e2πift),然后在整个时间轴上进行积分,得到了频域信号。

    第一个公式中的指数项应用欧拉公式,可以有如下表示:

    e − 2 π i f t = c o s ( 2 π f t ) − i ⋅ s i n ( 2 π f t ) e^{-2 \pi i f t}=cos(2 \pi f t) - i\cdot sin(2 \pi f t) e2πift=cos(2πft)isin(2πft)

    上述公式有一个实部和一个虚部。所以,傅里叶变换的实质是将原始的时域信号乘上一个由固定频率构成的余弦与正弦函数后,再进行积分。如果积分的结果是一个比较大的值,则说明信号 x ( t ) x(t) x(t)在其频谱上有一个主要的分量 f f f;如果积分的结果是一个比较小的值,那么就说频率 f f f并不是信号 x ( t ) x(t) x(t)的主要组成频率;如果积分值是0,则说明信号中不存在频率 f f f

    注意到,积分是在整个时间轴上进行的。这也就是说,对于任一频谱分量 f f f,无论它是出现在某个特定的时间段,还是存在于信号的整个生命周期中,其对最终结果的影响是一样的。

    这回答了本文提出的第一个问题。

    总结一下,傅里叶变换仅能够识别出信号中的频谱分量有哪些,但无法进一步判断其出现的时间信息

    数学家们当然意识到了这个问题,并且也提出了一些解决方案,短时傅里叶变换就是其中的一种。

    2、短时傅里叶变换(The Short Time Fourier Transform, STFT)

    上一篇文章中已经提到,「短时」其实就是选择一个较短的时间窗口对信号进行观察,这个窗口是如此的短以至于我们可以认为观察到的信号是平稳的。

    信号平稳了,就可以对其应用傅里叶变换了。
    在这里插入图片描述

    图2.1

    上图的信号是通过程序生成的,这个人工信号从0时刻开始,每250ms区间内都是平稳的。如果选择的观察窗口不大于250ms,那么将窗口沿时间轴从0开始向右移动时,在某些情况下窗口内观察到的信号一定会出现平稳的情况。

    有读者可能要问了,如果选择的时间窗口是250ms,那么当把它移动到200-450ms区间上时,覆盖的信号是不平稳的,这该怎么办呢?答案:没有任何办法,只能对这部分信号也应用傅里叶变换。这里需要注意的是,窗口信号的平稳性是应用短时傅里叶变换的前提假设,如果不满足这个假设,那么得到的变换结果可能没有参考价值。

    这其实也就回答了第一篇遗留的两个问题:如果窗口信号不平稳,那也只能照样对其进行傅里叶变换;要求窗口“尽可能的小”其实没有一个统一的标准,这要视信号的频率而定。

    相较于基础的FT,短时傅里叶变换的公式中会多一个表示信号观察时间窗口函数 w ( ⋅ ) w(·) w() w ( ⋅ ) w(·) w()的取值必须(尽可能)保证在其区间内的信号是平稳的。

    由于「时间窗口」是其显著的特征,因此有人把STFT称为「窗口傅里叶变换」。

    3、从例子入手

    用到了python(3)的第三方库包括:numpy、matplotlib、scipy,需要自行安装。

    3.1 对典型的非平稳信号进行STFT

    现在,我们给出图2.1所示信号的生成代码:

    import numpy as np
    import matplotlib.pyplot as plt
    
    plt.rcParams['font.sans-serif'] = ['Arial Unicode MS'] 
    
    x = np.linspace(0, 1000, 1000)
    y1 = np.sin(2*np.pi*50*x[:250])
    y2 = np.sin(2*np.pi*30*x[250:250*2])
    y3 = np.sin(2*np.pi*10*x[250*2:250*3])
    y4 = np.sin(2*np.pi*5*x[250*3:])
    y = np.concatenate((y1, y2, y3, y4))
    
    plt.plot(x, y)
    plt.xlabel('时间[ms]')
    plt.ylabel('振幅')
    plt.show()
    

    由此易知,它的频率每250ms变化一次,从0时刻开始,其频率分别为50Hz、30Hz、10Hz和5Hz。

    接下来,我们对这个信号应用STFT并对结果进行可视化:

    import scipy
    # 应用STFT
    f, t, Zxx = scipy.signal.stft(y, fs=1000, nperseg=250)
    # 对结果进行可视化
    plt.pcolormesh(t*1000, f, np.abs(Zxx), vmin=0, shading='gouraud')
    plt.ylim(top=100)
    plt.ylabel('Frequency [Hz]')
    plt.xlabel('Time [ms]')
    plt.show()
    

    效果如下:
    在这里插入图片描述

    图2.2

    尽管这还是一个二维图像,但是通过颜色的深浅读者仍可以了解到第三个纬度的信息。从图中可以看出,信号的频率非常明显地分为4个部分,且各个频率分量在时间上的持续情况也被区分了出来

    读者需要注意,上图中的坐标轴上的数值并非是时间和频率的真实数值,而是经过了一定比例的缩放(通过调整scipy.signal.stft中的fs参数的值,你可以看到纵轴的变化)。

    3.2 对典型的平稳信号进行STFT

    在本小节,我们构造一个平稳的信号,并对它应用STFT以观察效果。代码如下:

    x = np.linspace(0, 1000, 1000)
    y1 = np.sin(2*np.pi*50*x)
    y2 = np.sin(2*np.pi*30*x)
    y3 = np.sin(2*np.pi*10*x)
    y4 = np.sin(2*np.pi*5*x)
    y = y1 + y2 + y3 + y4
    # 原始信号图
    plt.plot(x, y)
    plt.show()
    # 频域图
    f, t, Zxx = scipy.signal.stft(y, fs=1000, nperseg=250)
    plt.pcolormesh(t*1000, f, np.abs(Zxx), vmin=0, shading='gouraud')
    plt.ylim(top=100)
    plt.ylabel('Frequency [Hz]')
    plt.xlabel('Time [ms]')
    plt.show()
    

    结果如下:

    在这里插入图片描述

    图2.3

    在这里插入图片描述

    图2.4

    从频域图上可以看出,这个信号在整个时间轴上都有4个频率值。

    3.3 窗口对结果的影响

    接下来,我们将使用图2.1中的信号来观察一下不同的窗口宽度对变换结果的影响。

    本节选用了高斯窗口函数,在调用时需要额外传入其标准差以控制窗口宽度:

    f, t, Zxx = scipy.signal.stft(y, window=('gaussian', 80), fs=1000)
    

    我们知道,标准差越大,整个高斯函数的图像越扁平,对应到STFT中,则体现为在时间轴上的跨度越大。按照上面的猜想,时间跨度越大,则频率分辨率越好、时间分辨率越差。为此,我们设置候选的标准差为[10,40,80],分别查看一下其转换结果:

    f, t, Zxx = scipy.signal.stft(y, window=('gaussian', 10), fs=1000)  # 将此处的10分别换成40、80即可调整标准差
    plt.pcolormesh(t*1000, f, np.abs(Zxx), vmin=0, shading='gouraud')
    plt.ylim(top=100)
    plt.ylabel('Frequency [Hz]')
    plt.xlabel('Time [ms]')
    
    

    在这里插入图片描述

    图2.5,标准差为10

    在这里插入图片描述

    图2.6,标准差为40

    在这里插入图片描述

    图2.7,标准差为80

    从上面三个图可以看出,随着标准差的增大(窗口逐渐变宽),各个频率分量在频率轴上逐渐「聚焦」,即分辨频率的能力增加,不确定性降低。

    读者可以自行尝试继续增加窗口宽度,会发现频谱中各条频率“横线”逐渐横向拉长,彼此交叉。这说明了随着窗口的变宽,覆盖的时间范围增加,导致其在时间上的区分度降低,即时间分辨率降低。

    上述图像直观地展现了STFT所面临的选择困境。任何人在使用STFT时都要考虑:到底选择什么样的窗口函数?较窄的窗口函数具有较好的时间分辨率但较差的频率分辨率;较宽的窗口函数则恰恰相反。此外,较宽的窗口函数还会面临截取的信号不平稳的问题。

    因此,如何选择窗口函数就成了一个具体问题具体分析、见仁见智的问题了。至此,终于可以引出我们的小波变换了——它在一定程度上就是来解决这个选择困境的。

    展开全文
  • 用matlab实现快速傅里叶变换及其逆变换,另外对图像的频域进行高通滤波
  • 第14章:傅里叶变换

    千次阅读 2022-01-23 06:20:05
    第14章:傅里叶变换一、理论基础:二、Numpy实现傅里叶变换:1. 实现傅里叶变换:2. 逆傅里叶变换:3. 高通滤波示例:三、OpenCV实现傅里叶变换:1. 实现傅里叶变换:2. 实现逆傅里叶变换:3. 低通滤波示例: 图像...

    图像处理一般分为空间域处理和频率域处理。

    空间域:

    空间域处理是直接对图像内的像素点进行处理。空间域处理主要划分为灰度变换和空间滤波两种形式。灰度变换是对图像内的单个像素进行处理,比如调节对比度和处理域值等。空间滤波涉及图像的质量改变,比如图像平滑处理。空间域处理的计算简单方便,运算速度更快。

    频率域:

    频率域处理是先将图像变换的频率域,然后在频率域对图像进行处理,最后在通过反变换将图像从频率域变换到空间域。傅里叶变换是应用最广泛的一种频域变换,它能够将图像从空间域变换到频率域,而你傅里叶变换能够将频率域信息变换到空间域内。傅里叶变换在图像处理领域有着非常重要的作用。

    下面从理论基础、基本实现、具体应用等角度对傅里叶变回进行简单的介绍。

    一、理论基础:

    傅里叶变换非常抽象,很多人在工程中用了很多年的傅里叶变换也没有彻底了解傅里叶变回到底是怎么回事。为了更好的说明傅里叶变换,我们先看生活中的一个例子。

    下表所示是某饮料的配方,该配方是一个以时间形式表示的表格,表格很长,这里仅仅截取了其中一部分。该表中记录冲时刻“00:00”开始到某个特定时间“00:11”内的操作。

    image-20211206104821991

    分析表格发现,该配方:

    • 每隔1分钟放一块冰糖。
    • 每隔2分钟放3粒红豆。
    • 每隔3分钟放2粒绿豆。
    • 每隔4分钟放4块西红柿。
    • 每隔5分钟放1杯纯净水。

    上述文字是从操作频率的角度对配方进行说明。

    在数据处理过程中,经常使用图表的形式表述信息。如果从时域的角度,该配方表可以表示为如下图形式。

    image-20211206104741398

    如果从频率(周期)的角度表示,这个配方表可以表示为如下图形式,图中横坐标是周期(频率的倒数),纵坐标是配料的份数。

    image-20211206105207723

    对于函数,同样可以将其从时域变换到频域。下图是一个频率为5(1秒5个周期)、振幅为1的正弦曲线。

    image-20211206114814578

    如果从频率的角度考虑,则可以将其绘制为如下图所示的频域图。图中横坐标是频率,纵坐标是振幅。

    image-20211206115520854

    上述两图是等价的,他们是一个函数的两种不同表示方式。可以通过频域表示得到对应的时域表示,也可以通过时域表示得到对应的频域表示。

    法国数学家傅里叶指出,任何周期函数都可以表示为不同频率的正弦函数和的形式。在今天看来,这个理论是理所当然的,但是这个理论难以理解,在但是遭受了很大的质疑。

    下面我们来看傅里叶变换的具体过程。例如,周期函数的曲线如下图左上角所示。该周期函数可以表示为:

    • y = 3 * np.sin(0.8 * x) + 7 * np.sin(0.5 * x) +2 * np.sin(0.2 * x)

    因此,该函数可以看成是由下列三个函数的和构成的:

    • y1 = 3 * np.sin(0.8 * x)
    • y2 = 7 * np.sin(0.5 * x)
    • y3 = 2 * np.sin(0.2 * x)

    上述三个函数对应的函数曲线分别如下图右上角、左下角、右下角的图所示。

    image-20211206135247642

    如果从频域的角度考虑,上述三个正弦函数可以分别表示为下图中的三根柱子,图中横坐标是频率,纵坐标是振幅。

    image-20211206135428503

    通过以上分析可知,图中左上角的曲线可以表示为上图所示的频域图。

    从左上角的时域函数图形,构造出上图所示频域图像的过程,就是傅里叶变换。

    左上角的时域函数图形,与上图的频域图形表示的是完全相同的信息。傅里叶变换就是从频域的角度完整地表述时域信息。

    除了上述的频率和振幅外,还要考虑时间的问题。例如,饮料配方为了控制风味,需要严格控制加入配料的时间。上表中“00:00”时刻的操作,在更精细的控制下,实际上如下表所示。

    image-20211206141619089

    如果加入配料的时间发生了变换,饮料的风味就会发生变化。所以,在实际处理过程中,还要考虑时间差。这个时间差,在傅里叶变换中就是相位。相位表述的时域时间差相关的信息。

    例如,函数

    • y = 3 * np.sin(0.8 * x) + 7 * np.sin(0.5 * 2 + 2) + 2 * np.sin(0.2 * x + 2)

    可以看成下列三个函数的和的形式:

    • y1 = 3 * np.sin(0.8 * x)
    • y2 = 7 * np.sin(0.5 * x + 2)
    • y3 = 2 * np.sin(0.2 * x + 3)

    上述的四个函数分别对应的函数曲线为下图中的左上角、右上角、左下角、右下角。

    image-20211206142355153

    ​ 在本例中,如果把横坐标看成开始时间,则构成函数y的三个正弦函数并不都是从0时刻开始的,他们之间存在时间差。如果直接使用没有时间差的函数,则无法构成上图左上角所示的函数,而是会构成前面我们所说到的那个周期函数y = 3 * np.sin(0.8 * x) + 7 * np.sin(0.5 * x) +2 * np.sin(0.2 * x)。所以,相差是傅里叶变换中一个非常重要的条件。

    ​ 上面分别用饮料配方和函数的例子介绍了时域与频域转换的可行性,以帮助我们来理解傅里叶变换。

    在图形处理过程中,傅里叶变换就是将图像分解成正弦分量和余弦分量两部分,即将图像从空间域转换到频率域。 数字图像经过傅里叶变换后,得到的频域是复数。因此,显示傅里叶变换的结果需要使用实数图像(real image)加虚数图像(complex image),或者幅度图像(magnitude image)加相位图像(phase image)的形式。

    ​ 因为幅度图像中包含了原图我们所需要的大部分信息,所以图像处理过程中,通常仅使用幅度图像。当然,如果希望在频域内对图像进行处理,在通过逆傅里叶变换得到修改后的空域图像。就必须同时保留幅度图像和相位图像。

    对图像进行傅里叶变换后,我们会得到图像中低频和高频的信息。低频信息对应图像中变化缓慢的灰度分量。高频信息对应图像内变换越来越快的灰度分量,是由灰度的尖锐过渡造成的。例如,在一幅大草原的图像中,低频信息就对应着广袤的颜色趋于一致的草原等细节信息,而高频信息则对应着狮子的轮廓等各种边缘及噪声信息。

    傅里叶变换的目的,就是为了将图像从空域转换到频域,并在频域内实现对图像的特定处理,然后再对经过处理的频域图像进行逆傅里叶变换得到空域图像。 傅里叶变换在图像处理领域发挥着非常关键的作用,可以实现图像增强、图像去噪、边缘检测、特征提取、图像压缩和加密等。

    二、Numpy实现傅里叶变换:

    Numpy模块提供了傅里叶变换功能,Numpy模块中的fft2()函数可以实现图像的傅里叶变换。本节介绍如何用Numpy模块实现图像的傅里叶变换,以及在频域内过滤图像的低频信息,保留高频信息,实现高通滤波。

    1. 实现傅里叶变换:

    Numpy模块提供了实现傅里叶变换的函数numpy.fft.fft2(),它的语法格式是:

    • 返回值=numpy.fft.fft2(原始图像)

    这里需要注意的是,参数“原始图像”的类型是灰度图像,函数的返回值是一个复数数组(complex ndarray)。

    ​ 经过该函数的处理,就能得到图像的频谱信息。此时,图像频谱中的零频率分量位于频谱图像(频域图像)的左上角,为了便于观察,通常会使用numpy.fft.fftshift()函数将零频率成分移动到频域图像的中心位置,如图所示。

    image-20211207111412008

    函数numpy.fft.fftshift()的语法格式是:

    • 返回值=numpy.fft.fftshift(原始频谱)

    ​ 使用该函数处理后,图像频谱中的零频率分量会被移到频域图像的中心位置,对于观察傅里叶变换后频谱中的零频率部分非常有效。

    ​ 对图像进行傅里叶变换后,得到的是一个复数数组。为了显示为图像,需要将它们的值调整到[0,255]的灰度空间内,使用的公式为:

    • 像素新值=20*np.log(np.abs(频谱值))

    示例1:实现傅里叶变换

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    
    img = cv2.imread('../lena.bmp', 0)
    
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    magnitude_spectrum = 20 * np.log(np.abs(fshift))
    
    plt.subplot(1, 2, 1)
    plt.imshow(img, cmap='gray')
    plt.title('original')
    plt.axis('off')
    
    plt.subplot(1, 2, 2)
    plt.imshow(magnitude_spectrum, cmap='gray')
    plt.title('result')
    plt.axis('off')
    plt.show()
    

    image-20211207112001219

    2. 逆傅里叶变换:

    ​ 需要注意的是,如果在傅里叶变换过程中使用了numpy.fft.fftshift()函数移动零频率分量,那么在逆傅里叶变换过程中,需要先使用numpy.fft.ifftshift()函数将零频率分量移到原来的位置,再进行逆傅里叶变换。

    image-20211207112310325

    函数numpy.fft.ifftshift()是numpy.fft.fftshift()的逆函数,其语法格式为:

    • 调整后的频谱=numpy.fft.ifftshift(原始频谱)

    numpy.fft.ifft2()函数是numpy.fft.fft2()的逆函数。用来实现逆傅里叶变换,返回空域的复数数组。该函数的语法格式为:

    • 返回值=numpy.fft.ifft2(频域数据)

    函数numpy.fft.ifft2()的返回值仍旧是一个复数数组(complex ndarray)。

    逆傅里叶变换得到的空域信息是一个复数数组,需要将该信息调整至[0,255]灰度空间内,使用的公式为:

    • iimg=np.abs(逆傅里叶变换结果)

    示例1:实现逆傅里叶变换

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    
    img = cv2.imread('../boat.512.tiff', 0)
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    magnitude_spectrum = 20 * np.log(np.abs(fshift))
    
    ishift = np.fft.ifftshift(fshift)
    iimg = np.fft.ifft2(ishift)
    iimg = np.abs(iimg)
    
    plt.subplot(131)
    plt.imshow(img, cmap='gray')
    plt.title('img')
    plt.axis('off')
    
    plt.subplot(132)
    plt.imshow(magnitude_spectrum, cmap='gray')
    plt.title('fft2')
    plt.axis('off')
    
    plt.subplot(133)
    plt.imshow(iimg, cmap='gray')
    plt.title('iimg')
    plt.axis('off')
    plt.show()
    

    image-20211207114029406

    3. 高通滤波示例:

    在一幅图像内,同时存在着高频信号和低频信号。

    • 低频信号对应图像内变化缓慢的灰度分量。例如,在一幅大草原的图像中,低频信号对应着颜色趋于一致的广袤草原。
    • 高频信号对应图像内变化越来越快的灰度分量,是由灰度的尖锐过渡造成的。如果在上面的大草原图像中还有一头狮子,那么高频信号就对应着狮子的边缘等信息。

    滤波器能够允许一定频率的分量通过或者拒绝其通过,按照其作用方式可以划分为低通滤波器和高通滤波器。

    • 允许低频信号通过的滤波器称为低通滤波器。低通滤波器使高频信号衰减而对低频信号放行,会使图像变模糊。
    • 允许高频信号通过的滤波器称为高通滤波器。高通滤波器使低频信号衰减而让高频信号通过,将增强图像中尖锐的细节,但是会导致图像的对比度降低。

    傅里叶变换可以将图像的高频信号和低频信号分离。 例如,傅里叶变换可以将低频信号放置到傅里叶变换图像的中心位置,如前图所示,低频信号位于右图的中心位置。可以对傅里叶变换得到的高频信号和低频信号分别进行处理,例如高通滤波或者低通滤波。在对图像的高频或低频信号进行处理后,再进行逆傅里叶变换返回空域,就完成了对图像的频域处理。通过对图像的频域处理,可以实现图像增强、图像去噪、边缘检测、特征提取、压缩和加密等操作。

    ​ 例如,在下图所示,左图original是原始图像,中间的图像result是对左图original进行傅里叶变化后得到的结果,右图则是对result进行高通滤波后的结果。将傅里叶变换结果图像result中的低频分量值都替换为0(处理为黑色),就屏蔽了低频信号,只保留高频信号,实现高通滤波。

    image-20211207144328404

    要将上图中右图中间的像素值都置零,需要先计算其中心位置的坐标,然后选取以该坐标为中心,上下左右各30个像素大小的区域,将这个区域内的像素值置零。该滤波器的实现方法为:

    • rows,cols=img.shape

      crow,ccol=int(rows/2),int(cols/2)

      fshift[crow-30:crow+30,ccol-30:ccol+30]=0

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    
    img = cv2.imread('../boat.512.tiff', 0)
    
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    
    rows, cols = img.shape
    crow, ccol = int(rows / 2), int(cols / 2)
    fshift[crow - 30: crow + 30, ccol - 30: ccol + 30] = 0
    
    ishift = np.fft.ifftshift(fshift)
    iimg = np.fft.ifft2(ishift)
    iimg = np.abs(iimg)
    
    plt.subplot(121)
    plt.imshow(img, cmap='gray')
    plt.title('original')
    plt.axis('off')
    
    plt.subplot(122)
    plt.imshow(iimg, cmap='gray')
    plt.title('iimg')
    plt.axis('off')
    
    plt.show()
    

    image-20211207131139733

    三、OpenCV实现傅里叶变换:

    OpenCV提供了函数cv2.dft()和cv2.idft()来实现傅里叶变换和逆傅里叶变换。

    1. 实现傅里叶变换:

    OpenCV中使用函数cv2.dft()进行傅里叶变换,语法格式为:

    • 返回结果=cv2.dft(原始图像,转换标识)

    在使用该函数时,需要注意参数的使用规范:

    • 对于参数“原始图像”,要首先使用np.float32()函数将图像转换成np.float32格式。
    • “转换标识”的值通常为“cv2.DFT_COMPLEX_OUTPUT”,用来输出一个复数阵列。

    函数cv2.dft()返回的结果与使用Numpy进行傅里叶变换得到的结果是一致的,但是它返回的值是双通道的,第1个通道是结果的实数部分,第2个通道是结果的虚数部分。

    ​ 经过函数 cv2.dft()的变换后,我们得到了原始图像的频谱信息。此时,零频率分量并不在中心位置,为了处理方便需要将其移至中心位置,可以用函数numpy.fft.fftshift()实现。例如,如下语句将频谱图像 dft 中的零频率分量移到频谱中心,得到了零频率分量位于中心的频谱图像dftshift。

    • dftShift=np.fft.fftshift(dft)

    ​ 经过上述处理后,频谱图像还只是一个由实部和虚部构成的值。要将其显示出来,还要做进一步的处理才行。

    函数cv2.magnitude()可以计算频谱信息的幅度。该函数的语法格式为:

    • 返回值=cv2.magnitude(参数1,参数2)
      • 参数1:浮点型x坐标值,也就是实部。
      • 参数2:浮点型y坐标值,也就是虚部,它必须和参数1具有相同的大小(size值的大小,不是value值的大小)。

    函数cv2.magnitude()的返回值是参数1和参数2的平方和的平方根,公式为:

    image-20211207155205911

    I表示原始图像,dst表示目标图像。

    得到频谱信息的幅度后,通常还要对幅度值做进一步的转换,以便将频谱信息以图像的形式展示出来。简单来说,就是需要将幅度值映射到灰度图像的灰度空间[0,255]内,使其以灰度图像的形式显示出来。

    这里使用的公式为:

    • result=20*np.log(cv2.magnitude(实部,虚部))

    示例1:

    import cv2
    import numpy as np
    
    img = cv2.imread('../lena.bmp', 0)
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)
    
    result = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
    
    print(dft)
    print(dft_shift)
    print(result)
    

    image-20211207161323111

    示例2:用OpenCV函数对图像进行傅里叶变换,并展示其频谱信息。

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    
    img = cv2.imread('../lena.bmp', 0)
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)
    result = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
    plt.subplot(121)
    plt.imshow(img, cmap='gray')
    plt.title('original')
    plt.axis('off')
    
    plt.subplot(122)
    plt.imshow(result, cmap='gray')
    plt.title('result')
    plt.axis('off')
    
    plt.show()
    

    image-20211207162454629

    2. 实现逆傅里叶变换:

    在OpenCV中,使用函数cv2.idft()实现逆傅里叶变换,该函数是傅里叶变换函数cv2.dft()的逆函数。其语法格式为:

    • 返回结果=cv2.idft(原始数据)

    ​ 对图像进行傅里叶变换后,通常会将零频率分量移至频谱图像的中心位置。如果使用函数numpy.fft.fftshift()移动了零频率分量,那么在进行逆傅里叶变换前,要使用函数numpy.fft.ifftshift()将零频率分量恢复到原来位置。

    注意,在进行逆傅里叶变换后,得到的值仍旧是复数,需要使用函数cv2.magnitude()计算其幅度。

    示例:

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    
    img = cv2.imread('../boat.512.tiff', 0)
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)
    
    rst = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
    
    ishift = np.fft.ifftshift(dft_shift)
    iimg = cv2.idft(ishift)
    iimg = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])
    
    plt.subplot(131)
    plt.imshow(img, cmap='gray')
    plt.title('original')
    plt.axis('off')
    
    plt.subplot(132)
    plt.imshow(rst, cmap='gray')
    plt.title('rst')
    plt.axis('off')
    
    plt.subplot(133)
    plt.imshow(iimg, cmap='gray')
    plt.title('iimg')
    plt.axis('off')
    
    plt.show()
    

    image-20211207165841677

    3. 低通滤波示例:

    ​ 前面讲过,在一幅图像内,低频信号对应图像内变化缓慢的灰度分量。例如,在一幅大草原的图像中,低频信号对应着颜色趋于一致的广袤草原。低通滤波器让高频信号衰减而让低频信号通过,图像进行低通滤波后会变模糊。

    ​ 例如,在下图中,左图original是原始图像,中间的图像result是对original进行傅里叶变换后得到的结果,右图是低通滤波后的图像。将傅里叶变换结果图像result中的高频信号值都替换为0(处理为黑色),就屏蔽了高频信号,只保留低频信号,从而实现了低通滤波。
    在这里插入图片描述

    在实现低通滤波时,可以专门构造一个如下图左图所示的掩码图像,用它与原图的傅里叶变换频谱图像进行与运算,就能将频谱图像中的高频信号过滤掉。

    image-20211207171556251

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    
    img = cv2.imread('../boat.512.tiff', 0)
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)
    
    rows, cols = img.shape
    crow, ccol = int(rows / 2), int(cols / 2)
    mask = np.zeros((rows, cols, 2), np.uint8)
    mask[crow - 30: crow + 30, ccol - 30: ccol + 30] = 1
    
    f_shift = dft_shift * mask
    print(f_shift)
    
    # mask = np.zeros((rows, cols), np.uint8)
    # mask[crow - 30: crow + 30, ccol - 30: ccol + 30] = 1
    # f_shift = cv2.bitwise_and(dft_shift, dft_shift, mask=mask)
    # print(f_shift)
    
    ishift = np.fft.ifftshift(f_shift)
    iimg = cv2.idft(ishift)
    iimg = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])
    
    plt.subplot(121)
    plt.imshow(img, cmap='gray')
    plt.title('original')
    plt.axis('off')
    
    plt.subplot(122)
    plt.imshow(iimg, cmap='gray')
    plt.title('iimg')
    plt.axis('off')
    
    plt.show()
    

    image-20211207174941537

    展开全文
  • 资源是本人在计算机图像处理课程的结课作业做的一个图像傅里叶变换,可以完成把图像从空间域转换为频率域
  • 在利用空间调制傅里叶变换光谱仪对远距离目标进行光谱遥测时, 大气湍流扰动引起的入射光场的波前畸变会影响干涉图像和复原光谱的质量。根据大气湍流扰动对光场的相位调制作用, 建立了大气湍流的随机相位模型与光场在...
  • (2)计算离散图像的傅里叶变换; (3)掌握图像的傅里叶频谱图及离散傅里叶变换性质; (4)掌握 MATLAB中的傅立叶变换函数; (5)实现数字图像的傅立叶变换与反变换; (6)了解相位谱和幅值谱的意义,并且分别...
  • matlab离散傅里叶变换平滑代码令人敬畏的地理空间 数据库 -PostgreSql空间扩展。 -一组PostgreSQL函数,在创建矢量切片源时很有用。 -SQLite空间扩展。 -Neo4j的空间实用程序库。 -Oracle数据库空间扩展。 -MySql...
  • 文章目录傅里叶变换落地:离散傅里叶变换(DFT)1 傅里叶变换中连续到离散的演化1.1 由傅里叶变换(FT)演化出离散时间傅里叶变换(DTFT)1.2 离散时间傅里叶变换(DTFT)演化出离散傅里叶变换(DFT)1.3 傅里叶级数(FS)...
  • 在检测频率复杂的信号时,光学傅里叶变换的高带宽、快速性优势凸显,其在光谱和高速射频信号的实时检测中都有重要应用。基于色散的光学傅里叶变换方法包括色散展宽和时间透镜方法。基于光纤色散方法,对比分析了基于...
  • 16 分数阶傅里叶变换的程序FRFT 17 魔方模拟器与规划求解 18 隐马尔可夫模型工具箱 HMM 19 图理论工具箱GrTheory 20 自由曲线拟合工具箱ezyfit 21 分形维数计算工具箱FracLab 2.2 22 For-Each 23 PlotPub 24 ...
  • 傅里叶变换和线性空间

    千次阅读 2016-04-13 14:56:13
    如果是有限维的空间,这个线性变换就是矩阵,无限维空间,这个线性变换就是傅里叶变换,而且是线性空间中一个特殊的正交变换(三角函数是正交的),更浅薄一点的理解,傅里叶变换就是一个矩阵与向量之间的乘法,说...
  • 图像离散傅里叶变换,C++实现,可直接运行,调用OpenCV。
  • 快速傅里叶变换v9.0 IP核指南 ——Vivado设计套件 介绍:Xilinx FFT IP核是一种计算DFT的有效方式。 特点:•前向变换(FFT)和反向变换(IFFT)在复数空间,并且可以在运行的同时进行选择配置 •变换点数...
  • 傅里叶变换公式推导(一)

    千次阅读 2021-01-17 17:27:18
    而这神奇的函数变换规律,就来源于傅里叶变换,学习傅里叶变换,让我们透过现象看本质。 下面的图是由不同的正弦波所构成的矩形脉冲,它就像不同的齿轮相互嵌套旋转所形成。 二、完备正交函数集 我们先从空间中的一...
  • 此文档为个人原创首发,详细介绍了傅里叶变换与N维空间向量正交分解的对应关系,欢迎阅读,如有兴趣可与本作者联系,共同学习
  • 图像处理之图像傅里叶变换

    千次阅读 2022-08-21 09:30:11
    傅里叶变换是在以时间为自变量的“信号”与频率为自变量的“频谱”函数之间的某域研究中较复杂的问题在频域中变得简单起来,从而简化其分析过程;当自变量“时间”或“频率”为连续形式和离散形式的不同组合时,就...
  • 结果表明该法可以有效抑制传统傅里叶变换法处理光载波干涉条纹图时边缘效应所造成的较大误差,在基于空间域相位调制技术的波面干涉测量中,对空间载波条纹图进行处理,可以使相位的计算精度达到3.3 mrad。
  • 基于并行光学向量矩阵乘法器原理的离散傅里叶变换方法,利用光传播的高速和低损耗特点,提出以相位空间光调制器为核心变换矩阵的全光并行离散傅里叶变换方法,并通过实验进行了验证。实验结果显示,所提出的全光并行...
  • 点击上方蓝字关注我吧傅里叶变换是将图像从空间域转换到频率域,首先把图像波段转换成一系列不同频率的二维正弦波傅里叶图像;然后,在频率域内对傅里叶图像进行滤波、掩膜等各种操作,减少或者消除部分高频或低频...
  • 傅里叶变换 ~ 四种可能形式

    千次阅读 2019-10-10 15:49:04
    连续时间,离散频率的傅里叶变换——傅里叶级数 连续时间,连续频率的傅里叶变换——连续傅里叶变换 离散时间,连续频率的傅里叶变换——序列的傅里叶变换 离散时间,离散频率的傅里叶变换——离散傅里叶变换 ...
  • 由于用傅里叶变换轮廓术进行三维面形测量时, 测得的变形光场是空间有限函数, 故离散傅里叶变换时先要进行周期拓展, 如果拓展周期选择不当, 拓展后的条纹将不连续, 对之进行傅里叶变换会产生频谱泄漏。 文章从理论上...
  • 这篇文章介绍了首先介绍点的概念,从简单的点到复杂的点然后再从点到函数,接着引出傅里叶变换,介绍了它的优缺点,根据缺点提出改进的措施。
  • 最详细的图像傅里叶变换

    千次阅读 多人点赞 2022-05-04 10:41:50
    物理效果:傅里叶变换实现了将信号从空间域到频率域的转换 信号分析: 一维傅里叶变换(将杂乱的信号由时域转化到频域中)一维傅里叶变化是将信号分解为正弦波的和的形式 时域 横轴是时间,纵轴是振幅 频域 横轴是频率...
  • 傅里叶变换

    2022-05-16 17:12:10
    傅里叶变换是将按时间或空间采样的信号与按频率采样的相同信号进行关联的数学公式。在信号处理中,傅里叶变换可以揭示信号的重要特征(即其频率分量)。 对于包含 n 个均匀采样点的向量 x,其傅里叶变换定义为 ...
  • 傅里叶变换 傅氏变换的目的是讲函数整体从空域变换到频域,以便于作分析。它本身是一种线性变换。 F(μ)=∫−∞+∞f(t)∗e−2πμtdtF(\mu)=\int_{-\infty}^{+\infty} f(t)*e^{-2\pi\mu t}dtF(μ)=∫−∞+∞​f(t)∗e...
  • 采用坐标扩展变换,突破了快速傅里叶变换计算过程中输入屏和衍射屏空间尺度必须相同、抽样格点必须等间隔的限制,使上述问题得到解决,可以得到更加详细的聚焦光束光场分布。同时,采用分两步计算的思想,避免了计算...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 16,210
精华内容 6,484
关键字:

空间傅里叶变换

友情链接: myWebSite.zip