%pylab
%reset
from brian import *
set_global_preferences(useweave=True)
#set_global_preferences(usegpu=True)
from brian.hears import *

duration = .1*second
samplerate = 44.1*kHz
sound1 = tone(1*kHz, duration,samplerate = samplerate)
sound2 = whitenoise(duration,samplerate = samplerate)
taxis = arange(0*ms,duration,1./samplerate)
sound = sound1+sound2
sound = sound.ramp()
#sound = loadsound('test.wav')

#sound = sound.extended(20*ms)
#sound = sound.repeat(5)
#sound.level = 50*dB

figure()
subplot(311)
plot(taxis,sound1)
subplot(312)
plot(taxis,sound2)
subplot(313)
plot(taxis,sound)
xlabel('time axis (s)')
show()


#sound = powerlawnoise(200*ms, 1, samplerate=44100*Hz)
#sound = whitenoise(duration,samplerate = samplerate)
sound = harmoniccomplex(1000*Hz,100*ms)
sound = sound.ramp(when='both', duration=10.0 * msecond, envelope=lambda t :sin(pi*t/2)**2)
#sound = sound.extended(20*ms)
#sound = sound.repeat(5)
sound.level = 50*dB
figure()
subplot(211)
sound.spectrogram()
subplot(212)
plot(sound)
show()

cf = erbspace(20*Hz, 20*kHz, 3000)
basilar_membrane = Gammatone(sound, cf)
half_wave_rect = lambda x: 3*clip(x, 0, Inf)**(1.0/3.0)
hair_cell1 = FunctionFilterbank(basilar_membrane,half_wave_rect)
hair_cell = LowPass(hair_cell1, 800*Hz)

out_hair_cell = hair_cell.process()


imshow(out_hair_cell.T, origin='lower left', aspect='auto', vmin=0,extent=(0, sound.duration/ms,cf[0], cf[-1]))
ylabel('Frequency (Hz)')
xlabel('Time (ms)')
show()

# Leaky integrate-and-fire model with noise and refractoriness
tau = 1*ms
eqs = '''
dv/dt = (I-v)/(tau)+0.2*xi*(2/(tau))**.5 : 1
I : 1
'''
anf = FilterbankGroup(hair_cell1, 'I', eqs, reset=0, threshold=1, refractory=2*ms)

M = SpikeMonitor(anf)
run(sound.duration)
figure()
raster_plot(M)
show()

sound.level = 50*dB
linear_path = Gammatone(sound, cf)

bandpass_nonlinear1 = Gammatone(sound, cf+10)
func_compression = lambda x: sign(x)*abs(x)**0.3
non_linear_path = FunctionFilterbank(bandpass_nonlinear1, func_compression)

total_path = linear_path+non_linear_path
output = total_path.process()

figure()
imshow(output.T, origin='lower left', aspect='auto', vmin=0,extent=(0, sound.duration/ms,cf[0], cf[-1]))
ylabel('Frequency (Hz)')
xlabel('Time (ms)')
show()

sound.level = 100*dB
cf = erbspace(20*Hz, 2*kHz, 300)
param_drnl = {'lp_nl_cutoff_m':1.1}
drnl_filter=DRNL(sound, cf, type='human', param=param_drnl)
output = drnl_filter.process()

imshow(output.T, origin='lower left', aspect='auto', vmin=0,extent=(0, sound.duration/ms,cf[0], cf[-1]))
ylabel('Frequency (Hz)')
xlabel('Time (ms)')
show()

nsounds = 3
sound = Sound(randn(4000, nsounds), samplerate=16*kHz)

nchannels = 100
cf = erbspace(20*Hz, 10*kHz, nchannels)
cf = tile(cf,nsounds)

sound = Repeat(sound,nchannels)

gf = Gammatone(sound, cf)

output = gf.process()

imshow(output.T, origin='lower left', aspect='auto', vmin=0,extent=(0, sound.duration/ms,cf[0], cf[-1]))
ylabel('Frequency (Hz)')
xlabel('Time (ms)')
show()
example: compute online rms value

Filterbank.process() method allows us to pass an optional function f(input, running) 
process() will first call running = f(output, 0) for the first buffered segment input
It will then call running = f(output, running) for each subsequent segment


def sum_of_squares(input, running):
    return running+sum(input**2, axis=0)
nsamples = 100
rms = sqrt(gf.process(sum_of_squares)/nsamples)
print rms