傅里叶光学&MATLB编程系列【2】频谱的基本概念&空域时域到频域坐标的变换

0.前言

   这篇文章是整个系列第二篇,重点在于了解频谱的概念、空域坐标、时域坐标,以及空时域到频域坐标的变换,针对坐标问题,给出了相关得示例。下一篇文章讲述傅里叶变换的基本性质以及对应仿真代码的实现,并将会对一些光现象进行分析和解释。如果你是小白,建议从第一篇文章开始阅读
傅里叶光学&MATLAB编程系列【1】基本函数


1. 傅里叶级数

1.1 来源介绍

\qquad 对数学不感兴趣的读者可以跳过该部分,毕竟听到级数这两个字,我心里都难免一惊,这一切都得从大一时一本叫做高等数学的书说起,在这本书里面有很多让我们十分喜爱(痛恨 )的人,比如傅里叶和泰勒,他们的核心思想就是说满足一定条件的函数均可分解,傅里叶说可以分解为很多个正弦函数,泰勒说可以分解成幂函数
\qquad 傅里叶认为,如果函数 f ( x ) f(x) f(x)是一个周期为 T T T的函数,且满足一定条件,那么它就可以分解为傅里叶级数

f ( x ) = a 0 + ∑ n = 1 ∞ [ a n s i n ( 2 π n x T ) + b n c o s ( 2 π n x T ) ] (1) f(x)=a_0+\sum_{n=1}^\infty[ a_nsin(\frac{2\pi nx}{T}) +\\ b_n cos(\frac{2\pi nx}{T})] \tag{1} f(x)=a0+n=1[ansin(T2πnx)+bncos(T2πnx)](1)

\qquad 在上面公式中的系数 a 0 、 a n a_0 、 a_n a0an以及 b n b_n bn叫傅里叶系数。因为正弦、余弦函数是正交完备的,即不同周期的两个正余弦函数相乘之后在一个周期内积分为0,上式两边同时在一个周期内积分就可以消掉后面的求和部分,得到

a 0 = 1 T ∫ − T 2 − T 2 f ( x ) d x (2) a_0 =\frac{1}{T} \int_{-\frac{T}{2}} ^{-\frac{T}{2}}f(x)dx \tag{2} a0=T12T2Tf(x)dx(2)

同时剩上 s i n ( 2 π n x T ) sin(\frac{2\pi nx}{T}) sin(T2πnx)之后在一个周期内积分就可以消掉 b n b_n bn得到

a n = 1 T ∫ − T 2 − T 2 f ( x ) s i n ( 2 π n x T ) d x (3) a_n =\frac{1}{T} \int_{-\frac{T}{2}} ^{-\frac{T}{2}}f(x) sin(\frac{2\pi nx}{T})dx \tag{3} an=T12T2Tf(x)sin(T2πnx)dx(3)

同时剩上 c o s ( 2 π n x T ) cos(\frac{2\pi nx}{T}) cos(T2πnx)之后在一个周期内积分就可以消掉 a 0 a_0 a0 a n a_n an得到

b n = 1 T ∫ − T 2 − T 2 f ( x ) c o s ( 2 π n x T ) d x (4) b_n =\frac{1}{T} \int_{-\frac{T}{2}} ^{-\frac{T}{2}}f(x) cos(\frac{2\pi nx}{T})dx \tag{4} bn=T12T2Tf(x)cos(T2πnx)dx(4)

1.2傅里叶级数到频谱

\qquad 傅里叶级数到频谱的变换需要应用一个熟悉的欧拉公式

e j x = c o s ( x ) + j s i n ( x ) (5) e^{jx}=cos(x)+jsin(x) \tag{5} ejx=cos(x)+jsin(x)(5)

x x x变成 − x -x x得到

e − j x = c o s ( x ) − j s i n ( x ) (6) e^{-jx}=cos(x)-jsin(x) \tag{6} ejx=cos(x)jsin(x)(6)

将式(5)和式(6)相加、相减之后除于2得到

c o s ( x ) = e j x + e − j x 2 s i n ( x ) = e j x − e − j x 2 j (7) cos(x)=\frac{e^{jx}+e^{-jx}}{2} \\ sin(x)=\frac{e^{jx}-e^{-jx}}{2j} \tag{7} cos(x)=2ejx+ejxsin(x)=2jejxejx(7)

