基于python(scipy和numpy)的快速傅里叶变换FFT

东方耀

基于python(scipy和numpy)的快速傅里叶变换FFT

1. import numpy as np
   
2. from scipy.fftpack import fft, ifft
   
3. import matplotlib.pyplot as plt
   
4. from matplotlib.pylab import mpl
   
5. # 基于python(scipy)的快速傅里叶变换FFT
   
6. mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
   
7. mpl.rcParams['axes.unicode_minus'] = False  # 显示负号
   
8. 
   
9. # fft简单滤波的步骤：
   
10. # 1、将原信号进行FFT
    
11. # 2、去掉需要滤波的频率
    
12. # 3、逆变换得到信号数据
    
13. 
    
14. # FFT——离散傅里叶变换（DFT）的快速算法。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的
    
15. # 采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性
    
16. # 由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果
    
17. 
    
18. # 采样点选择1400个，因为设置的信号频率分量最高为600赫兹，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400赫兹（即一秒内有1400个采样点，一样意思的）
    
19. # 一般实际应用中保证采样频率为信号最高频率的2.56～4倍
    
20. N = 1400
    
21. x_time = np.linspace(0, 1, N, endpoint=True)
    
22. # x_time = np.arange(0, 1, 1/N)
    
23. 
    
24. # 设置需要采样的信号，频率分量有200，400和600
    
25. y = 7 * np.sin(2 * np.pi * 200 * x_time) + 5 * np.sin(2 * np.pi * 400 * x_time) + 3 * np.sin(2 * np.pi * 600 * x_time)
    
26. 
    
27. print("变换之前，看看原始值：", y[:5])
    
28. fft_y = fft(y)  # 快速傅里叶变换
    
29. yreal = fft_y.real               # 获取实数部分
    
30. yimag = fft_y.imag               # 获取虚数部分
    
31. 
    
32. print("快速傅里叶变换后的长度=%d;看看前面的值=" % len(fft_y), fft_y[:5])
    
33. # 变换之后的结果数据长度和原始采样信号是一样的
    
34. # 每一个变换之后的值是一个复数，为a+bj的形式  复数具有模和角度
    
35. # “振幅谱”“相位谱”，它其实就是通过对快速傅里叶变换得到的复数结果进一步求出来的
    
36. # FFT得到的复数的模（即绝对值）就是对应的“振幅谱”，复数所对应的角度，就是所对应的“相位谱”
    
37. 
    
38. 
    
39. x = np.arange(N)  # 频率个数
    
40. half_x = x[range(int(N / 2))]  # 取一半区间
    
41. 
    
42. abs_y = np.abs(fft_y)  # 取复数的绝对值，即复数的模(双边频谱)   振幅谱
    
43. angle_y = np.angle(fft_y)  # 取复数的角度                     相位谱
    
44. # 我们发现，振幅谱的纵坐标很大，而且具有对称性
    
45. # 关键：关于振幅值很大的解释以及解决办法——归一化和取一半处理
    
46. normalization_y = abs_y / N  # 归一化处理（双边频谱）
    
47. normalization_half_y = normalization_y[range(int(N / 2))]  # 由于对称性，只取一半区间（单边频谱）
    
48. 
    
49. plt.figure(figsize=(16, 8))
    
50. 
    
51. plt.subplot(231)
    
52. plt.plot(x_time, y)
    
53. plt.title('原始波形')
    
54. plt.xlabel('time')
    
55. # plt.plot(x_time[:50], y[:50])
    
56. # plt.title('原始波形(前50组样本)')
    
57. 
    
58. plt.subplot(232)
    
59. plt.plot(x, fft_y, 'black')
    
60. plt.title('双边振幅谱(未求振幅绝对值)', fontsize=9, color='black')
    
61. 
    
62. plt.subplot(233)
    
63. plt.plot(x, abs_y, 'r')
    
64. plt.title('双边振幅谱(未归一化)', fontsize=9, color='red')
    
65. 
    
66. plt.subplot(234)
    
67. plt.plot(x, angle_y, 'violet')
    
68. plt.title('双边相位谱(未归一化)', fontsize=9, color='violet')
    
69. 
    
70. plt.subplot(235)
    
71. plt.plot(x, normalization_y, 'g')
    
72. plt.title('双边振幅谱(归一化)', fontsize=9, color='green')
    
73. 
    
74. plt.subplot(236)
    
75. plt.plot(half_x, normalization_half_y, 'blue')
    
76. plt.title('单边振幅谱(归一化)', fontsize=9, color='blue')
    
77. 
    
78. # plt.savefig('fft.png', dpi=200)
    
79. plt.show()
    

复制代码

东方耀

numpy中也有fft的实现：
from numpy import fft
fft_y = fft.fft(y)

东方耀

东方耀 发表于 2021-1-26 08:51
numpy中也有fft的实现：
from numpy import fft
fft_y = fft.fft(y)

1. import numpy as np
   
2. import pandas as pd
   
3. import math
   
4. import matplotlib.pyplot as plt
   
5. 
   
6. # 基于python(numpy)的快速傅里叶变换FFT
   
7. x = np.linspace(0, 2 * np.pi, 50)
   
8. # print(x)
   
9. wave = np.cos(x)
   
10. transformed = np.fft.fft(wave)  # 傅里叶变换
    
11. 
    
12. # plt.plot(transformed) #绘制变换后的信号
    
13. 
    
14. plt.plot(np.fft.ifft(transformed))  # 反变换
    
15. 
    
16. plt.show()

复制代码