【fourier】傅里叶变换



2018年01月06日    Author:Guofei

文章归类: 0x53_复分析与积分变换    文章编号: 92601

版权声明:本文作者是郭飞。转载随意,但需要标明原文链接,并通知本人
原文链接:https://www.guofei.site/2018/01/06/fourier.html


傅里叶展开

函数系

考察函数系:
$1, \cos t, \sin t,\cos 2t, \sin 2t, …, \cos nt, \sin nt,…$
这个函数系有一个性质 “正交性” ,任意两个不同的函数乘积在$[-\pi,\pi]$上的积分都是0.

卷积
$\int_{-\infty}^{+\infty}f_1(\tau)f_2(t-\tau)d\tau$

傅里叶级数的三角形式

狄利克雷条件
  1. 连续或只有有限个第一类间断点。2. 只有有限个极值点。

周期为$T$的周期函数$f(t)$,如果满足狄利克雷条件,就可以表示为:
$f(t)=\dfrac{a_0}{2}+\sum\limits_{n=1}^{+\infty}(a_n\cos n wt+b_n \sin nwt)$
其中,$w=\dfrac{2\pi}{T}$
$a_0=\dfrac{2}{T}\int_{-T/2}^{T/2} f(t) dt$
$a_n=\dfrac{2}{T}\int_{-T/2}^{T/2} f(t)\cos nwt dt$
$b_n=\dfrac{2}{T}\int_{-T/2}^{T/2} f(t)\sin nwt dt$
(被称为 欧拉-傅里叶公式
(可以由函数系的 正交性 推导出来)

性质

如果$f(t)$是偶函数,那么$b_n=0$
如果$f(t)$是奇函数,那么$a_n=0$

傅里叶级数的复指数形式

考虑到
$\cos nwt=\dfrac{1}{2}(e^{inwt}+e^{-inwt})$
$\sin nwt=\dfrac{1}{2}(e^{inwt}-e^{-inwt})$
得到
$f(t)=\dfrac{a_0}{2}+\sum\limits_{n=1}^{+\infty}(\dfrac{a_n-ib_n}{2}e^{inwt}+\dfrac{a_n+ib_n}{2}e^{-inwt})$
所以傅里叶变换也可以写为:
$f(t)=\sum\limits_{n=-\infty}^{+\infty}c_ne^{inwt}$
其中,$c_n=\dfrac{1}{T}\int_{-T/2}^{T/2}f(t)e^{inwt}dt(n=0,\pm1,\pm2,…)$

三角形式

令$A_0=a_0/2,A_n=\sqrt{a_n^2+b_n^2},\cos \theta_n=a_n/A_n,\sin\theta_n=\dfrac{-b_n}{A_n}$
那么,
$f(t)=A_0+\sum\limits_{n=1}^{+\infty}A_0\cos(nwt+\theta_n)$
$A_n$称为 振幅 ,$\theta_n$表示 相位

复指数形式中,$\mid c_n \mid=\mid c_{-n} \mid =A_n/2$就是振幅谱,$\arg c_n=-\arg c_{-n}=\theta_n$就是相位谱

傅里叶变换

连续傅里叶变换

当T增大时,$w=2\pi/T$越来越小,意味着频率的间隔越来越小;当$T\to+\infty$,频谱将变成一个连续值,现在我们分析这种情况
$f(t)=\lim\limits_{T\to+\infty}f_T(t)$
$=\dfrac{1}{2\pi}\int_{-\infty}^{+\infty}[\int_{-\infty}^{+\infty}f(\tau)e^{-iwt}d\tau]e^{iwt}dw$

傅里叶变换
$F(w)=\int_{-\infty}^{+\infty} e^{-iwt} f(t) dt$
傅里叶逆变换
$f(t)=\dfrac{1}{2 \pi} \int_{-\infty}^{+\infty} e^{iwt}F(w)dw$

有教材的另一种写法

傅里叶变换
$F(w)=\int_{-\infty}^{+\infty} f(t) e^{-i2\pi wt} dt$
傅里叶逆变换
$f(t)=\dfrac{1}{2 \pi} \int_{-\infty}^{+\infty} e^{i2\pi wt}F(w)dw$

$F(w)$是频谱密度函数,$\mid F(w)\mid$是振幅谱,$\arg F(w)$是相位谱

th

如果: $f(t)$是周期函数
那么,
$f(t)$可以表示为$f(t)=\sum\limits_{n=-\infty}^\infty C_n e^{-i\pi n t/p}$

$\delta$函数

$\delta$函数(单位冲击函数,狄拉克函数,Dirac)
满足这两个条件的函数: 1) \(\delta(t)= \left \{ \begin{array}{c} 0& t\neq 0 \\ \infty & t=0\end{array}\right.\) 2) $\int_{-\infty}^{+\infty}\delta(t)dt=1$

