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

人在驴途

傅里叶变换是一种数学变换，它可以将一个函数或信号转换为另一个函数或信号，它可以将时域信号转换为频域信号，也可以将频域信号转换为时域信号。
在很多的领域都有广泛的应用，例如信号处理、通信、图像处理、计算机科学、物理学、生物学等。
它最大的功能是能够分析和提取信号的特征，将复杂的信号分解为简单的信号。
有人甚至说傅里叶变换是一种可以让我们看透世界本质的变换，将纷繁世界的表象中的本质显现出来。
1. 简介

从数学的角度来推导傅里叶变换的话，会涉及到很多的基本的数学概念，比如正交基，级数，正弦余弦函数等等。
这里主要是为了看懂和使用傅里叶变换，了解傅里叶变换中的主要概念即可，
Scipy库会帮助我们处理其中的细节。
举例来说，对于现实中遇到的如下一个原始信号（一般都是非周期性且无规律的）：

这样的信号，可以是声音，可以是气温变化，可以是心电图，甚至可以是股票涨跌情况等等。
原始的信号变化复杂，而且大部分情况都看不出什么规律。
这时，傅里叶变换就能发挥作用了，它能够将复杂无序的信号转换为一系列简单规则的信号。
比如上面的信号转换之后变成下面6个简单的信号：

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

Scipy中处理傅里叶变换有2个子模块：fft和fftpack。
fftpack即将被淘汰，所以尽量使用fft模块。
fft模块中，傅里叶变换的主要函数有：
函数名说明fft计算一维离散傅里叶变换ifft计算一维离散傅里叶逆变换fft2计算二维离散傅里叶变换ifft2计算二维离散傅里叶逆变换fftn计算N维离散傅里叶变换ifftn计算N维离散傅里叶逆变换2.1. 变换示例

创建一个复合的信号：

1. import numpy as np
2. import matplotlib.pyplot as plt
4. # 生成 0~8pi 之间100个点
5. x = np.linspace(0, 8*np.pi, 100)
7. # 随便生成6个不同的正弦信号
8. y1 = np.sin(x)
9. y2 = 4*np.sin(2*x)
10. y3 = 2*np.sin(4*x)
11. y4 = 8*np.sin(0.3*x)
12. y5 = 6*np.sin(0.8*x)
13. y6 = 0.5*np.sin(5*x)
15. y = y1 + y2 + y3 + y4 + y5 + y6
17. plt.plot(x, y)
18. plt.show()

复制代码

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

1. from scipy import fft as spfft
3. fft_result = spfft.fft(y)
4. print(fft_result.shape)
5. # 运行结果
6. (100,)
8. fft_result[1:11]
9. # 运行结果（显示前10个，总共100个）
10. array([ 99.548317    -0.j        , 273.43743482-274.69892934j,
11.        -10.91586207 +66.24943608j, 167.24328363-151.29814008j,
12.        -57.93543536 +49.80593353j, -30.50244642 +56.21303872j,
13.        -19.99077861 +37.24642604j, -13.54997057 +21.30562857j,
14.         35.119202  -161.91923934j, -17.17382544 +41.98143084j])
17. fft_result[-10:]
18. # 运行结果（显示后10个，总共100个）
19. array([-13.61273919 -29.66475047j, -17.17382544 -41.98143084j,
20.         35.119202  +161.91923934j, -13.54997057 -21.30562857j,
21.        -19.99077861 -37.24642604j, -30.50244642 -56.21303872j,
22.        -57.93543536 -49.80593353j, 167.24328363+151.29814008j,
23.        -10.91586207 -66.24943608j, 273.43743482+274.69892934j])

复制代码

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

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

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

• 频谱图：各个频率的波的振幅信息
  
• 相位图：各个频率的波的相位信息，相位就是波在特定时间所处的位置
  

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

1. fig = plt.figure(figsize=[8,4])
2. ax1 = fig.add_subplot(121)
3. data = np.abs(fft_result)
4. ax1.plot(data)
5. ax1.set_title("双边频谱图")
7. ax2 = fig.add_subplot(122)
8. data = np.abs(fft_result[:50])
9. ax2.plot(data)
10. ax2.set_title("单边频谱图")
12. plt.show()

复制代码

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

1. fig = plt.figure(figsize=[8,4])
2. ax1 = fig.add_subplot(121)
3. data = np.angle(fft_result)
4. ax1.plot(data)
5. ax1.set_title("双边相位图")
7. ax2 = fig.add_subplot(122)
8. data = np.angle(fft_result[:50])
9. ax2.plot(data)
10. ax2.set_title("单边相位图")
12. plt.show()

复制代码

从这两个图中就能得变换后各个频率的信号的频率，振幅，相位信息。
2.2. 逆变换示例

逆变换就是将变换后的信号还原为原始信号。
因为傅里叶变换并没有任何信息的损失，所以逆变换之后可以看出信号的波形没有任何改变。

1. data = spfft.ifft(fft_result)
3. # 逆变换之后，忽略虚部的数字
4. plt.plot(x, np.real(data))
5. plt.show()

复制代码

3. 应用

最后，利用傅里叶变换，我们试着做一个改变声音效果的小例子。
首先，读取一段音频，Scipy就可以直接读取wav文件。

1. import scipy.io.wavfile as wav
3. # 读取音频，返回采样率和采样的数值
4. rate, all_samples = wav.read("/path/to/fft-test.wav")
5. print(rate)
6. print(all_samples[1000:1010])
7. # 运行结果
8. 16000
9. [122 133 149 151 165 151 160 159 155 151]
11. plt.plot(all_samples)
12. plt.show()

复制代码

音频文件是网上随便找的一段英语对话。
接着，对读取的信号做傅里叶变换，观察变换后的结果：

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

复制代码

注意，这里用了 rfft函数，不是之前的fft函数，
两者的区别在于fft的结果是复数，会形成对称的双边图，就像上一节介绍的那样。
而rfft的结果是实数，且是单边的，如上图所示。
这两个函数根据实际情况选择使用，都可以对信号进行傅里叶变换。
因为我后面的处理不需要双边的信息，所以这里用 rfft 函数来做傅里叶变换。
3.1. 处理一

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

1. new_data = dd.copy()
2. # 20000HZ以上频率的数据设为0
3. new_data[20000:] = 0
5. fig = plt.figure(figsize=[8,4])
6. ax1 = fig.add_subplot(121)
7. ax1.plot(np.abs(new_data))
9. ax2 = fig.add_subplot(122)
10. ax2.plot(np.abs(dd))
12. plt.show()

复制代码

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

1. new_data = spfft.irfft(new_data)
2. new_data = np.rint(new_data)
3. new_data = new_data.astype("int16")
5. wav.write("/path/to/fft-test-1.wav", rate, new_data)

复制代码

转换后的音频和原音频比，听起来声音更加闷一些，模糊一些。
3.2. 处理二

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

1. new_data = dd.copy()
2. # new_data[20000:] = 0
3. new_data[:20000] = 0
5. fig = plt.figure(figsize=[8,4])
6. ax1 = fig.add_subplot(121)
7. ax1.plot(np.abs(new_data))
9. ax2 = fig.add_subplot(122)
10. ax2.plot(np.abs(dd))
12. plt.show()

复制代码

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

1. new_data = spfft.irfft(new_data)
2. new_data = np.rint(new_data)
3. new_data = new_data.astype("int16")
5. 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)

来源:https://www.cnblogs.com/wang_yb/p/17806842.html
免责声明：由于采集信息均来自互联网，如果侵犯了您的权益，请联系我们【E-Mail:cb@itdo.tech】 我们会及时删除侵权内容，谢谢合作！