# 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(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) ```