离散傅里叶变换
离散傅里叶变换及逆变换
这里直接上离散傅里叶变换的数学公式:
用D(x)代表傅里叶变换输出，E(x)代表原始波形数据
用ID(x)代表傅里叶逆变换输出
(x=0,1,2,3…)
用N代表数据量，也就是采样点数量
$D(k) = \sum_{n=0}^{N-1}E(n)(cos\frac{2\pi kn}{N} - jsin\frac{2\pi kn}{N})$
$ID(k) = \frac{1}{N}\sum_{n=0}^{N-1}D(n)(cos\frac{2\pi kn}{N} + jsin\frac{2\pi kn}{N})$
其中的”j”代表虚数单位
python代码如下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
| import math
pi = 3.1415926
def DFT(origin):
N= len(origin)
output = []
for i in range(N):
temp = 0
for j in range(N):
temp += origin[j] * (
(math.cos(2 * math.pi * i * j / N) - math.sin(2 * math.pi * i * j / N) * 1j)
)
output.append(temp)
return output
def IDFT(origin):
N= len(origin)
output = []
for i in range(N):
temp = 0
for j in range(N):
temp += origin[j] * (
(math.cos(2 * math.pi * i * j / N) + math.sin(2 * math.pi * i * j / N) * 1j)
)
output.append(temp / N)
return output
|
样例：
滤波
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
| import numpy as np
import math
import matplotlib.pyplot as plt
start_X = -2*np.pi
end_x = 2*np.pi
slice_num = 200
X = np.linspace(start_X,end_x,slice_num,endpoint=True)
wave1 = 10*np.sin(X) + 5*np.append(np.random.rand(len(X)-100),np.zeros(100))
wave1 = X + wave1 * 1j
wave1_real=np.real(wave1)
wave1_imag=np.imag(wave1)
plt.title('original data')
plt.plot(np.real(wave1), np.imag(wave1))
plt.show()
pi = 3.1415926
def DFT(origin):
N= len(origin)
output = []
for i in range(N):
temp = 0
for j in range(N):
temp += origin[j] * (
(math.cos(2 * math.pi * i * j / N) - math.sin(2 * math.pi * i * j / N) * 1j)
)
output.append(temp)
return output
def IDFT(origin):
N= len(origin)
output = []
for i in range(N):
temp = 0
for j in range(N):
temp += origin[j] * (
(math.cos(2 * math.pi * i * j / N) + math.sin(2 * math.pi * i * j / N) * 1j)
)
output.append(temp / N)
return output
dft_wave = DFT(wave1_imag)
idft_wave = IDFT(dft_wave)
dft_new_wave = dft_wave.copy()
for i in range(len(dft_new_wave)):
if(dft_new_wave[i].__abs__() <= 200):
dft_new_wave[i]=0
idft_wave = IDFT(dft_new_wave)
fig = plt.figure()
plt.title('DFT filter data')
plt.plot(np.real(wave1), np.real(idft_wave))
plt.show()
|
结果:
