【scipy 基础】--傅里叶变换

发布时间 2023-11-03 09:14:29作者: wang_yb

傅里叶变换是一种数学变换,它可以将一个函数或信号转换为另一个函数或信号,它可以将时域信号转换为频域信号,也可以将频域信号转换为时域信号。
在很多的领域都有广泛的应用,例如信号处理、通信、图像处理、计算机科学、物理学、生物学等。

它最大的功能是能够分析和提取信号的特征,将复杂的信号分解为简单的信号。
有人甚至说傅里叶变换是一种可以让我们看透世界本质的变换,将纷繁世界的表象中的本质显现出来。

1. 简介

从数学的角度来推导傅里叶变换的话,会涉及到很多的基本的数学概念,比如正交基,级数,正弦余弦函数等等。
这里主要是为了看懂和使用傅里叶变换,了解傅里叶变换中的主要概念即可,
Scipy库会帮助我们处理其中的细节。

举例来说,对于现实中遇到的如下一个原始信号(一般都是非周期性且无规律的):
image.png
这样的信号,可以是声音,可以是气温变化,可以是心电图,甚至可以是股票涨跌情况等等。

原始的信号变化复杂,而且大部分情况都看不出什么规律。
这时,傅里叶变换就能发挥作用了,它能够将复杂无序的信号转换为一系列简单规则的信号。
比如上面的信号转换之后变成下面6个简单的信号
image.png

另外,傅里叶变换是可逆的,它也能够将变换后的多个简单的信号还原成原始的复杂信号。

2. fft 模块

Scipy中处理傅里叶变换有2个子模块:fftfftpack
fftpack即将被淘汰,所以尽量使用fft模块。

fft模块中,傅里叶变换的主要函数有:

函数名 说明
fft 计算一维离散傅里叶变换
ifft 计算一维离散傅里叶逆变换
fft2 计算二维离散傅里叶变换
ifft2 计算二维离散傅里叶逆变换
fftn 计算N维离散傅里叶变换
ifftn 计算N维离散傅里叶逆变换

2.1. 变换示例

创建一个复合的信号:

import numpy as np
import matplotlib.pyplot as plt

# 生成 0~8pi 之间100个点
x = np.linspace(0, 8*np.pi, 100)

# 随便生成6个不同的正弦信号
y1 = np.sin(x)
y2 = 4*np.sin(2*x)
y3 = 2*np.sin(4*x)
y4 = 8*np.sin(0.3*x)
y5 = 6*np.sin(0.8*x)
y6 = 0.5*np.sin(5*x)

y = y1 + y2 + y3 + y4 + y5 + y6

plt.plot(x, y)
plt.show()

image.png

下面通过一维傅里叶变换,看看得到的结果:

from scipy import fft as spfft

fft_result = spfft.fft(y)
print(fft_result.shape)
# 运行结果
(100,)

fft_result[1:11]
# 运行结果(显示前10个,总共100个)
array([ 99.548317    -0.j        , 273.43743482-274.69892934j,
       -10.91586207 +66.24943608j, 167.24328363-151.29814008j,
       -57.93543536 +49.80593353j, -30.50244642 +56.21303872j,
       -19.99077861 +37.24642604j, -13.54997057 +21.30562857j,
        35.119202  -161.91923934j, -17.17382544 +41.98143084j])


fft_result[-10:]
# 运行结果(显示后10个,总共100个)
array([-13.61273919 -29.66475047j, -17.17382544 -41.98143084j,
        35.119202  +161.91923934j, -13.54997057 -21.30562857j,
       -19.99077861 -37.24642604j, -30.50244642 -56.21303872j,
       -57.93543536 -49.80593353j, 167.24328363+151.29814008j,
       -10.91586207 -66.24943608j, 273.43743482+274.69892934j])

观察快速傅里叶转换之后得到的结果:

  1. 转换的结果中,每个元素都是复数
  2. 前10个元素中,除去第一个,剩下的元素和后10个元素相比,实部相同,虚部相反

傅里叶变换之后,可以得到两个描述信号的图形:

  1. 频谱图:各个频率的波的振幅信息
  2. 相位图:各个频率的波的相位信息,相位就是波在特定时间所处的位置

频谱信息可以通过 np.abs 方法计算得出:

fig = plt.figure(figsize=[8,4])
ax1 = fig.add_subplot(121) 
data = np.abs(fft_result)
ax1.plot(data)
ax1.set_title("双边频谱图")

ax2 = fig.add_subplot(122) 
data = np.abs(fft_result[:50])
ax2.plot(data)
ax2.set_title("单边频谱图")

