---

import numpy as np
from math import *
import matplotlib.pyplot as plt

def signal(x):
    return 1.0+sin(2.0*pi*x)+2.0*cos(4.0*pi*x)+0.5*cos(6.0*pi*x)

def filter(time, signal, fl, fh):
  n = len(signal)
  freq = np.fft.fftfreq(n, time[1] - time[0])
  spectr = np.fft.fft(signal)
  for i in range(n):
    if not fl < abs(freq[i]) < fh:
      spectr[i] = 0+0j
  return np.fft.ifft(spectr)

f=float(input('Опорная частота сигнала '))
T=float(input('Временной интервал '))
n=int(input('Число временных отчетов '))
t=np.linspace(0,T,n)
u=[0]*n
for i in range(n):
    u[i]=signal(f*t[i])

print('Расчёт БПФ...')
spec = np.fft.fft(u)
freq = np.fft.fftfreq(n, T / n)

print('Расчёт ФНЧ...')
sig_lf = filter(t, u, 0, 1.5*f)
print('Расчёт ПФ...')
sig_pb = filter(t, u, 1.5*f, 2.5*f)
print('Расчёт ФВЧ...')
sig_hf = filter(t, u, 2.5*f, 3.5*f)

print('Расчёт БПФ...')
spec_lf = np.fft.fft(sig_lf)
spec_pb = np.fft.fft(sig_pb)
spec_hf = np.fft.fft(sig_hf)
plt.plot(t, u)
plt.axis(xmin=0, xmax = 1/f)
plt.show()

plt.plot(freq[0:n//2], (np.hypot(spec.real, spec.imag)/n*2.0)[0:n//2])
plt.axis(xmin=0, xmax=4*f)
plt.show()

plt.plot(t, sig_lf.real, label='ФНЧ')
plt.plot(t, sig_pb.real, '--', label='ПФ')
plt.plot(t, sig_hf.real, '-.',  label='ФВЧ')
plt.legend(loc='best')
plt.axis(xmin=0, xmax=1/f)
plt.show()


plt.plot(freq[0:n//2], (np.hypot(spec_lf.real, spec_lf.imag)/n*2.0)[0:n//2], '-', label='ФНЧ')
plt.plot(freq[0:n//2], (np.hypot(spec_pb.real, spec_pb.imag)/n*2.0)[0:n//2], '--', label='ПФ')
plt.plot(freq[0:n//2], (np.hypot(spec_hf.real, spec_hf.imag)/n*2.0)[0:n//2], '-.',  label='ФВЧ')
plt.legend(loc='best')
plt.axis(xmin=0, xmax=4*f)
plt.show()