梦家blog 傅里叶变换找周期 python代码 傅里叶变换

在达姆经过Pro.Zoubir的DSP洗礼之后，回来上Money老师的DSP，依然听着一脸懵逼。。嗐，很多原因就是学了忘吧，于是打算记录一下。这篇blog就专门用来记录关于傅里叶变换的一些知识点与坑。

文章目录

• 1 四种傅里叶变换与DFT及FFT的关系
• 2 傅里叶变换中时域与频域信号强度对应关系的问题

• 2.1 问题描述：
• 2.1 连续信号
• 2.2 离散信号
• 2.3 实验检验

• 2.3.1 对一个两种频率的正弦信号进行fft
• 2.3.2 功率谱估计（周期图法）
• 2.3.3 整周期采样的影响
• 2.3.4我们来计算一下当我们采取1.2s，400Hz的采样频率的时候，我们可以准确得到哪些频率点的值

1 四种傅里叶变换与DFT及FFT的关系

这里我就直接用一张手写的笔记汇总了

基本思路就是从最原始的傅里叶级数，一步一步朝着实际使用里的离散和有限长度进行转化。

这个图片方向不对也是很难受，但我不会旋转过来，有看到的大哥求指导。

2 傅里叶变换中时域与频域信号强度对应关系的问题

2.1 问题描述：

大部分时候，我拿到信号之后就会去做傅里叶变换，做功率谱估计，去看看有什么频率，那个频率比较强（强度相对关系）。但是这个傅里叶变换的具体幅值（强度的绝对值）很多时候我都不知道对不对。所以这里主要针对这个问题进行分析。

2.1 连续信号

讨论连续信号，即讨论FS和FT，两个的傅里叶正变换的公式如下：
$FS: X(k\Omega_0)=\frac{1}{T}\int_{-T/2}^{T/2}x(t)e^{-jk\Omega_0 t}\\ FT: X(j\Omega)=\int_{-\infty}^{\infty}x(t)e^{-j\Omega t}d\Omega$
对于FS，其幅值为对应谐波的幅值，对于FT，由于频谱是连续的，所以对应的是频谱密度的概念
这个感觉就没什么好说的了，并不是很重要，因为只存在理论中

2.2 离散信号

这一part就比较重要了，主要关注实际使用的DFT。
$x(n) = \frac{1}{N}\sum^{N-1}_{k=0}X(k)e^{j\frac{2\pi}{N}nk}, n=0,1,2,...,N-1$
注意在DFT的反变换中，分解的各个频率分量前面还有一个1/N的系数，所以为了真实地反应信号的幅值，在做完DFT后要除以信号长度。
Note：DFT的来源为DFS在时域、频域都只取一个周期进行分析。
Note：由于DFT得到的频谱也是离散的，所以类比FS，得到的信号幅值也应该反映了该谐波信号的幅值。
Note：对应的功率谱分析则应该对应了该谐波的信号功率

2.3 实验检验

2.3.1 对一个两种频率的正弦信号进行fft

# 生成正弦周期信号
t = 1 # signal time in second
fs = 400 # sample frequency in Hz
n = np.arange(0, t, 1/fs)
N = len(n)
f1 = 4 # 4Hz
f2 = 25 # 25Hz
x = 2*np.sin(2 * np.pi * f1 * n) + np.sin(2 * np.pi * f2 * n)

X_w = np.fft.fft(x)/N # 傅里叶变换
X_w = np.abs(X_w)

1. 在计算fft的过程中，需要对计算得到的fft结果除以信号长度N，才可以得到正确的幅值。
2. 对于实信号，在对信号进行分解的时候，会分别投影到共轭的一对基向量上，反映到频谱上就是偶对称，因此，如果想要得到正确的幅值，需要把负半轴上的信号和正半轴上的信号相加（简单操作就是得到的幅值*2）才是真正的幅值。

2.3.2 功率谱估计（周期图法）

还是对于上面的信号，对其进行功率谱估计

1. 对于正弦信号，其功率是$A^2/2$(这里的功率求的时候是一个周期的能量除以一个周期的时间，既一个周期的平均功率，而不是瞬时功率)，也就是两个信号分别是2和1。也可以看出与图中的对应。
2. 同样地，在计算fft的时候也进行的除以信号长度的动作（实际代码用的是同一个fft的结果）
3. 对于实信号，也同样存在着上面偶对称的问题。