将式(7)代入式(1)得到

f ( x ) = a 0 + ∑ n = 1 ∞ [ 1 2 ( a n − j b n ) e j 2 π n x T + 1 2 ( a n + j b n ) e − j 2 π n x T (8) f(x)=a_0+\sum_{n=1}^\infty[ \frac{1}{2}(a_n-jb_n) e^{j\frac{2\pi nx}{T}}+\\ \frac{1}{2}(a_n+jb_n) e^{-j\frac{2\pi nx}{T}} \tag{8} f(x)=a0+n=1[21(anjbn)ejT2πnx+21(an+jbn)ejT2πnx(8)

c 0 = a 0 , c n = 1 2 ( a n − j b n ) , c n ∗ = 1 2 ( a n + j b n ) c_0=a_0,c_n= \frac{1}{2}(a_n-jb_n),c_n^*=\frac{1}{2}(a_n+jb_n) c0=a0,cn=21(anjbn),cn=21(an+jbn),上式的傅里叶级数就可以改写成

f ( x ) = ∑ n = − ∞ ∞ c n e j 2 π n x T (9) f(x)=\sum_{n=-\infty}^\infty{ c_n e^{j\frac{2\pi nx}{T}} } \tag{9} f(x)=n=cnejT2πnx(9)

将上式两边乘上 e j − 2 π n x T e^{j\frac{-2\pi nx}{T}} ejT2πnx,并在 [ − T 2 , T 2 ] [-\frac{T}{2},\frac{T}{2}] [2T,2T]积分得到

c n = 1 T ∫ − T 2 − T 2 f ( x ) e j − 2 π n x T d x (10) c_n=\frac{1}{T} \int_{-\frac{T}{2}} ^{-\frac{T}{2}}f(x) e^{j\frac{-2\pi nx}{T}}dx \tag{10} cn=T12T2Tf(x)ejT2πnxdx(10)

看到上面的公式,是不是觉得很熟悉,其实就是傅里叶变换的基本定义式。其中, 1 T \frac{1}{T} T1叫做 f ( x ) f(x) f(x)的基频,可以用 Δ u \Delta u Δu表示, u u u叫做函数 f ( x ) f(x) f(x)的谐频。

2. 频谱的概念

  一个物理量可以在空间域 x x x f ( x ) f(x) f(x)描述,可以在时间域 t t t f ( t ) f(t) f(t)来描述
相应地,将其变换到频域内, 也可以在空间频率域或者时间频率域 c n c_n cn来描述。这里的 c n c_n cn与公式(10)中的一致,由于 c n c_n cn表示频率为的谐波成分的复振幅,所以将 c n c_n cn按频率的分布图形称为 f ( x ) f(x) f(x)的频谱。由于 c n c_n cn一般是复函数,所以 c n c_n cn的模值 c n c_n cn随频率的分布图叫做 f ( x ) f(x) f(x)的振幅频谱,而 c n c_n cn的辐角随的分布图叫 做 f ( x ) f(x) f(x)的位相频谱。

3. 空域、时域及对应的坐标

3.1 一维离散空域坐标

  空域坐标用 x x x表示。在空域的描述是函数 f ( x ) f(x) f(x)随着空域坐标 x x x的变化,这里可以参考第一篇文章的相关内容:
矩形函数为例分析,需要在一个长为无穷大,宽为10mm的黑屏上开出一条宽为1mm的缝,应该如何描述黑屏的透过率 f ( x , y ) f(x,y) f(x,y)随空间坐标 x 、 y x、y xy的变换呢?在分析这一类问题时,首先需要建立坐标系,然后确定主要参数。比如,以x坐标表示宽,y坐标表示长,屏中心为原点建立坐标系,如下图所示:
在这里插入图片描述

  1. 分析问题:观察上图,我们可以发现屏的透过率不随坐标 y y y变化
    所以这个问题就变成了一维的问题
  1. 确定物的尺寸L:根据1可知,只分析x方向就行了。=物体的宽=度L的取值范围为[-5,5]mm。
  1. 矩形函数的宽度W:矩形孔的宽度W的取值范围为 [-0.5,0.5]mm。
  1. 采样数目M:M是所需要的采样点数,也就是创建的向量长度,这里不做多余考虑,取M为500😭事实上,采样的点数不能过少,根据奈奎斯特采样定律,采样太少,会发生混叠,信号失真,矩孔函数信号部分10个点以上,即 M ≥ 100 M\geq100 M100即可
  1. 采样间隔 Δ x = L / m = 10 / 500 = 0.02 m m \Delta x=L/m=10/500=0.02mm Δx=L/m=10/500=0.02mm;
  1. 坐标x: x ∈ [ − L / 2 : Δ x : L / 2 − Δ x x\in[-L/2:\Delta x:L/2-\Delta x x[L/2:Δx:L/2Δx]

对上面的例子分析过后就可以利用matlab绘制图形了,代码如下:

clc;
clear;
close all;%matlab的素质三连

L=10;%物体宽10mm
wx = 1;%缝宽1mm
M=500;%采样点数500
dx=L/M;%采样间隔
x=-L/2:dx:L/2-dx;%一维空间坐标
fx=(abs(x)<wx/2);
figure(1)
plot(x,fx);

在这里插入图片描述

3.2 二维离散空域坐标

   同样考虑上面提到的问题,如果我想要在一个长为50mm,宽为10mm的黑屏上开出一条长为5mm,宽为1mm的缝,应该如何描述黑屏的透过率 f ( x , y ) f(x,y) f(x,y)随空间坐标 x 、 y x、y xy的变换呢?
   一维到二维的问题很简单就可以依靠matlab提供的函数解决,因为一般初等函数都是可以分离变量的,即

f ( X , Y ) = f ( X ) ∗ f ( Y ) (11) f(X,Y)=f(X)*f(Y) \tag{11} f(X,Y)=f(X)f(Y)(11)

也就是可以分开考虑x和y。在建立坐标系之后,得到一维坐标 ( x , y ) (x,y) (x,y),然后用MATLAB提供的meshgrid函数就可以得到二维的坐标 ( X , Y ) (X,Y) (X,Y),再用二维坐标求 f ( X ) f(X) f(X) f ( Y ) f(Y) f(Y),相乘得到 f ( X , Y ) f(X,Y) f(X,Y)。代码如下:

clc;
clear;
close all;

long = 50;%50mm
wide = 10;%10mm
wx=5;%孔长5mm
wy=1;%孔宽1mm

M=500;%采样点数
dx=long/M;%x方向采样间隔
dy=wide/M;%y方向采样间隔

x=-long/2:dx:long/2-dx;
y=-wide/2:dy:wide/2-dy;

[X,Y] = meshgrid(x,y);
fX=abs(X)<wx/2;
fY=abs(Y)<wy/2;
fXY = fX.*fY;
mesh(x,y,fXY);
colorbar;

在这里插入图片描述

3.1 一维离散时域坐标

   时域坐标t(单位:s)通常用来描述信号得传输,无论是光信号还是通讯信号得秒速都离不开。根据傅里叶提出的分解思想,任何函数的基本模型都可以用正余弦函数描述,假设信号频率为 ω \omega ω(单位:Hz),振幅为 A m Am Am,初位相为 φ \varphi φ那么这个信号可以描述为

f ( t ) = A m ∗ s i n ( 2 π ω t + φ ) (12) f(t)=Am*sin(2 \pi \omega t+\varphi) \tag{12} f(t)=Amsin(2πωt+φ)(12)

在描述这个信号之前,通常需要考虑几个问题。

  1. 确定信号的时间长度。因为计算机的信息存储能力和计算能力有限,不可能直接描述信号在所有时间内的变化,通常情况下都只是取某一段时间 [ t 0 , t 1 ] [t_0,t_1] [t0t1]内的信号进行描述。
  2. 确定采样频率。这个问题可以用奈奎斯特采样定律解决,采样频率fs(单位:Hz)至少是 2 ω 2\omega 2ω才不会发生信号混叠现象
  3. 确定信号坐标。根据采样频率fs可以确定采样时间间隔为 Δ t = 1 / f s \Delta t=1/fs Δt=1/fs,时间坐标可以为 t = t 0 : Δ t : t 1 − Δ t t = t_0:\Delta t:t_1-\Delta t t=t0:Δt:t1Δt。其中, t 0 t_0 t0表示开始时刻, t 1 t_1 t1表示结束时刻。
  4. 描述信号。根据式(12)在matlab内绘图描述与观察。

举个例子,假设信号频率为 ω = 10 M \omega=10M ω=10M(单位:Hz),振幅为 A m = 1 Am=1 Am=1,初位相为 φ = 0 \varphi=0 φ=0,按照上述步骤分析:

  1. 确定信号的时间长度。取一段时间[0,1] μ s \mu s μs内的信号进行描述。
  2. 确定采样频率取采样频率fs=300M(单位:Hz),是 ω \omega ω的30倍,满足奈奎斯特采样定律,
  3. 确定信号坐标。根据采样频率fs可以确定采样时间间隔为 Δ t = 1 / f s \Delta t=1/fs Δt=1/fs,时间坐标可以为 t = 0 : Δ t : 1 e − 6 − Δ t t =0:\Delta t:1e^{-6}-\Delta t t=0:Δt:1e6Δt
  4. 描述信号。根据式(12)在matlab内绘图,代码如下:
clc;
clear;
close all;

t0 = 0;%开始时刻
t1 = 1e-6;%结束时刻
w = 10e6;%信号频率
fs = 300e6;%采样频率;
dt = 1/fs;%采样时间间隔
t = t0:dt:t1-dt;%时间坐标;

Am = 1;%信号幅度
varphi = 0;%初相位
ft = Am*sin(2*pi*w*t);
figure(1)
plot(t,ft);

在这里插入图片描述

4. 空域、时域坐标到频域坐标

4.1 空间频域坐标

  在分析信号频谱的时候,频域坐标是相当重要的,它会直接告诉你信号的组成成分。空域坐标用 x x x表示,其对应的频域坐标可以用 f x fx fx表示。如果知道了空域的坐标,其频域的坐标就可以计算出来,如果空域坐标

x ∈ [ − L / 2 : Δ x : L / 2 − Δ x ] (13) x\in[-L/2:\Delta x:L/2-\Delta x] \tag{13} x[L/2:Δx:L/2Δx](13)

其中L为信号的长度, Δ x \Delta x Δx为采样间隔,那么对应的频域坐标为

f x ∈ [ − 1 / ( 2 Δ x ) : 1 / L : 1 / ( 2 Δ x ) ] (14) fx\in[-1/(2\Delta x):1/L:1/(2\Delta x)] \tag{14} fx[1/(2Δx):1/L:1/(2Δx)](14)

其中的 f x fx fx单位为lp/mm,在光学中喜欢用cyc/mm描述,对应的坐标为

f x ∈ [ − π / Δ x : 2 π / L : π / Δ x ] (15) fx\in [-\pi / \Delta x : 2 \pi /L: \pi / \Delta x] \tag{15} fx[π/Δx:2π/L:π/Δx](15)

4.2 时间频域坐标

  时间频域坐标其实和空间频域坐标的分析基本一致,如果知道了时域坐标

t = t 0 : Δ t : t 1 − Δ t (16) t = t_0:\Delta t:t_1-\Delta t \tag{16} t=t0:Δt:t1Δt(16)

其中 ( t 1 − t 0 ) (t_1-t_0) (t1t0)为信号的长度,记作 t 10 t_{10} t10 Δ x \Delta x Δx为采样间隔,那么对应的频域坐标为

f t ∈ [ − 1 / ( 2 Δ t ) : 1 / t 10 : 1 / ( 2 Δ t ) ] (17) f t\in[-1/(2\Delta t):1/t_{10}:1/(2\Delta t)] \tag{17} ft[1/(2Δt):1/t10:1/(2Δt)](17)

当然,因为 Δ t = 1 / f s \Delta t=1/fs Δt=1/fs,也可以将其描述为

f t ∈ [ − f s / 2 : 1 / t 10 : f s / 2 ] (18) f t\in[-fs/2:1/t_{10}:fs/2] \tag{18} ft[fs/2:1/t10:fs/2](18)

下一篇文章讲述傅里叶变换的基本性质以及对应仿真代码的实现,将会告诉如何利用这一部分学习到的知识对信号和光现象进行分析和解释。

5、欢迎各位批评指正,我将在第一时间解答

  • 14
    点赞
  • 62
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

fightandstrive

创作不易,你的打赏,最大动力。

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值