%pylab inline
%pylab --no-import-all
pylab.rcParams['figure.figsize'] = 14, 4 

from IPython.display import Image
from scipy.io import wavfile
import Audio

display(Image(filename='data/beatles.png', width=400))

# fs will be the global "clock" of all discrete-time systems in this notebook
fs, data = wavfile.read("data/iff.wav")
# bring the 16bit wav samples into the [-1, 1] range
data = data / 32767.0
Audio.Audio(data=data, rate=fs, embed=True)

display(Image(filename='data/bd.png', width=600))

plot(data)
xlabel("sample")
ylabel("amplitude")

s = abs(np.fft.fftpack.fft(data[10000:40000]));
s = s[0:len(s)/2]
plot(np.linspace(0,1,len(s))*(fs/2), s)
xlabel("frequency (Hz)")
ylabel("magnitude")

from scipy import signal

w, h = signal.freqz(1, [1, 0, 0, 0, 0, 0, 0, 0, 0, -.99**9])
plot(w, abs(h))

class guitar:
    def __init__(self, pitch=110, fs=24000):
        # init the class with desired pitch and underlying sampling frequency
        self.M = np.round(fs / pitch) # fundamental period in samples
        self.R = 0.9999               # decay factor
        self.RM = self.R ** self.M 
        self.ybuf = np.zeros(self.M)  # output buffer (circular)
        self.iy = 0                   # index into out buf
        self.xbuf = 0                 # input buffer (just one sample)
        
    def play(self, x):
        y = np.zeros(len(x))
        for n in range(len(x)):
            t = x[n] - self.R * self.xbuf + self.RM * self.ybuf[self.iy]
            self.ybuf[self.iy] = t
            self.iy = (self.iy + 1) % self.M
            self.xbuf = x[n]
            y[n] = t
        return y

# create a 2-second signal
d = np.zeros(fs*2)
# impulse in zero (string plucked)
d[0] = 1

# create the A string
y = guitar(110, fs).play(d)
Audio.Audio(data=y, rate=fs, embed=True)


s = abs(np.fft.fftpack.fft(y));
s = s[0:len(s)/2]
plot(np.linspace(0,1,len(s))*(fs/2), s)

from scipy import signal

class guitar:
    def __init__(self, pitch=110, fs=24000):
        # init the class with desired pitch and underlying sampling frequency
        self.M = np.round(fs / pitch) # fundamental period in samples
        self.R = 0.9999               # decay factor
        self.RM = self.R ** self.M 
        self.ybuf = np.zeros(self.M)  # output buffer (circular)
        self.iy = 0                   # index into out buf
        self.xbuf = 0                 # input buffer (just one sample)
        # 6th-order Butterworth, keep 5 harmonics:
        self.bfb, self.bfa = signal.butter(6, min(0.5, 5.0 * pitch / fs))
        self.bfb *= 1000              # set a little gain 
        # initial conditions for the filter. We need this because we need to
        # filter on a sample-by-sample basis later on
        self.bfs = signal.lfiltic(self.bfb, self.bfa, [0])
        
    def play(self, x):
        y = np.zeros(len(x))
        for n in range(len(x)):
            # comb filter
            t = x[n] - self.R * self.xbuf + self.RM * self.ybuf[self.iy]
            self.ybuf[self.iy] = t
            self.iy = (self.iy + 1) % self.M
            self.xbuf = x[n]
            # lowpass filter, keep filter status for next sample
            y[n], self.bfs = signal.lfilter(self.bfb, self.bfa, [t], zi=self.bfs)
        return y

y = guitar(110, fs).play(d)
Audio.Audio(data=y, rate=fs, embed=True)

s = abs(np.fft.fftpack.fft(y[10000:30000]));
s = s[0:len(s)/2]
plot(np.linspace(0,1,len(s))*(fs/2), s)

def amplify(x):
    TH = 0.9           # threshold
    y = copy(x)
    y[y >  TH] =  TH
    y[y < -TH] = -TH
    return y

x = np.linspace(-2, 2, 100)
plot(x, amplify(x))
xlabel("input")
ylabel("output")

display(Image(filename='data/specgram.png', width=800))

class feedback:
    SPEED_OF_SOUND = 343.0 # m/s
    def __init__(self, max_distance_m = 5, fs=24000):  
        # init class with maximum distance
        self.L = ceil(max_distance_m / self.SPEED_OF_SOUND * fs);
        self.xbuf = zeros(self.L)     # circular buffer
        self.ix = 0
        
    def get(self, x, distance):
        d = ceil(distance / self.SPEED_OF_SOUND * fs)    # delay in samples
        self.xbuf[self.ix] = x
        x = self.xbuf[(self.L + self.ix - d) % self.L]
        self.ix = (self.ix + 1) % self.L
        return x / float(distance)

g = guitar(110)    # the A string
f = feedback()     # the feedback channel

# the "coupling loss" between air and string is high. Let's say that
#  it is about 80dBs
COUPLING_LOSS = 0.0001

# John starts 3m away and then places the guitar basically against the amp
#  after 1.5 seconds
START_DISTANCE = 3 
END_DISTANCE = 0.05

N = fs * 5         # play for 5 seconds
y = zeros(N)     
x = [1]            # the initial plucking
# now we create each sample in a loop by processing the guitar sound
# thru the amp and then feeding back the attenuated and delayed sound
# to the guitar
for n in range(N):
    y[n] = amplify(g.play(x))
    x = [COUPLING_LOSS * f.get(y[n], START_DISTANCE if n < (1.5 * fs) else END_DISTANCE)]
       
        
Audio.Audio(data=y, rate=fs, embed=True)