%pylab inline
from __future__ import print_function
from __future__ import division

linspace(0, 1, 10)

wt = linspace(0, 6 * pi, 20000)

oscillation = cos(wt)
plot(oscillation);

f = 3
w = 2 * pi * f

phase_t = linspace(0, w, 50000)
impulse = cos(phase_t) + cos(2*phase_t)
impulse /= 2
plot(impulse);

impulse = cos(phase_t) + cos(2*phase_t) + cos(3*phase_t) + cos(4*phase_t)
impulse /= 4
plot(impulse);

N = 100
phase = linspace(0, 6*pi, 50000)
harmonics = arange(N) + 1
impulse = zeros_like(phase)
for harmonic in harmonics:
    impulse += cos(harmonic * phase)
    
impulse /= N
plot(impulse);

plot(impulse)
xlim((0, 5000))

N = 100
phase = linspace(0, 6*pi, 50000)
harmonics = arange(N) + 1
impulse = zeros_like(phase)
for harmonic in harmonics:
    impulse += sin(harmonic * phase)
    
impulse /= N
plot(impulse);

comb = [1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] * 25
comb [-1] = 1 # unholy trick to make some things simpler later...
# plot(comb)
stem(comb)
ylim(0, 1.1)
xticks(())
title("Dirac Comb")

from scipy.special import jn
jn(0, 0), jn(0, 1)

jn(1,0), jn(7,0.3), jn(1, 1)

x = linspace(0,10, 500)
plot(x, jn(1,x))
title('Bessel function of the first type, order 1');

x = linspace(0,10, 500)
plot(x, jn(1,x))
plot(x, jn(1,x) * comb);

plot(x, jn(1,x) * comb);

sampled = jn(1,x) * comb
samples = []
for i in range(len(comb)):
    if comb[i] == 1:
        samples.append(sampled[i])

plot(samples, 'o');

len(samples)

x = linspace(0,10, 500)
x_sampled = linspace(0,10, 26)
plot(x_sampled, samples)
plot(x, jn(1,x));

plot(x_sampled, samples, 'x-')
plot(x, jn(1,x));
xlim((0.5, 3))
ylim((0.3, 0.65))

from scipy.interpolate import interp1d
interpf = interp1d(linspace(0,10, len(samples)), samples)
sqe = (jn(1,x) - interpf(x))**2
plot(interpf(x))
plot(jn(1,x))

twinx()
plot(sqe, 'r')
axis(ymax= sqe.max())
ylabel('squared error', color='r',fontsize=18);

sqe = (jn(1,x) - interpf(x))**2
plot(interpf(x))
plot(jn(1,x))
ylim((0.3, 0.65))

twinx()
axis(ymax= sqe.max())
gca().set_ylabel('squared error', color='r',fontsize=18)
plot(sqe, 'r')
xlim((50,150));

MSE = sqe.mean()
print(MSE)

x = linspace(0, 2*pi, 10)
x0 = linspace(0, 2*pi, 100)
plot(x, sin(x))
plot(x0, sin(x0))

phs = linspace(0, 10 * 2 * pi, 300)
plot(sin(phs))

phs = linspace(0, 10 * 2 * pi, 300)
plot(phs, sin(phs))
phs = linspace(0, 10 * 2 * pi, 12)
plot(phs, sin(phs), 'o--')

phs = linspace(0, 10 * 2 * pi, 300)
plot(phs, sin(phs))
phs = linspace(0, 10 * 2 * pi, 21)
plot(phs, sin(phs), 'o--')

phs = linspace(0, 10 * 2 * pi, 300)
plot(phs, cos(phs))
phs = linspace(0, 10 * 2 * pi, 11)
plot(phs, cos(phs), 'o--')

x = linspace(0, 2*pi, 300)
f = sin(x)
plot(f);

f = sin(3 * x)
f2 = (f*2).astype(int)
plot(f) 
plot(f2); 

