DFT处理合体正弦波

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
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,ifft

# 样本采样率（每秒钟采几个样）
sr=100
# 样本采样时间间隔
T_s=1/sr
# 样本采样的时间长度为5s，样本采样的个数为500，最小频率=1/最大周期=1/5s=0.2Hz
t=np.arange(0,5,T_s)
y_0=3*np.sin(2*np.pi*t)
y_1=np.sin(2*np.pi*4.3*t)
y_2=0.5*np.sin(2*np.pi*7*t)
y_3=y_0+y_1+y_2

# fig,ax=plt.subplots(2,2,figsize=(6,6),sharex='col')
# fig.subplots_adjust(hspace=0.2,wspace=0.3)

# ax[0,0].plot(t,y_0)
# ax[0,1].plot(t,y_1)
# ax[1,0].plot(t,y_2)
# ax[1,1].plot(t,y_3)

# ax[1,0].set_xlabel('Time')
# ax[1,1].set_xlabel('Time')
# ax[0,0].set_ylabel('Amplitude')
# ax[1,0].set_ylabel('Amplitude')

# ax[0,0].set_title('$y_0=3\sin(2\pi\cdot t)$')
# ax[0,1].set_title('$y_1=\sin(2\pi\cdot 4.3t)$')
# ax[1,0].set_title('$y_2=0.5\sin(2\pi\cdot 7t)$')
# ax[1,1].set_title('$y_3=y_0+y_1+y_2$')

# plt.show()

y_3_fft=fft(y_3)
N=len(t)
# 采样数据的idx
n=np.arange(N)
# 总的采样时间
T=N/sr
# 频度区间
freq=n/T
# 实际幅度
y_3_fft_norm=y_3_fft/N*2 # 还有一个问题是左图中虽然有明显的三个振幅，但是这三个振幅对应的值却与原来函数 y_0,y_1,y_2不对应，这是因为离散傅立叶内部公式实现上的原因导致，细节不用纠结，记住这一步就行了。除以 N 是因为 scipy 包中封装的离散傅立叶变换公式为了和傅立叶变换公式保持一致，所以内部没有除以 N；乘以 2 是因为由于复数的引入，同一个振幅被分配至两个共轭复数上。

fig,ax=plt.subplots(1,2,figsize=(10,4),dpi=150)
ax[0].stem(freq,np.abs(y_3_fft),'b',markerfmt=' ',basefmt='-b')
ax[0].set_xlabel('Freq(Hz)')
ax[0].set_ylabel('FFT Amplitude |X(freq)|')
ax[1].stem(freq,np.abs(y_3_fft_norm),'b',markerfmt=' ',basefmt='-b')
ax[1].set_xlim(0,10)
ax[1].set_xlabel('Freq(Hz)')
ax[1].set_ylabel('FFT Amplitude |X(freq)|')
ax[1].annotate(r'$y_0=3\sin(2\pi\cdot t)$',xy=(1,2.8),xytext=(+30,-30),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=.2'))
ax[1].annotate(r'$y_1=\sin(2\pi\cdot 4.3t)$',xy=(4.2,.4),xytext=(+10,50),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
ax[1].annotate(r'$y_2=0.5\sin(2\pi\cdot 7t)$',xy=(7,.4),xytext=(+10,30),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
plt.show()

应用——信号降噪

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
65
66
67
68
69
70
71
72
73
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,ifft
import pandas as pd

# 采样频率
sr=1000
# 采样时间间隔
ts=1/sr
# 样本采样点
t=np.arange(0,2,ts)
# 原始信号
f_clean=5
freq=10
f_clean+=2*np.sin(2*np.pi*freq*t+3)
freq=30
f_clean+=5*np.sin(2*np.pi*freq*t)
# 信号噪音
f_noise=f_clean+3*np.random.randn(len(t))
# 绘制信号
# fig,ax=plt.subplots(figsize=(12,3))
# ax.plot(t,f_noise,linewidth=.5,color='c',label='Noisy')
# ax.plot(t,f_clean,linewidth=.5,color='r',label='Clean')
# ax.set_xlabel('Sampling Time')
# ax.set_ylabel('Amplitude')
# ax.legend()
# plt.show()

X=fft(f_noise)
N=len(X)
# 采集数据的idx
n=np.arange(N)
# 总的采样时间
T=N/sr
freq=n/T

# fig,ax=plt.subplots(2,1,figsize=(10,8))
# ax[0].stem(freq,np.abs(X),'b',markerfmt=' ',basefmt='-b')
# ax[0].set_xlabel('Freq(Hz)')
# ax[0].set_ylabel('FFT Amplitude |X(freq)|')
# # 实际幅度
# X_norm=X/N*2
# ax[1].stem(freq,np.abs(X_norm),'b',markerfmt=' ',basefmt='-b')
# ax[1].set_xlim(-1,40)
# ax[1].set_xlabel('Freq(Hz)')
# ax[1].set_ylabel('FFT Amplitude |X(freq)|')
# ax[1].annotate(r'5(Constant terms)',xy=(0,5),xytext=(+10,50),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
# ax[1].annotate(r'$2\sin(2\pi\cdot 10t+3)$',xy=(10,2),xytext=(+10,30),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
# ax[1].annotate(r'$5\sin(2\pi\cdot 30t)$',xy=(30,4),xytext=(+10,30),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
# plt.show()

# 对DFT结果进行噪音去除，过滤除去低振幅项，并可视化
freq_clean=pd.Series(X).apply(lambda x: x if np.abs(x)>1000 else 0).to_numpy()
# fig,ax=plt.subplots(figsize=(10,4))
# ax.stem(freq,np.abs(freq_clean),'b',markerfmt=' ',basefmt='-b')
# ax.set_xlim(-1,100)
# ax.set_xlabel('Freq(Hz)')
# ax.set_ylabel('FFT Amplitude |X(freq)|')
# ax.annotate(r'5(Constant terms)',xy=(0,5*N/2),xytext=(+10,50),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
# ax.annotate(r'$2\sin(2\pi\cdot 10t+3)$',xy=(10,2*N/2),xytext=(+10,30),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
# ax.annotate(r'$5\sin(2\pi\cdot 30t)$',xy=(30,4*N/2),xytext=(+10,30),textcoords='offset points',arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=-.5'))
# plt.show()

# 将过滤后的信号进行傅里叶逆变换并可视化
ix=ifft(freq_clean)
# 绘制信号
fig,ax=plt.subplots(figsize=(12,3))
ax.plot(t,ix,linewidth=.5,color='c',label='IFFT')
ax.plot(t,f_clean,linewidth=.5,color='r',linestyle='-.',label='Raw signal',alpha=0.7)
ax.set_xlabel('Sampling Time')
ax.set_ylabel('Amplitude')
ax.legend()
plt.show()

实现DFT

X_k=\sum_{n=0}^{N-1}x_n\cdot e^{-i2\pi kn/N}是DFT实现的算法

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
from scipy.fft import fft

def DFT_slow(x):
	x=np.asarray(x,dtype=float) # ensure the data type
	N=x.shape[0] # get the x array length
	n=np.arange(N) # 1d array
	k=n.reshape((N,1)) # 2d array,10x1,aka column array
	M=np.exp(-2j*np.pi*k*n/N)
	return np.dot(M,x) # [a,b]\cdot [c,d]=ac+bd, it is a sum

x=np.random.random(1024)
print(np.allclose(DFT_slow(x),fft(x)))

https://blog.manyacan.com/archives/2041/ 离散傅立叶变换的Python实现

2311082236