import numpy as np
import math
from matplotlib.pyplot import *
import scipy.signal
b=[1]
a=[1,-1]
[w,h] = scipy.signal.freqz(b,a)
figure(figsize=(6,5))
plot(np.log10(w/(2*math.pi)),20*np.log10(np.abs(h)),'r')
xlabel('log(f/fe)')
ylabel('GdB')
grid()
                

figure(figsize=(6,5))
plot(np.log10(w/(2*math.pi)),np.angle(h),'r')
xlabel('log(f/fe)')
ylabel('phi')
grid()
                

b=[0.5,0.5]
a=[1,-1]
[w,h] = scipy.signal.freqz(b,a)
figure(figsize=(6,5))
plot(np.log10(w/(2*math.pi)),20*np.log10(np.abs(h)),'r')
xlabel('log(f/fe)')
ylabel('GdB')
grid()
                

figure(figsize=(6,5))
plot(np.log10(w/(2*math.pi)),np.angle(h),'r')
xlabel('log(f/fe)')
ylabel('phi')
grid()
                

figure(figsize=(6,5))
plot(w/(2*math.pi),20*np.log10(np.abs(h)),'r')
xlabel('f/fe')
ylabel('GdB')
grid()
                

[r,p,k] = scipy.signal.residuez(b,a)
                

def signal(t):
    return np.cos(2*math.pi*t)\
    +0.5*np.cos(2*2*math.pi*t-math.pi/3)\
    +0.2*np.cos(4*2*math.pi*t-math.pi/4)
                

fe = 20.0
te=1.0/fe
t = np.arange(start=0.0,stop=3.0,step=te)
x = signal(t)
figure(figsize=(12,5))
plot(t,x,".-")
xlabel('t')
ylabel('x')
grid()
                

a=[1.0,-1.0]
b=[0.5*te,0.5*te]
y = scipy.signal.lfilter(b,a,x)
figure(figsize=(12,5))
plot(t,y,".-")
xlabel('t')
ylabel('y')
grid()
                

def integre(t):
    return np.sin(2*math.pi*t)/(2*math.pi)\
    +0.5*(np.sin(2*2*math.pi*t-math.pi/3)-np.sin(-math.pi/3))/(2*2*math.pi)\
    +0.2*(np.sin(4*2*math.pi*t-math.pi/4)-np.sin(-math.pi/4))/(4*2*math.pi)
xi = integre(t)
figure(figsize=(12,5))
plot(t,xi,".-")
xlabel('t')
ylabel('xi')
grid()
                

zi = scipy.signal.lfiltic(b,a,y=[integre(-te)],x=[signal(-te)])
[y,zf] = scipy.signal.lfilter(b,a,x,axis=-1,zi=zi)
figure(figsize=(12,5))
plot(t,y,".-")
xlabel('t')
ylabel('y')
grid()
                

a=[1]
b=[1,-1]
[w,h] = scipy.signal.freqz(b,a)
figure(figsize=(6,5))
plot(np.log10(w/(2*math.pi)),20*np.log10(np.abs(h)),'r')
xlabel('log(f/fe)')
ylabel('GdB')
grid()
                

figure(figsize=(6,5))
plot(np.log10(w/(2*math.pi)),np.angle(h)/np.pi,'r')
xlabel('log(f/fe)')
ylabel(r'$\phi/\pi$')
grid()
                

fe = 100.0
te=1.0/fe
t = np.arange(start=0.0,stop=2.0,step=te)
x = signal(t)
def derive(t):
    return -np.sin(2*math.pi*t)*2*math.pi\
    -0.5*np.sin(2*2*math.pi*t-math.pi/3)*2*2*math.pi\
    -0.2*np.sin(4*2*math.pi*t-math.pi/4)*4*2*math.pi
a=[te]
b=[1,-1]
y = scipy.signal.lfilter(b,a,x)
figure(figsize=(12,5))
plot(t,y,".-")
xlabel('t')
ylabel('y')
axis([0,2,-20,20])
grid()                   
                    


xd = derive(t)
figure(figsize=(12,5))
plot(t,xd,".-")
xlabel('t')
ylabel('xd')
axis([0,2,-20,20])
grid() 
                    

import random
def signal(t):
    return np.cos(2*math.pi*t)\
    +0.5*np.cos(2*2*math.pi*t-math.pi/3)\
    +0.2*np.cos(4*2*math.pi*t-math.pi/4)\
    +0.05*random.uniform(-1.0,1.0)
fe = 100.0
te=1.0/fe
t = np.arange(start=0.0,stop=2.0,step=te)
n = t.size
x = np.zeros(n)
for k in range(n):
    x[k] = signal(te*k)
figure(figsize=(12,5))
plot(t,x,".-")
xlabel('t')
ylabel('x')
grid()
                    

y = scipy.signal.lfilter(b,a,x)
figure(figsize=(12,5))
plot(t,y,".-")
xlabel('t')
ylabel('y')
axis([0,2,-20,20])
grid() 
                    

P=3
h = scipy.signal.firwin(numtaps=2*P+1,cutoff=[0.1],nyq=0.5,window='hann')
b_pb=h
a_pb=[1.0]
[w,hf] = scipy.signal.freqz(b_pb,a_pb)
figure(figsize=(6,5))
plot(w/(2*math.pi),20*np.log10(np.abs(hf)),'r')
xlabel('f/fe')
ylabel('GdB')
grid()
                    

x1 = scipy.signal.convolve(x,h,mode='valid')
n1 = x1.size
t1 = (np.arange(n1)+P)*te
y = scipy.signal.lfilter(b,a,x1)
figure(figsize=(12,5))
plot(t1,y,".-")
xlabel('t')
ylabel('y')
axis([0,2,-20,20])
grid()
                    

hc = np.zeros(h.size)
for k in range(1,h.size):
    hc[k] = (h[k]-h[k-1])/te
y = scipy.signal.convolve(x,hc,mode='valid')
figure(figsize=(12,5))
plot(t1,y,".-")
xlabel('t')
ylabel('y')
axis([0,2,-20,20])
grid()
                    

fe = 1000.0
te=1.0/fe
t = np.arange(start=0.0,stop=2.0,step=te)
n = t.size
x = np.zeros(n)
for k in range(n):
    x[k] = signal(te*k)
figure(figsize=(12,5))
plot(t,x,".-")
xlabel('t')
ylabel('x')
grid()
                     

P=30
h = scipy.signal.firwin(numtaps=2*P+1,cutoff=[0.01],nyq=0.5,window='hann')
b_pb=h
a_pb=[1.0]
[w,hf] = scipy.signal.freqz(b_pb,a_pb)
figure(figsize=(6,5))
plot(w/(2*math.pi),20*np.log10(np.abs(hf)),'r')
xlabel('f/fe')
ylabel('GdB')
grid()
                      

x1 = scipy.signal.convolve(x,h,mode='valid')
n1 = x1.size
t1 = (np.arange(n1)+P)*te
figure(figsize=(12,5))
plot(t1,x1,".-")
xlabel('t')
ylabel('y')
grid()
                      

x2 = x1[0::10]
t2 = t1[0::10]
figure(figsize=(12,5))
plot(t2,x2,".-")
xlabel('t')
ylabel('y')
grid()
                      

y2 = scipy.signal.lfilter(b,a,x2)
figure(figsize=(12,5))
plot(t2,y2,".-")
xlabel('t')
ylabel('y')
axis([0,2,-20,20])
grid()