%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division

sig = random.random(2048)*2 -1

Pxx, freqs, times, im = specgram(sig, interpolation='nearest');
colorbar()

Pxx.shape

2048/256

plot(Pxx[:,10])

filtered = (sig + r_[0,sig[:-1]])
Pxx2, freqs, times, im = specgram(filtered);
colorbar()

plot(Pxx2[:,10])

plot(Pxx[:,10])
plot(Pxx2[:,10])

plot(Pxx[:,10])
plot(Pxx2[:,10]/4.0)

from scipy.signal import freqz
frq, resp = freqz([1,1])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()


plot(frq, angle(resp))
title('Phase Response')
grid()

filtered2 = (sig + 0.5*r_[0,sig[:-1]])
Pxx3, freqs, times, im = specgram(filtered2);

plot(Pxx3[:,10])

plot(Pxx[:,10])
plot(Pxx2[:,10]/4.0)
plot(Pxx3[:,10]/(1.5**2))

from scipy.signal import freqz
frq, resp = freqz([1,1])
semilogy(frq, abs(resp))

frq, resp = freqz([1,0.5])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()

frq, resp = freqz([1,1])
semilogy(frq, abs(resp))

frq, resp = freqz([1,2])
semilogy(frq, abs(resp))


frq, resp = freqz([1,0.5])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()

frq, resp = freqz([1,-1])
semilogy(frq, abs(resp))
title('Frequency Response')
grid()

from scipy.signal import lfilter
filtered4 = lfilter([1],[1, 1], sig)
specgram(filtered4);

frq, Y = freqz([1], [1,1])
semilogy(frq, abs(Y))


frq, Y = freqz([1], [1,0.5])
semilogy(frq, abs(Y))

frq, Y = freqz([1], [1,-0.5])
semilogy(frq, abs(Y))

frq, Y = freqz([1], [1,-0.5])
plot(frq, angle(Y))

from scipy.signal import tf2zpk

tf2zpk([1,1],[1])

tf2zpk([1],[1, 1])

tf2zpk([2, 2],[1])

def PoleZeroPlot(b, a):
    (zeros,poles,gain) = tf2zpk(b, a)
    angle = np.linspace(-np.pi,np.pi,50)
    cirx = np.sin(angle)
    ciry = np.cos(angle)
    figure()
    plot(poles.real, poles.imag, 'x', zeros.real, zeros.imag, 'o', cirx,ciry, 'k-')
    grid()
    
    xlim((-2, 2))
    xlabel('Real')
    ylim((-1.5, 1.5))
    ylabel('Imag')
    gcf().set_figwidth(5)
    return (zeros,poles,gain)

PoleZeroPlot([1,1],[1])

PoleZeroPlot([1],[1,1])

PoleZeroPlot([1],[1,1.1])

PoleZeroPlot([1],[1.1,1.1])

b = [1,0,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))

b = [1,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))

b = [1,0,0,0,0,0,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp))

twinx()

plot(frq, angle(resp), 'r')

b = [1,0,0,0,0,0,0,-1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))

b = [1,0,0,0,0,0,0,1]
a = [1]
frq, resp = freqz(b,a)
plot(frq, abs(resp))

b = [1,0,0,0,0,0,0,-1]
a = [1,0,0,0,0,0,0,1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp));

b = [1]
a = [1,0,0,0,0,0,0,1]
frq, resp = freqz(b,a, worN=8192)
semilogy(frq, abs(resp));

b = [1, 1, 1]
a = [1]
frq, resp = freqz(b,a)
semilogy(frq, abs(resp))

R = 0.8
theta_c = 1.0
b = [1, -2*R*cos(theta_c), R**2]
a = [1]

frq, resp = freqz(b,a)
semilogy(frq, abs(resp))

R = 1
theta_c = 2.0
b = [1, -2*R*cos(theta_c), R**2]
a = [1]

frq, resp = freqz(b,a)
plot(frq, abs(resp))

R = 1
theta_c = 1.5
b = [1, -2*R*cos(theta_c), R**2]
a = [1]

frq, resp = freqz(b,a)
plot(frq, abs(resp))

PoleZeroPlot(b,a)

R = 1
theta_c = 2.0
b = [1]
a = [1, -2*R*cos(2.0), R**2]

frq, resp = freqz(b,a)
plot(frq, abs(resp))

poles, zeros, k = PoleZeroPlot(b,a)

abs(zeros)

R = 0.9
theta_c = 1.5
b = [1]
a = [1, -2*R*cos(theta_c), R**2]

frq, resp = freqz(b,a)
plot(frq, abs(resp))

poles, zeros, k = PoleZeroPlot(b,a)

# low shelf-filter

Fs = 44100
f0 = 10000.0
dBgain = 30.0
S = 1.0 # shelf slope
# -----------------------
A  = 10**(dBgain/40)

w0 = 2*pi*f0/Fs
alpha = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 ) 
       
