import numpy as np import matplotlib pyplot as plt from scipy import f

 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 matplotlib.pyplot as plt
from scipy import fftpack, signal
def get_fft(signal_):
fft = fftpack.fft(signal_)
return abs(fft)*2/n_Sample
def get_fft_with_n_h(n, fft):
new_fft = [0, ] * len(fft)
for i in range(n):
new_fft[i] = fft[i]
new_fft[-i] = fft[-i]
return new_fft
def get_rms_with_fft(fft):
return np.sqrt(sum([x**2 for x in fft]))/2
def get_n_for_rms(fft, accuracy):
sqrt_ = np.sqrt(sum([x**2 for x in fft]))/2
for i in range(len(fft)):
new_fft = get_fft_with_n_h(i, fft)
if get_rms_with_fft(new_fft)>= (1-accuracy)*sqrt_:
return i
n_Sample = 100
f_Signal = 50
f_sampling = f_Signal*n_Sample
plt.subplot(211)
t = np.linspace(0, 1/f_Signal, n_Sample)
square_100 = signal.square(2*np.pi*50*t)
fft_100 = get_fft(square_100)
rms_graph_100 = [get_rms_with_fft(get_fft_with_n_h(i, fft_100)) for i in range(n_Sample)]
plt.plot(rms_graph_100),plt.grid()
plt.axhline((1-0.01)*get_rms_with_fft(fft_100),color='r')
plt.axvline(get_n_for_rms(fft_100, 0.01), color='r')
plt.axhline((1-0.001)*get_rms_with_fft(fft_100),color='g')
plt.axvline(get_n_for_rms(fft_100, 0.001), color='g')
plt.xlim(0,100),plt.ylim(0.8,1.1),plt.title('Зависимость СКЗ от гармоник, 100 гармоник')
plt.subplot(212)
n_Sample *= 10
f_sampling *= 10
t = np.linspace(0, 1/f_Signal, n_Sample)
square_1000 = signal.square(2*np.pi*50*t)
fft_1000 = get_fft(square_1000)
rms_graph_1000 = [get_rms_with_fft(get_fft_with_n_h(i, fft_1000)) for i in range(n_Sample)]
plt.plot(rms_graph_1000),plt.grid()
plt.axhline((1-0.01)*get_rms_with_fft(fft_1000),color='r')
plt.axvline(get_n_for_rms(fft_1000, 0.01), color='r')
plt.axhline((1-0.001)*get_rms_with_fft(fft_1000),color='g')
plt.axvline(get_n_for_rms(fft_1000, 0.001), color='g')
plt.xlim(0,1000),plt.ylim(0.8,1.1),plt.title('Зависимость СКЗ от гармоник, 1000 гармоник')
plt.show()