import numpy as np
import matplotlib.pyplot as plt

SincPoints=101
GraficoCol=2
GraficoRow=3
samples=1024
freq=1

plt.figure(1)

# kaiser window
window = np.kaiser(SincPoints,8.6)
plt.subplot(GraficoRow,GraficoCol,1)
plt.title('Kaiser window')
plt.plot(window)

# sinc function
x = np.linspace(-4, 4, SincPoints)
y=np.sinc(x)
plt.subplot(GraficoRow,GraficoCol,2)
plt.title('Sinc function')
plt.plot(x, y)

# windowed sinc function
yw=y*window
print (yw)
plt.subplot(GraficoRow,GraficoCol,3)
plt.title('Windowed Sinc function')
plt.plot(x, yw)


# sin function
xs=np.arange(0.0,100,0.1)
ys=np.sin(xs*0.2)
plt.subplot(GraficoRow,GraficoCol,4)
plt.title('Example Sin function')
plt.plot(xs, ys)


# convolute sin with sinc
print ("yw")
print (len(yw))
print ("sin")
print (len(ys))
c=np.zeros(len(ys),dtype=float)

for i in range(len(ys)):
    for i1 in range(len(yw)):
        ix=i+i1-int(len(yw)/2)
        if ((ix>0)and(ix<len(ys))):
            c[ix]+=ys[i]*yw[i1]
           
for i in range(int(len(yw)/2)):
    c[i]=0
    c[len(c)-i-1]=0
massimo=np.max(c)
minimo=np.min(c)
            
print ("massimo")
print (massimo)
print ("minimo")
print (minimo)

plt.subplot(GraficoRow,GraficoCol,5)
plt.title('Example Sin convoluted with Sinc')
plt.plot(xs, c)    
   
          
# filtro
# frequenza della sinusoide per il test
freq=np.arange(0.1,10,0.1)
FilterResponse=np.zeros(len(freq),dtype=float)

for f in range(len(FilterResponse)):
    
    # sinusoide 
    xsin=np.arange(0.0,100,0.1)
    ysin=np.sin(xs*freq[f])
    print(freq[f])
    # colvoluzione della sinc windowed function
    c=np.zeros(len(ysin),dtype=float)
    for i in range(len(ysin)):
        for i1 in range(len(yw)):
            ix=i+i1-int(len(yw)/2)
            if ((ix>0)and(ix<len(ysin))):
                c[ix]+=ysin[i]*yw[i1]
              
    for i in range(int(len(yw)/2)):
        c[i]=0
        c[len(c)-i-1]=0
    massimo=np.max(c)
    minimo=-np.min(c)
    if (minimo>massimo):
        FilterResponse[f]=minimo
    else:
        FilterResponse[f]=massimo
        
    
plt.subplot(GraficoRow,GraficoCol,6)
plt.title('Filter response')
plt.plot(freq, FilterResponse)     
    



plt.show()