以上定义并不是严格定义,而是一种直观描述,$\delta$函数并不是经典意义上的函数,而是一种广义函数。

性质1
$f(t)$在R上有界,在$t=0$处连续,那么,
$\int_{-\infty}^{+\infty}\delta(t)f(t)dt=f(0)$

性质1又称为取样(sifting)特性
也可以写为 $\int_{-\infty}^{+\infty}f(t)\delta(t-t_0)dt=f(t_0)$

性质2
$\delta(t)=\delta(-t)$ (偶函数)
性质3
设$u(t)$为单位阶跃函数,即\(u(t)=\left\{\begin{array}{l}1,t>0\\0,t<0 \end{array}\right.\),那么有,
$\int_{-\infty}^t\delta(t)dt=u(t),\dfrac{d[u(t)]}{dt}=\delta(t)$
性质4
$\mathscr F(\delta(t))=1,\mathscr F^{-1}(1)=\delta(t)$

关于性质4,用$F(w)=1$画出形象示意图如下:

# dirac函数的傅里叶性质
# F(1)=\delta(t)的形象表示
import numpy as np
import matplotlib.pyplot as plt
w_list = np.arange(start=-5, stop=5, step=np.pi / 100)  # step定为一个无理数,最大程度模拟连续型w
x = np.arange(start=-100, stop=100, step=0.01)
y_sum = np.zeros_like(x)
for w in w_list:
    y = 1 * np.cos(w * x)
    y_sum += y
plt.plot(x, y_sum)
plt.show()

complexanalysis1

傅里叶变换的性质

  1. 线性性质 若$F(w)=\mathscr F(f(t)),G(w)=\mathscr F(g(t))$,那么,
    $\mathscr F(af(t)+bg(t))=aF(w)+bG(w)$
    $\mathscr F^{-1}(aF(t)+bG(t))=af(t)+bg(t)$
  2. 位移性质 若$F(w)=\mathscr F(f(t)),t_0,w_0$是实常数,那么
    $\mathscr F(f(t-t_0))=e^{-iwt_0}F(w)$
    $\mathscr F^{-1}(F(w-w_0))=e^{iw_0t}f(t)$
  3. 相似性 若$F(w)=\mathscr F(f(t)),a$是非0常数,那么
    $\mathscr F(f(at))=\dfrac{1}{\mid a\mid}F(\dfrac{w}{a})$
  4. 微分性质
    若$\lim\limits_{\mid t\mid \to +\infty}f(t)=0$,那么
    $\mathscr F(f^{(n)}(t))=(iw)^n \mathscr F(f(t))$
    $\dfrac{d^n F(w)}{dw^n}=(-i)^n \mathscr F(t^nf(t))$
  5. 积分性质
    设$g(t)=\int_{-\infty}^tf(t)dt,\lim\limits_{t\to+\infty}g(t)=0$,那么,
    $\mathscr F[f(t)]=\mathscr F[g’(t)]=iw\mathscr F[g(t)]$
  6. 帕塞瓦尔(Parseval)等式
    设$F(w)=\mathscr F(f(t))$,那么
    $\int_{-\infty}^{+\infty}f^2(t)=\dfrac{1}{2\pi}\int_{-\infty}^{+\infty} \mid F(w)\mid^2dw$

离散傅里叶变换

如果:
${ f_0,f_1,…f_{N-1} }$满足$\sum\limits_{n=0}^{N-1}|f_n|<\infty$

傅里叶变换: $X(k)=F(f_n)=\sum\limits_{n=0}^{N-1} f_n e^{-i \dfrac{2 \pi k}{N} n}$

傅里叶逆变换: $f_n=\dfrac{1}{N} \sum\limits_{k=0}^{N-1} X(k) e^{i \dfrac{2 \pi k}{N} n}$

冲击函数

离散系统中,也可以定义冲击函数 \(\delta(t)= \left \{ \begin{array}{c} 0& t\neq 0 \\ 1 & t=0\end{array}\right.\)

然后,可以得到与连续系统一致的性质

  • $\sum\limits_{x=-\infty}^\infty \delta(x)=1$
  • $\sum\limits_{x=-\infty}^\infty f(x)\delta(x)=f(0)$
  • $\sum\limits_{x=-\infty}^\infty f(x)\delta(x-x_0) = f(x_0)$

卷积

卷积是一种算子,定义为
$f(t) \ast h(t)=\int_{-\infty}^\infty f(\tau)h(t-\tau)d\tau$

直观理解:把一个函数做翻转,然后划过另一个函数。滑动过程中每个位移处执行乘积之和。

TH1:卷积可交换,就是 $f \ast g = g \ast f$

定理:如果$f(t),h(t)$的傅立叶变换是$F(u),H(u)$,那么卷积 $f(t) \ast h(t)$ 的傅立叶变换是 $F(u)H(u)$

  • 用各自的定义去推导,用到一个积分交换,不难。
  • 换句话说,两边可以相互获得 $f(t)\ast h(t) \Leftrightarrow F(u)H(u)$
  • 定理还有另一半。$f(t)h(t)\Leftrightarrow F(u)\ast H(u)$

取样

取样。用计算机处理之前,需要把连续函数转换成离散的序列。用语言描述取样,实际上是每隔一个间隔,取函数值。
数学表示: $f_k=\int_{-\infty}^\infty f(t) \delta(t-k\Delta T) dt=f(k\Delta T)$

如果取样后的结果记为$\tilde f$,对应傅立叶变换为$\tilde F$,

  • $\tilde f=f(t)s_{\Delta T}(t)$
  • 根据卷积定理$\tilde F=F(u)\ast S(u)$
  • $S(u)$事先可以算出来$S(u)=\dfrac{1}{\Delta T}\sum\limits_{n=-\infty}^\infty \delta(u-\dfrac{n}{\Delta T})$
取样定理
如果以超过函数最高频率的两倍的取样率来获取样本,连续的带限函数可以完全冲它的样本集中恢复。

还有一些结论:连续函数 $f$ 的取样是 $\tilde f$,对应的傅立叶变换 $F$ 的取样是 $\tilde F$,
那么把取样后的函数当成离散数列,恰好对应离散傅立叶变换 $\tilde F(u)=\sum\limits_{n\to -\infty}^\infty f_n \exp{(-j2\pi un\Delta T)}$
(其实不是恰好,而是说离散傅立叶变换就是这么定义的)

二维傅立叶变换

推广

定义连续二维冲击函数
\(\delta(t,z)=\left \{ \begin{array}{cc}\infty & t=z=0\\ 0 & o/w \end{array}\right.\)
并且
$\int_{-\infty}^\infty \int_{-\infty}^\infty \delta(t,z)dtdz=1$

然后有一系列类似一维的性质,不多说。

二维傅立叶变换:
$F(u,v)=\int_{-\infty}^\infty \int_{-\infty}^\infty f(t,z)\exp{(-j2\pi (ut+vz))} dtdz$
二维傅立叶逆变换:
$f(t,z)=\int_{-\infty}^\infty \int_{-\infty}^\infty F(u,v) \exp{(j2\pi (ut+vz))} dudv$

二维傅立叶变换也有对应的 取样定理,不多说。

关于采样,(在数字图像处理中),同样衍生处 混淆 的概念。

  • 一维混淆。包括空间混淆和时间混淆。都是欠取样造成的。
    • 空间混淆。例如线装特征中的锯齿、伪高光、原图像不存在的模式
    • 时间混淆,典型例子是,电影中车轮倒转的现象。

(《数字图像处理》冈萨雷斯版,提供了详尽的例子,可以去参考)

平移和旋转

$f(x,y)e^{j2\pi(u_0x/M+v_0y/N)} \Leftrightarrow F(u-u_0,v-v_0)$

$f(x-x_0,y-y_0)\Leftrightarrow F(u,v)e^{-j2\pi(x_0u/M+y_0v/N)}$

做极坐标变换。$x=r\cos \theta, y=r\sin\theta, u=w\cos \phi, v=w\sin \phi_0$

得到 $f(r,\theta+\theta_0) \Leftrightarrow F(w,\phi+\phi_0)$

也就是,如果 $f$ 旋转某角度,$F$ 也旋转同样的角度。

周期性

对称性

实用案例

import numpy as np

sample_rate = 500  # 采样率为 500Hz
t_max = 5  # 时间为5秒
t = np.linspace(start=0, stop=t_max, num=sample_rate * t_max, endpoint=False)  # 时间轴
# 创建信号:频率为 5 的信号 + 频率为 20 的信号
y = np.sin(2 * np.pi * 5 * t) + 0.8 * np.sin(2 * np.pi * 20 * t)

y_fft = np.fft.fft(y)
y_fft_magnitude = np.abs(y_fft)

画图展示 fft 和原始序列的关系

import matplotlib.pyplot as plt

fig = plt.figure(2)
axes = fig.subplots(nrows=2, ncols=1, sharex=True, sharey=False)

axes[0].plot(y)
axes[1].plot(y_fft_magnitude)

plt.show()

fft

计算周期(完整代码)

# %%计算最强的周期
len_y = len(y)
# 频率轴
freq = np.fft.fftfreq(len_y, 1 / sample_rate)

# 找到最大幅度对应的频率,只取左半边
peak_frequency = freq[np.argmax(y_fft_magnitude[:len_y // 2])]

# 周期是频率的倒数
period = 1 / peak_frequency
# 打印结果
print(f"主要周期: {period}秒")

# %%找到所有的周期。
all_possible = np.argsort(y_fft_magnitude[:len_y // 2])[::-1]

pred_i = all_possible[0]
for i in all_possible:
    # 如果强度比前一个小太多(这里是小一半),就停了
    if y_fft_magnitude[i] < y_fft_magnitude[pred_i] / 2:
        break
    print(f"周期:{1 / freq[i]},对应强度 {y_fft_magnitude[i]}")

结果:

主要周期: 0.2秒
周期:0.2,对应强度 1250.0
周期:0.05,对应强度 999.9999999999999

  • 如果 y 的均值不为0,那么 y_fft_magnitude[0] 代表的是直流通量。它会比较大,导致上面的代码报错。解决方法有2种
    • 跳过 all_possible[0]
    • y = y - y.mean() 使其均值归零

参考文献

张筑生:数学分析新讲
李红:《复变函数与积分变换》高等教育出版社
“十五”国家规划教材《复变函数与积分变换》高等教育出版社
钟玉泉:《复变函数论》高等教育出版社
冈萨雷斯:《数字图像处理》


您的支持将鼓励我继续创作!