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

最新推荐文章于 2022-10-30 23:23:36 发布

yingfly

最新推荐文章于 2022-10-30 23:23:36 发布

阅读量2k

点赞数

分类专栏：数据分析 #scipy 文章标签： python 信号处理傅立叶分析

本文链接：https://blog.csdn.net/yingfly/article/details/118145656

版权

数据分析同时被 2 个专栏收录

6 篇文章 1 订阅

订阅专栏

#scipy

5 篇文章 1 订阅

订阅专栏

文章目录

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

1 傅里叶级数

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

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

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

$\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}$ ，
其中，系数
$\large \displaystyle c_n =\frac 1 T \int_0^T f(t) \mathrm e^{-j\frac{2n\pi} T t}dt$

2 傅里叶变换

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

$\large \displaystyle\mathcal Ff(s)= \int_{n=-\infty}^{+\infty}f(t) e^{j 2\pi s t}\mathrm dt$

其逆变换：

$\large \displaystyle\mathcal F^{-1}f(t)= \int_{n=-\infty}^{+\infty}g(s) e^{-j 2\pi s t}\mathrm ds$

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

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

在实践中发现A与B 的选择是：
$\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}$
斯坦福教材¹选的定义是A=1， B=- $2\pi$ 。

3 离散傅里叶变换DFT

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

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

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

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

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

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

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

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

$t_0=0, t_1=\frac 1{2B}, t2=\frac 2 {2B},...,t_{N-1}=\frac {N-1}{2B}$

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

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

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

$\displaystyle \sum_{n=0}^{N-1}\delta(t-t_n)$ , 即

$\displaystyle f_{discrete}(t)=f(t)\sum_{n=0}^{N-1}\delta(t-t_n)$

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

$\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}$

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

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

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

$\frac {2B}{1/L}=2BL=N$

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

$s_0=0, s_1=\frac 1 L, ..., s_{N-1}=\frac {N-1}L$

重点：

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

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

已知的函数很简单 $f(t)=\cos(\pi 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)

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

# 对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,
所以简单的余弦函数 𝑐𝑜𝑠(𝜋𝑡) 未能被理想的复现，变成了多个谐波的叠加。

前6个谐波叠加与N/2个谐波叠加的图形与原始图形对比：

df_tmp = df_fft.loc[1:6]
yy=0
yy2=0
for i in range(1,7):
    yy += df_tmp.loc[i,'幅值归一化']*np.cos(df_tmp.loc[i,'角频率']*t+df_tmp.loc[i,'辐角'])
for i in range(1,int(N/2)):
    yy2 += df_fft.loc[i,'幅值归一化']*np.cos(df_fft.loc[i,'角频率']*t+df_fft.loc[i,'辐角'])    

plt.figure(dpi=100)
plt.plot(t,y,label='原始')
plt.plot(t,yy,label='前6项')
plt.plot(t,yy2,'--',label='前N/2项')
plt.legend()

图形对比
从上图可以看出，FFT后幅值较大前6项吻合性稍差，而前N/2项除端点外，已几乎与原曲线重合。

也可以看出，𝑐𝑜𝑠(𝜋𝑥) 与𝑐𝑜𝑠(2𝜋𝑥)在[0,1]上并不正交，
即函数组 { 𝑐𝑜𝑠(𝜋𝑥),𝑐𝑜𝑠(2𝜋𝑥),𝑐𝑜𝑠(3𝜋𝑥),…,𝑐𝑜𝑠(𝑛𝜋𝑥)…}的正交是在[0,2]区间上。

The Fourier Transform and its Applications. Brad Osgood. Electrical Engineering, Stanford University. ↩︎

yingfly

关注

0
点赞
踩
9

收藏

觉得还不错? 一键收藏
1
评论
傅里叶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,c
复制链接

扫一扫

专栏目录