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
65
66
67
68
69
70
71
72
73
74
75
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack
n_Sample = 1000
f_Signal = 50
f_sampling = f_Signal*n_Sample
def get_rms(signal):
return abs(np.sqrt(sum(signal**2)/len(signal)))
def get_max(signal):
return abs(max(signal))
def get_deal_sin(u_RMS=50):
x_signal = np.linspace(0.0, f_Signal ** -1, n_Sample)
y_signal = u_RMS * np.sqrt(2) * np.sin(f_Signal * 2.0 * np.pi * x_signal)
y_fft = fftpack.fft(y_signal) * 2 / n_Sample
x_fft = [x_ * (f_sampling / len(y_fft)) for x_ in range(len(y_fft) // 2)]
return ((x_signal, y_signal),(x_fft, y_fft))
def get_from_spectrum(u_RMS=50, harmonic_coefficient = 0.1, only_odd = False, title_show = True):
(x_signal, y_signal),(x_fft,y_fft) = get_deal_sin(u_RMS)
for n in list(range(2,50)):
if not(n%2) and only_odd:
continue
y_fft[n] = y_fft[1]*harmonic_coefficient**(n-1)
y_fft[n*-1] = y_fft[-1]*harmonic_coefficient**(n-1)
new_signal = fftpack.ifft(y_fft * n_Sample / 2).real
amp = get_max(new_signal)
rms = get_rms(new_signal)
crest_f = amp / rms
if crest_f < 3 and n < 49:
continue
if crest_f < 3:
return get_from_spectrum(u_RMS, harmonic_coefficient+0.01, only_odd)
title = 'Только нечетные гармоники' if only_odd else 'Все гармоники'
print('\n'+title+':') if title_show else None
if rms > 51.0 or rms < 49.0:
print('Текущее СКЗ {}'.format(rms))
print('Изменение СКЗ затравки синуса с {} до {}\n'.format(rms, 50*u_RMS/rms))
return get_from_spectrum(50*u_RMS/rms, harmonic_coefficient, only_odd, title_show = False)
print('Изменено {} гармоник'.format(n-1 if not only_odd else str((n-1)//2) + ' нечетных'))
print('Коэффициент гармоник: {}'.format(harmonic_coefficient))
print('Крест-фактор: {}'.format(crest_f))
print('СКЗ: {}'.format(rms))
print('Амплитуда (max): {}\n\n'.format(amp))
plt.subplot(2,1,1)
plt.plot(x_signal,new_signal)
plt.grid(True)
plt.title(title)
plt.subplot(2,1,2)
plt.stem(x_fft,abs(y_fft[:len(y_fft)//2]))
plt.axis([0, 50*50, 0, abs(y_fft[1])+1])
plt.grid(True)
plt.show()
break
if __name__ == '__main__':
get_from_spectrum(only_odd=False)
get_from_spectrum(only_odd=True)