plt.show()

image.png

相位信息可以通过 np.angle 方法计算得出:

fig = plt.figure(figsize=[8,4])
ax1 = fig.add_subplot(121) 
data = np.angle(fft_result)
ax1.plot(data)
ax1.set_title("双边相位图")

ax2 = fig.add_subplot(122) 
data = np.angle(fft_result[:50])
ax2.plot(data)
ax2.set_title("单边相位图")

plt.show()

image.png

从这两个图中就能得变换后各个频率的信号的频率,振幅,相位信息。

2.2. 逆变换示例

逆变换就是将变换后的信号还原为原始信号。

因为傅里叶变换并没有任何信息的损失,所以逆变换之后可以看出信号的波形没有任何改变。

data = spfft.ifft(fft_result)

# 逆变换之后,忽略虚部的数字
plt.plot(x, np.real(data))
plt.show()

image.png

3. 应用

最后,利用傅里叶变换,我们试着做一个改变声音效果的小例子。

首先,读取一段音频,Scipy就可以直接读取wav文件。

import scipy.io.wavfile as wav

# 读取音频,返回采样率和采样的数值
rate, all_samples = wav.read("/path/to/fft-test.wav")
print(rate)
print(all_samples[1000:1010])
# 运行结果
16000
[122 133 149 151 165 151 160 159 155 151]

plt.plot(all_samples)
plt.show()

image.png
音频文件是网上随便找的一段英语对话。

接着,对读取的信号做傅里叶变换,观察变换后的结果:

dd = spfft.rfft(all_samples)
plt.plot(np.abs(dd))
plt.show()

image.png
注意,这里用了 rfft函数,不是之前的fft函数,
两者的区别在于fft的结果是复数,会形成对称的双边图,就像上一节介绍的那样。
rfft的结果是实数,且是单边的,如上图所示。

这两个函数根据实际情况选择使用,都可以对信号进行傅里叶变换。
因为我后面的处理不需要双边的信息,所以这里用 rfft 函数来做傅里叶变换。

3.1. 处理一

第一种处理,我尝试把频率20000HZ以上的信号都去掉,看看音频的效果有什么变化。

new_data = dd.copy()
# 20000HZ以上频率的数据设为0
new_data[20000:] = 0

fig = plt.figure(figsize=[8,4])
ax1 = fig.add_subplot(121) 
ax1.plot(np.abs(new_data))

ax2 = fig.add_subplot(122) 
ax2.plot(np.abs(dd))

plt.show()

image.png

处理之后的信号逆变换为原始信号,再保存为wav音频文件,看音频的变化效果。

new_data = spfft.irfft(new_data)
new_data = np.rint(new_data)
new_data = new_data.astype("int16")

wav.write("/path/to/fft-test-1.wav", rate, new_data) 

转换后的音频和原音频比,听起来声音更加闷一些,模糊一些。

3.2. 处理二

这次与处理一相反,把20000HZ以下的信号去掉。

new_data = dd.copy()
# new_data[20000:] = 0
new_data[:20000] = 0

fig = plt.figure(figsize=[8,4])
ax1 = fig.add_subplot(121) 
ax1.plot(np.abs(new_data))

ax2 = fig.add_subplot(122) 
ax2.plot(np.abs(dd))

plt.show()

image.png

然后同样的把信号逆变换回去,并保存为wav文件。

new_data = spfft.irfft(new_data)
new_data = np.rint(new_data)
new_data = new_data.astype("int16")

wav.write("/path/to/fft-test-2.wav", rate, new_data) 

这次的声音听起来很遥远,像是长途电话的感觉。

4. 总结

本篇主要介绍了傅里叶变换是什么,以及如何使用Scipy库来进行傅里叶变换。
目的是了解和使用傅里叶变换,而不是从数学角度去推导傅里叶变换。

最后的小例子虽然简单,但见微知著,仅仅删除了一些频率的信号,声音效果就随之变化。
如果把不同声音的文件中,影响音调,音色的频率分析出来,就可以制作一个变声器,
把自己的声音变成男声女声,儿童老人等等。

进一步,用二维傅里叶变换的话,可以分析图像,把图像中的主要频率找出来,
深入下去既可以做图像修复,也可以做图像识别等等。

Scipy库中还有多维傅里叶变换,可以分析现实情况下更加复杂的信号。

5. 附录

文中用到的完整代码(ipynb格式)和一个音频文件,可以从下面的地址下载:
scipy-fft-sample.zip:
https://url11.ctfile.com/f/45455611-910542660-6d7e0f?p=6872
(访问密码: 6872)