大数据

Python科学计算——复杂信号FFT

FFT (Fast Fourier Transform, 快速傅里叶变换) 是离散傅里叶变换的快速算法,也是数字信号处理技术中经常会提到的一个概念。用快速傅里叶变换能将时域的数字信号转换为频域信号,转换为频域信号后我们可以很方便地分析出信号的频率成分。

单频信号FFT

# single frequency signal
sampling_rate = 2**14
fft_size = 2**12
t = np.arange(0, 1, 1.0/sampling_rate)
x = np.array(map(lambda x : x*1e3, t))
y = np.sqrt(2)*np.sin(2*np.pi*1000*t)
y = y + 0.005*np.random.normal(0.0,1.0,len(y))
# fft
ys = y[:fft_size]
yf = np.fft.rfft(ys)/fft_size
freq = np.linspace(0,sampling_rate/2, fft_size/2+1)
freqs = np.array(map(lambda x : x/1e3, freq))
yfp = 20*np.log10(np.clip(np.abs(yf),1e-20,1e100))

在对单频信号进行快速傅里叶变换的过程中,我们定义了两个常数: sampling_rate, fft_size,分别表示数字信号的采样频率和FFT的长度。由快速傅里叶变化的性质可知:采样频率 (sampling_rate) 确定的情况下,取波形中的 fft_size 个数据进行 FFT 变换时,若这 fft_size 个数据包含整数个周期, FFT 所计算的结果是精确的。即当被采样频率 f 满足如下公式时,FFT 的计算结果是精确的。

多频信号FFT

<1>双频信号FFT

# dual frequency signal
sampling_rate = 2**16
fft_size = 2**14
t = np.arange(0, 1, 1.0/sampling_rate)
x = np.array(map(lambda x : x*1000, t))
y = np.sqrt(2)*np.sin(2*np.pi*1000*t) + np.sqrt(2)*np.sqrt(2)*np.sin(2*np.pi*4000*t)
y = y + 0.005*np.random.normal(0.0,1.0,len(y))
# fft
ys = y[:fft_size]
yf = np.fft.rfft(ys)/fft_size
freq = np.linspace(0,sampling_rate/2, fft_size/2+1)
freqs = np.array(map(lambda x : x/1e3, freq))
yfp = 20*np.log10(np.clip(np.abs(yf),1e-20,1e100))

在对多频信号进行快速傅里叶变换的过程中,为了保证被采样数字信号的 FFT 计算结果精确,需要保证被采样数字信号中所有频率都是基频 (fo) 的整数倍,即满足如下公式:

相比于单频数字信号的快速傅里叶变换而言,多频数字信号的快速傅里叶变换有着更加苛刻的条件。除此之外,很多时候,我们并不知道被采样数字信号的所有频率成分,亦或我们只能在有限的时间段中对信号进行测量,无法知道在测量范围之外的信号是怎样的,因此只能对测量范围之外的信号进行假设。而快速傅里叶的假设很简单:测量范围之外的信号是所测量信号的重复。这就会引起在 fft_size 个点的采样范围内无法放下整数个所有被采样频率的波形而造成频谱泄漏

<2>频谱泄漏

当我们把双频信号FFT示例中的 fft_size 的值改为 2**12 时,这时,基频为 16Hz,不能被 1kHz整除,所以 1kHz 处发生了频谱泄露,而它能被 4kHz 整除,所以 4kHz 可以很好地被采样。

<2.1>频谱泄漏的原因

# 50Hz正弦波512点FFT采样
t = np.arange(0, 1.0, 1.0/8000)
x = np.sin(2*np.pi*50*t)[:512]
pl.plot(np.hstack([x,x,x]),linewidth=2)
pl.xlabel('Sampling point')
pl.ylabel('Amplitude[a.u.]')
pl.ylim(-1.5,1.5)
pl.show()

由于波形的前后不是连续的,出现波形跳变,而跳变处有着非常广泛的频谱,因此FFT的结果中出现了频谱泄漏。

<2.2>频谱泄漏的抑制

为了减小FFT所截取的数据段前后的跳变,可以对数据先乘以一个窗函数,使得其前后数据能平滑过渡。常用的hanning窗函数的定义如下:

pl.figure(figsize=(8,3))
pl.plot(signal.hann(512),linewidth=2)

512点 hann 窗函数

50Hz 正弦波与hann窗函数乘积之后的重复波形如下:

我们对频谱泄漏示例中的1kHz 和 4kHz 信号进行了 hann 窗函数处理,可以看出能量更加集中在 1kHz 和 4kHz,在一定程度上抑制了频谱泄漏。

<3>多点频信号FFT

以 1kHz 三角波为例,我们知道三角波信号中含有丰富的频率信息,它的傅里叶级数展开为:

<4>扫频信号FFT

当数字信号的频率随时间变化时,我们称之为扫频信号。以频率随时间线性变化的扫频信号为例,其数学形式如下:

其频率随时间线性变化,当我们在 [0,1] 的时间窗口对其进行采样时,其频率范围为 0~5kHz。当时间是连续时,扫频信号的频率也是连续的。但是在实际的处理中,是离散的点采样,因此时间是不连续的,这就使扫频信号的快速傅里叶变换问题退化为多点频信号快速傅里叶变换问题。其快速傅里叶变换得到的频谱图如下所示:

因上述扫频信号的频率随时间是线性变化的,因此功率谱是一条直线,即数字信号被采样频段范围内各个频率的功率是一样的。

<5>相位调制信号的FFT

以 50Hz 正弦信号相位调制到 1kHz 的信号为例,其信号形式如下:

它的时域波形,频率响应和相位响应如下图所示:

<6>FFT中的能量守恒

以扫频信号为例,当我们要探究FFT中的能量守恒时,我们要回归到信号最初的形式:

从频谱图上我们可以看出,上述数学形式的单频信号的功率为 -3dB。在扫频信号FFT示例中 5kHz 范围内信号的功率为 -40dB。它们之间的关系为:

当功率为 -3dB 的信号扩展到5kHz频谱范围的扫频信号时,其功率衰减到 -40dB,这也遵循了能量守恒定律。

Stay hungry, Stay foolish. — Steve Jobs