f = sin(x)
f2 = (f*2).astype(int)
f4 = (f*4).astype(int)
plot(f)
plot(f2)
plot(f4/3.0)
legend(['sine', '2-bit', '3-bit']);

2**2, 2**3

2**16

#integer representations
x = linspace(0, 2*pi, 100000)
f = sin(x)
N = 16 # number of bits
max_value = 2**(N-1) - 1
f16 = (f*(max_value)).astype(int16)
plot(f16);

plot(f16, 'x-' )
xlim((0, 50))
ylim((0, 80))

N = 5 # number of bits
max_value = 2**(N-1) - 1
f8 = (f*(max_value)).astype(int8)
plot(f8)

plot(f8, 'o')
xlim((0,20))
ylim((0, 10))
grid();

2**24

N = 16
20 * log10((2 ** (N - 1))/1)

def dynrange(N):
    return 20 * log10((2 ** (N - 1))/1)

print(dynrange(8), dynrange(16), dynrange(24))

20*log10(0.5)

20*log10(60/30)

22050 * 16

352800 * 60

21168000/(8 * 1024)

from scipy.io import wavfile

sr,audio = wavfile.read('passport.wav')
plot(audio)
print(audio.max(), audio.min(), sr)

2**16

audio.dtype

x = linspace(0, 2*pi, 50)
plot(sin(x), 'o')
plot(diff(sin(x)))

dpcm = diff(audio)
plot(dpcm)
dpcm.max(), dpcm.min()

log(max(dpcm.max(), abs(dpcm.min())))/log(2)

log(max(audio.max(), abs(audio.min())))/log(2)

x = linspace(-1,1, 100)
ylabel('out')
xlabel('in')
plot(x, x)
grid()

x = linspace(-1,1, 100)
bits = 8
mu = (2**bits) - 1

mu_shaping = sign(x)*log(1 + mu *abs(x))/log(1+mu)
plot(x, x, label= 'linear')
plot(x, mu_shaping, label= 'mu-law out')

legend(loc='best')
grid()

x = linspace(0, 2*pi, 100)
y = sin(x)

def mu_law(insig, nbits = 8):
    mu = (2**nbits) - 1
    return sign(insig)*log(1 + mu *abs(insig))/log(1+mu)

plot(x, y)
plot(x, mu_law(y))

x = linspace(-1,1, 100)
A = 87.5
x = linspace(-1,1, 100)
a_shaping = sign(x) * where(abs(x) < 1/A, A*abs(x)/(1+log(A)), (1 + log(A*abs(x)))/(1+log(A)))
plot(x, x,label = 'in')
plot(x, mu_shaping, label= 'mu-law out')
plot(x, a_shaping, label = 'A-law out')


legend(loc='best')

plot(x, label = 'in')
plot(x, mu_shaping, label= 'mu-law out')
plot(x, a_shaping, label = 'A-law out')

legend(loc='best')

xlim((0,0.5))
ylim((0,1))
     

samplehold = interp1d(linspace(0,10, len(samples)), samples, kind='zero')
new_x = linspace(0, 10, 500)
plot(new_x, samplehold(new_x))

from scipy.signal import butter, lfilter

b, a = butter(4, 0.2, 'low')
lopassed = lfilter(b, a, samplehold(new_x))
plot(lopassed)



b, a = butter(4, 0.1, 'low')
lopassed = lfilter(b, a, samplehold(new_x))
plot(lopassed)

b, a = butter(4, 0.05, 'low')
lopassed = lfilter(b, a, samplehold(new_x))
plot(lopassed, lw='4')

from scipy.signal import bessel

b, a = bessel(10, 0.05, 'low')
lopassed = lfilter(b, a, samplehold(new_x))
plot(lopassed)


plot(new_x*50, jn(1,new_x))
plot(lopassed)

plot(new_x*50, jn(1,new_x))
plot(lopassed[58:])