b0 =    A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
b1 =  2*A*( (A-1) - (A+1)*cos(w0)                   )
b2 =    A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
a0 =        (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
a1 =   -2*( (A-1) + (A+1)*cos(w0)                   )
a2 =        (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha

w, h = freqz([b0, b1, b2], [a0, a1, a2])
semilogy(w,abs(h))

Fs = 44100
f0 = 10000.0
dBgain = -20.0
S = 2.0 # shelf slope

A  = 10**(dBgain/40)

w0 = 2*pi*f0/Fs
alpha = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 ) 
    
    
b0 =    A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
b1 =  2*A*( (A-1) - (A+1)*cos(w0)                   )
b2 =    A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
a0 =        (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
a1 =   -2*( (A-1) + (A+1)*cos(w0)                   )
a2 =        (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha

w, h = freqz([b0, b1, b2], [a0, a1, a2])
plot(w,abs(h))

shelved = lfilter([b0, b1, b2], [a0, a1, a2], sig)
specgram(shelved);

PoleZeroPlot([b0, b1, b2], [a0, a1, a2]);

from scipy.signal import cheby1

b, a = cheby1(4, 0.5, 0.1)
b,a

w, h = freqz(b,a)
plot(w, np.abs(h), 'b')

PoleZeroPlot(b,a)

ripples = [2, 0.5, 0.1, 0.01]

for ripple in ripples:
    b, a = cheby1(4, ripple, 0.6)
    frq, resp = freqz(b,a)
    plot(frq, abs(resp))


legend(ripples)
title('Different ripple values for a Chebyshev I filter')
grid()

orders = [2,3,4,6]

for order in orders:
    b, a = cheby1(order, 2, 0.6)
    frq, resp = freqz(b,a)
    plot(frq, abs(resp))


legend(orders)
title('Different orders for a Chebyshev I filter')
grid()

from scipy.signal import cheby2

ripples = [12, 15, 20, 40]

for ripple in ripples:
    b, a = cheby2(4, ripple, 0.6)
    frq, resp = freqz(b,a)
    plot(frq, abs(resp))

legend(ripples)
title('Different ripple values for a Chebyshev II filter')
grid()

from scipy.signal import iirdesign

Wp = 0.5  # Cutoff frequency 
Ws = 0.6   # Stop frequency 
Rp = 0.1     # passband maximum loss (gpass)
As = 60      # stoppand min attenuation (gstop)
b,a = iirdesign(Wp, Ws, Rp, As, ftype='butter')
frq, resp = freqz(b,a)
plot(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')

title('Butterworth filter')
grid()

Wp = 0.5  # Cutoff frequency 
Ws = 0.6   # Stop frequency 
Rp = 1     # passband maximum loss (gpass)
As = 20      # stoppand min attenuation (gstop)

types = ['butter', 'ellip', 'cheby1', 'cheby2']

for t in types:
    b,a = iirdesign(Wp, Ws, Rp, As, ftype=t)
    frq, resp = freqz(b,a)
    plot(frq, abs(resp))
    
legend(types)
title('Filters')
grid()

from scipy.signal import firwin2

freqs = [0.0, 0.5, 1.0]
gains = [1.0, 1.0, 0.0]
order = 150

taps = firwin2(order, freqs, gains)

frq, resp = freqz(taps)
plot(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')


freqs = [0.0, 0.3, 0.5, 0.8, 1.0]
gains = [1.0, 0.2, 0.5, 0.1, 0.0]
order = 200

taps = firwin2(order, freqs, gains)

frq, resp = freqz(taps)
plot(frq, abs(resp))
twinx()
plot(frq, angle(resp), 'r')

shaped = lfilter(taps, [1], sig)
specgram(shaped);

from scipy.signal import remez
        
freqs = [0, 0.1, 0.2, 0.4, 0.45, 0.5]
gains = [0, 1, 0]

taps = remez(50, freqs, gains, type='bandpass')

frq, resp = freqz(taps)
plot(frq/(2*pi), abs(resp))
twinx()
plot(frq/(2*pi), angle(resp), 'r')

from scipy.signal import remez
        
freqs = [0, 0.1, 0.3, 0.4, 0.45, 0.5]
gains = [0, 1, 0]

taps = remez(16, freqs, gains, type='bandpass')

frq, resp = freqz(taps)
plot(frq/(2*pi), abs(resp))
twinx()
plot(frq/(2*pi), angle(resp), 'r')

from scipy.signal import freqz
g = 0.9
b = [-g, 0,0, 1]
a = [1, 0,0,-g]
frq, resp = freqz(b,a)
plot(frq, abs(resp))

plot(frq, angle(resp))
xlabel('freq')
ylabel('angle')

sig = sin(linspace(0, 2*pi*5, 1000, endpoint=False)) + sin(linspace(0, 2*pi*30, 1000, endpoint=False))
plot(sig)

from scipy.signal import lfilter
f = lfilter(b, a, sig)
plot(abs(rfft(f)))
plot(abs(rfft(sig)))
xlim((0,80))

plot(angle(rfft(f)))
plot(angle(rfft(sig)))

legend(['original', 'all-pass'])

stem(angle(rfft(f)),linefmt='b:')
stem(angle(rfft(sig)), 'g.-.')

legend(['original', 'all-pass'])
xlim(0,50)

angle(resp)[5], angle(resp)[30]

angle(rfft(f))[5] - angle(rfft(sig))[5]

angle(rfft(f))[30] - angle(rfft(sig))[30]

angle(rfft(f))[30] - angle(rfft(sig))[30] - (2 * pi)

plot(sig)
plot(f)
xlim((0, 200))
grid()