信号处理 傅里叶变换 频谱分析

傅里叶主要针对的是平稳信号，简单说就是具有一定周期性的信号，

因为傅里叶变换采取的是有限取样的方式，所以对于 取样长度 和 取样对象 有着一定的要求。

FFT的基本用法

详见 参考资料1 

 

import numpy as np
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
mpl.rcParams['axes.unicode_minus'] = False  # 显示负号

'''
个人想法，有待验证：
采样不够如140，无法进行正确的分解，采样超出14000，不影响分解；（实验可得）
也就是说，当不确定频率或者采样点数时，如果采样点数增加不影响分解结果，说明采样点数足够了，（使用时的困惑）
'''

# 采样点选择1400个，
# 因为设置的信号频率分量最高为600赫兹，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400赫兹（即一秒内有1400个采样点，一样意思的）
N = 1400
x = np.linspace(0, 1, N)

# 设置需要采样的信号，频率分量有200，400和600
# y = 7 * np.sin(200 * x) + 5 * np.sin(400 * x) + 3 * np.sin(600 * x)
y = 7 * np.sin(2 * np.pi * 200 * x) + 5 * np.sin(2 * np.pi * 400 * x) + 3 * np.sin(2 * np.pi * 600 * x)
y += np.random.randn(N) * 10            # 加入噪声
plt.plot(x, y)
plt.show()


fft_y = fft(y)  # 快速傅里叶变换

x = np.arange(N)  # 频率个数
half_x = x[range(int(N / 2))]  # 取一半区间

abs_y = np.abs(fft_y)  # 取复数的绝对值，即复数的模(双边频谱)
angle_y = np.angle(fft_y)  # 取复数的角度
normalization_y = abs_y / N  # 归一化处理（双边频谱）
normalization_half_y = normalization_y[range(int(N / 2))]  # 由于对称性，只取一半区间（单边频谱）

plt.subplot(231)
plt.plot(x, y)
plt.title('原始波形')

plt.subplot(232)
plt.plot(x, fft_y, 'black')
plt.title('双边振幅谱(未求振幅绝对值)', fontsize=9, color='black')

plt.subplot(233)
plt.plot(x, abs_y, 'r')
plt.title('双边振幅谱(未归一化)', fontsize=9, color='red')

plt.subplot(234)
plt.plot(x, angle_y, 'violet')
plt.title('双边相位谱(未归一化)', fontsize=9, color='violet')

plt.subplot(235)
plt.plot(x, normalization_y, 'g')
plt.title('双边振幅谱(归一化)', fontsize=9, color='green')

plt.subplot(236)
plt.plot(half_x, normalization_half_y, 'blue')
plt.title('单边振幅谱(归一化)', fontsize=9, color='blue')

plt.show()

加入噪声

取样要求

FFT对于取样时间有要求，N点FFT 进行 精确频谱分析 的要求是 N个取样点包含 整数个 取样对象 的波形。

即 N点FFT 能够精确计算频谱 对 取样对象 的要求是 n*Fs/N（n*采样频率/FFT长度），

对8000HZ和512点而言,完美采样对象的周期最小要求是8000/512=15.625HZ，

所以156.25的n为10，234.375的n为15 (15.625*15=234.375)

import numpy as np
import pylab as pl

sampling_rate = 8000    # 取样频率
fft_size = 512          # FFT处理的取样长度

t = np.arange(0, 1.1, 1.0 / sampling_rate)
# np.arange(起点，终点，间隔)产生1s长的取样时间，两个正弦波叠加，156.25HZ和234.375HZ
x = np.sin(2 * np.pi * 156.25 * t) + 2 * np.sin(2 * np.pi * 234.375 * t)
# N为512时发生频谱泄露，100/(8000/512)=6.4，非整数，
# N为5120时不发生频谱泄露，100/(8000/5120)=64，整数
# N为5121时发生频谱泄露，100/(8000/5121)=64.0125，非整数
# x = np.sin(2 * np.pi * 100 * t) + 2 * np.sin(2 * np.pi * 234.375 * t)

xs = x[: fft_size]   # 从波形数据中取样fft_size个点进行运算
xf = np.fft.rfft(xs) / fft_size # 利用np.fft.rfft()进行FFT计算，rfft()是为了更方便

# 对实数信号进行变换，由公式可知/fft_size为了正确显示波形能量
# rfft函数的返回值是 N/2+1 个复数，分别表示从0(Hz)到sampling_rate/2(Hz)的分。
# 于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：

freqs = np.linspace(0, sampling_rate//2, fft_size//2+1)
# xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e1000))
xfp = np.abs(xf)
# 最后我们计算每个频率分量的幅值，并通过 20*np.log10()将其转换为以db单位的值。
# 为了防止0幅值的成分造成log10无法计算，我们调用np.clip对xf的幅值进行上下限处理


pl.subplot(211)
pl.plot(t[:fft_size], xs)
pl.xlabel(u"时间(秒)")
pl.title(u"The Wave and Spectrum 156.25Hz234.375Hz")

pl.subplot(212)
pl.plot(freqs, xfp, '.-')
pl.xlabel(u"Hz")
pl.subplots_adjust(hspace=0.4)
pl.show()

频谱泄露

当 其中一个  正弦波 频率改为 100，取样长度为 512 时，发生频谱泄露

 第一个正弦波的频谱出现泄露，能量分散到其他频率上，无法 精确 计算该信号 的频谱特性

参考资料：

https://cloud.tencent.com/developer/article/1441550

https://blog.csdn.net/weixin_42018112/article/details/98334685　　Python频域信号处理       和下面这篇文章内容差不多，但有些地方讲得更透

https://zhuanlan.zhihu.com/p/31962169　　基于Python的数字信号处理初步   介绍了 fft 的用法，频谱泄露 和 窗函数