2.3.3 整周期采样的影响

看了上面两个情况，是不是觉得简直不能更简单了，那就Too young too naive了，下面我们来看如果我把上面的信号从1s，变成1.2s会发生什么。(因为这里信号比较简单，功率谱估计就是DFT的图平方，就只放一个DFT的频谱图了）

对照着上面看看，是不是发现形状也没那么好看了，而且幅值也变得奇奇怪怪。一个从1变成0.9几，一个却依然保持不变。那我们再来看看这个时候的信号在时域的样子。

似乎也没什么问题。

敲黑板 敲黑板 敲黑板，在上面的笔记中，当把离散非周期信号的DTFT进行频谱采样，或FS进行时间采样，转化成DFT时，都写的是要求每个周期采N个点，因此前提都是整周期采样。并且在进行DFT的时候是对所采样的信号进行了周期延拓，假设了采样的点为一个周期进行计算的。即：

上面的图为模拟的信号的真实在时间上的体现，下面的图为采样1.2s后进行周期延拓的结果，重点观察红色框框处，由于周期延拓，强行给信号加了一个新的为1.2s的周期，并且原来的4Hz和25Hz的信号的周期性也遭到了破坏。

2.3.4我们来计算一下当我们采取1.2s，400Hz的采样频率的时候，我们可以准确得到哪些频率点的值

信号的长度N = 480个点

所以对应的频率区间也有480个点

所以对应的频率为

$f = 400 \times \frac{n}{480}, n\in [0,479]$

对应地画在上面的得到的频谱图中即为：

上图生动地反应了FFT的栅栏效应，图中各横坐标以及图形表示了可以准确计算的频率。可以看到最接近4Hz的为4.167Hz，而25Hz则可以被准确地反应。因此，产生了上面做傅里叶变换时，4Hz的幅值发生了变形，而25Hz的依然保持了准确

并且，这一误差并不会随着你测量时间的增加而减小，下面四张图分别为上述信号在1.1s测量，1.2s测量，10.1s测量，10.2s测量下的傅里叶变换后绝对值，除以长度的结果。

• 随着时间从1.2s变为1.1s，25Hz处同样地出现了误差
• 时间从1s增加到10s，提高的是频率的分辨率，可以看到整个峰的宽度更窄了，但是峰值点的误差并没有改善

# 生成正弦周期信号
ts = [1.1, 1.2, 10.1, 10.2] # signal time in second
fig, axs = plt.subplots(2,2, figsize=(12,6))

for idx,t in enumerate(ts):
    fs = 400 # sample frequency in Hz
    n = np.arange(0, t, 1/fs)
    N = len(n)
    f1 = 4 # 4Hz
    f2 = 25 # 25Hz
    x = 2*np.sin(2 * np.pi * f1 * n) + np.sin(2 * np.pi * f2 * n)

    X_w = np.fft.fft(x)/N # 傅里叶变换
    X_w = np.abs(X_w)

    axs[idx//2, idx%2].plot(np.arange(0, fs, fs/N), X_w, 'r')
    axs[idx//2, idx%2].set_xlim(0,50)
    axs[idx//2, idx%2].set_ylim(0,1)
    axs[idx//2, idx%2].set_xlabel("Frequency[Hz]", fontsize=15)
    axs[idx//2, idx%2].set_ylabel("Amplitude", fontsize=15)
    axs[idx//2, idx%2].set_title("sample time t=%s"%ts[idx], fontsize=15)
    axs[idx//2, idx%2].grid(which='both', axis='y')

plt.tight_layout()

解决办法：

为了让fft的栅格准确地落在目标频率上，可以观察下式，那么需要采样频率fs/N后乘上某个整数Z，可以准确取到目标频率。因此，稍加变换，可以看到目标频率fs*测量时间t需要为某个整数Z。

因此，一个简单的规律可以总结为：

• 若目标频率是最小为为个位数，即1Hz，10Hz，101Hz等等，则采样时间应为1s的正整数倍。
• 若目标频率最小为小数点后一位，即0.1Hz，10.1Hz，101.1Hz等等，则采样时间应为10s的正整数倍。
• 以此类推
  $f_{signal} = Z \frac{f_s}{N}\\ f_{signal}*N*\frac{1}{f_s} = f_{signal}*t = Z$

所以，在进行傅里叶变换的时候，只有当采样点为整周期的时候，才能反正信号在该频率的真实幅值