傅里叶Fourier变换fft-python-scipy-幅值-辐角-相位(二）

文章目录

1 傅里叶级数
2 傅里叶变换
3 离散傅里叶变换DFT
重点：
举例1-尝试一个周期为2，采样时间为1的例子

重新梳理傅里叶变换时，在变换后的频率处理时，踩了一个大坑，记录一下。

1 傅里叶级数

傅里叶的发现总结成一般规律就是任意复杂的信号都能由一个个混合在一起的正弦函数的和来表示。

假设函数f(t)\Large f(t)f(t)定义在[0,T)上，则ω0=2πT\Large\omega_0 = \frac {2\pi} Tω0=T2π,
{sin⁡ω0t,cos⁡ω0t,...,sin⁡nω0t,cos⁡nω0t...\large\sin \omega_0 t,\cos \omega_0t, ...,\sin n\omega_0 t, \cos n\omega_0 t...sinω0t,cosω0t,...,sinnω0t,cosnω0t...} 是函数空间的标准正交基，
也可以写为复数形式，{ejω0t,...,ejnω0t,...\Large \displaystyle e^{j\omega_0 t}, ...,e^{jn\omega_0 t},...ejω0t,...,ejnω0t,...}

函数f(t)\Large f(t)f(t)的傅里叶级数，可以理解为在正交基上的投影之和

F(f(t))=∑n=−∞+∞Cnejnω0t=∑n=−∞+∞Cnej2nπTt\large \displaystyle \mathscr F(f(t))= \sum_{n=-\infty}^{+\infty}C_n e^{jn\omega_0 t} = \sum_{n=-\infty}^{+\infty} C_n \mathrm e^{j\frac {2n\pi} T t}F(f(t))=n=−∞∑+∞Cnejnω0t=n=−∞∑+∞CnejT2nπt，
其中，系数
cn=1T∫0Tf(t)e−j2nπTtdt\large \displaystyle c_n =\frac 1 T \int_0^T f(t) \mathrm e^{-j\frac{2n\pi} T t}dtcn=T1∫0Tf(t)e−jT2nπtdt

2 傅里叶变换

傅里叶变换，可以看作当周期趋于无穷大的情形下，傅里叶级数的扩展：

Ff(s)=∫n=−∞+∞f(t)ej2πstdt\large \displaystyle\mathcal Ff(s)= \int_{n=-\infty}^{+\infty}f(t) e^{j 2\pi s t}\mathrm dtFf(s)=∫n=−∞+∞f(t)ej2πstdt

其逆变换：

F−1f(t)=∫n=−∞+∞g(s)e−j2πstds\large \displaystyle\mathcal F^{-1}f(t)= \int_{n=-\infty}^{+\infty}g(s) e^{-j 2\pi s t}\mathrm dsF−1f(t)=∫n=−∞+∞g(s)e−j2πstds

傅里叶变换的定义是标准的，但不是唯一的。问题是把2π2\pi2π放在指数的哪里.
在T. W. Korner在他的书《傅里叶分析》中提供的总结，傅里叶变换可以写成以下形式：

Ff(s)=1AejBstf(t)\large\displaystyle\mathcal Ff(s)=\frac 1 A \mathrm e^{jBst}f(t)Ff(s)=A1ejBstf(t)

在实践中发现A与B 的选择是：
A=2πB=±1A=1B=±2πA=1B=±1\begin{aligned} \mathrm A &=\sqrt{2\pi}\quad &\mathrm B&= \pm 1\\ \mathrm A&=1 &\mathrm B&= \pm 2\pi\\ \mathrm A&=1 &\mathrm B&= \pm 1 \end{aligned} AAA=2π=1=1BBB=±1=±2π=±1
斯坦福教材¹选的定义是A=1， B=-2π2\pi2π。

3 离散傅里叶变换DFT

我们可以把离散傅里叶变换看作是一种操作，它接受N个数字作为输入，返回N个数字作为输出。

我们通常使用向量和“离散信号”表示法，并将n元组写成：
f=f[0],f[1],...,f[n−1]\displaystyle \pmb f = f[0], f[1],...,f[n-1]fff=f[0],f[1],...,f[n−1]

从信号f(t)f(t)f(t)及其傅里叶变换Ff(s)\mathcal Ff(s)Ff(s)开始，这两个函数都是连续变量。我们希望:

找到f(t)f(t)f(t)的离散版本这是f(t)f(t)f(t)的合理近似。
找到Ff(s)\mathcal Ff(s)Ff(s)的离散版本它是Ff(s)\mathcal Ff(s)Ff(s)的合理近似。
找到Ff(s)\mathcal Ff(s)Ff(s)与f(t)f(t)f(t)相关联的离散形式,这是Ff(s)\mathcal Ff(s)Ff(s)与f(t)f(t)f(t)相关联的合理近似。

假设f(t)f(t)f(t)在0<t<L0<t<L0<t<L之外时为零。我们也假设傅里叶变换Ff(s)\mathcal Ff(s)Ff(s)在0<s<2B0<s<2B0<s<2B之外时为零，或者实际上为零。

如果我们以每秒2B2B2B个样本的速率抽样，我们可以完美地从f(t)f(t)f(t)的样本中重建f(t)f(t)f(t)。因为f(t)f(t)f(t)定义在长度为LLL的区间内并且样本间隔是12B\displaystyle\frac 1{2B}2B1，这意味着我们想要的总数是:

N=L12B=2BL\displaystyle N=\frac L {\frac 1 2B}=2BLN=21BL=2BL (注意N为偶数)

均匀间隔的采样，这些点为：

t0=0,t1=12B,t2=22B,...,tN−1=N−12Bt_0=0, t_1=\frac 1{2B}, t2=\frac 2 {2B},...,t_{N-1}=\frac {N-1}{2B}t0=0,t1=2B1,t2=2B2,...,tN−1=2BN−1

知道f(tk)f(t_k)f(tk)就等于知道f(t)f(t)f(t)，这很合理，因此我们得到:

f(t)f(t)f(t)的离散版本是采样值列表f(t0),f(t1),…f(tN−1)f(t_0),f(t_1),…f (t_{N-1})f(t0),f(t1),…f(tN−1)。

接下来，表示f(t)f(t)f(t)的离散版本fdiscrete(tn)f_{discrete}(t_n)fdiscrete(tn)(采样值列表)在样本点上借助δ\deltaδ函数“连续”:

∑n=0N−1δ(t−tn)\displaystyle \sum_{n=0}^{N-1}\delta(t-t_n)n=0∑N−1δ(t−tn), 即

fdiscrete(t)=f(t)∑n=0N−1δ(t−tn)\displaystyle f_{discrete}(t)=f(t)\sum_{n=0}^{N-1}\delta(t-t_n)fdiscrete(t)=f(t)n=0∑N−1δ(t−tn)

这是我们之前考虑过的f(t)f(t)f(t)的采样形式，fdiscretef_{discrete}fdiscrete的傅里叶变换是

Ffdiscrete(s)=∑n=0N−1f(tn)Fδ(t−tn)=∑n=0N−1f(tn)e−2jπstn\displaystyle \mathcal Ff_{discrete}(s)=\sum_{n=0}^{N-1}f(t_n)\mathcal F\delta(t-t_n)=\sum_{n=0}^{N-1}f(t_n)e^{-2j\pi st_n}Ffdiscrete(s)=n=0∑N−1f(tn)Fδ(t−tn)=n=0∑N−1f(tn)e−2jπstn

这和我们想要的很接近，它是f(t)采样形式的连续傅里叶变换。

现在我们换个角度看频域的东西。
函数f(t)f(t)f(t)被限制为0≤t≤L0\le t\le L0≤t≤L，这决定了从它的频域重构Ff(s)\mathcal Ff(s)Ff(s)的采样率。采样间隔为1/L。

我们在频域从0到2B的区间内，在间隔为1/L1/L1/L的点上采样Ff(s)\mathcal Ff(s)Ff(s)。采样点的个数为

2B1/L=2BL=N\frac {2B}{1/L}=2BL=N1/L2B=2BL=N

和f(t)f(t)f(t)的采样点数目相同。Ff(s)\mathcal Ff(s)Ff(s)的采样点形式为m/Lm/Lm/L，共有NNN个采样点:

s0=0,s1=1L,...,sN−1=N−1Ls_0=0, s_1=\frac 1 L, ..., s_{N-1}=\frac {N-1}Ls0=0,s1=L1,...,sN−1=LN−1

重点：

默认情况下，采样时间是整周期的。
傅里叶变换描述的是针对周期函数或周期为无穷大的函数，
对于一个特定周期为T的函数，
如果采样时间为1.5个周期，那么傅里叶变换将以为1.5T为其周期，
这种情况下将很可能得不到想要的结果。

举例1-尝试一个周期为2，采样时间为1的例子

已知的函数很简单f(t)=cos⁡(πt)f(t)=\cos(\pi t)f(t)=cos(πt), 周期为2.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号# 采样时间
L = 1.0
# 采样点数
N = 128
# 采样点
t = np.linspace(0,L,N)
# 采样频率
Fs = N/L
# 函数为sin(4*np.pi*t)
y= np.cos(np.pi*t)

对 yyy 进行傅里叶变换，并作相应处理

# 对y进行傅里叶变换
y_fft = np.fft.fft(y)FFT得到的复数的模（即绝对值）就是对应的“振幅谱”，复数所对应的角度，就是所对应的“相位谱”
## 构造一个dataFrame
df_fft=pd.DataFrame(data=y_fft)
df_fft.columns=['变换值']
df_fft['幅值']=df_fft['变换值'].abs()
df_fft['幅值归一化']=df_fft['幅值']/N*2
df_fft.loc[0,'幅值归一化']=df_fft['幅值归一化'][0]/2
df_fft['辐角']=angle_y=np.angle(y_fft)
# 分辨率
dfs = 1/L
df_fft['频率']=np.arange(N)*dfs
df_fft['角频率']=df_fft['频率']*2*np.pi
pd.options.display.float_format = '{:.4f}'.format# 加入函数fft.fftfreq()频率
df_fft['函数直接计算频率']=np.fft.fftfreq(t.shape[-1])
# df_fft['函数直接计算频率*N']=df_fft['函数直接计算频率']*N
# 函数计算频率乘以采样频率
df_fft['函数直接计算频率*N/L']=df_fft['函数直接计算频率']*N/L
df_fft['函数直接计算角频率']=df_fft['函数直接计算频率*N/L']*2*np.pi

幅值最大的几个输出如下，

df_fft[df_fft.幅值归一化>0.1]变换值           幅值  幅值归一化   辐角  频率    角频率   直接频率    直接频率*N/L    函数直接计算角频率
1   1.3403-54.5962j 54.6126 0.8533  -1.5463 1.0000  6.2832  0.0078  1.0000      6.2832
2   1.0677-21.7332j 21.7594 0.3400  -1.5217 2.0000  12.5664 0.0156  2.0000      12.5664
3   1.0289-13.9489j 13.9868 0.2185  -1.4972 3.0000  18.8496 0.0234  3.0000      18.8496
4   1.0160-10.3159j 10.3658 0.1620  -1.4726 4.0000  25.1327 0.0312  4.0000      25.1327
5   1.0102-8.1902j  8.2522  0.1289  -1.4481 5.0000  31.4159 0.0391  5.0000      31.4159
6   1.0070-6.7887j  6.8629  0.1072  -1.4235 6.0000  37.6991 0.0469  6.0000      37.6991

可以看到，由于函数的周期为2，而采样时间为1，在这种采样情况下的频率分辨率为1Hz,
所以简单的余弦函数