.. currentmodule:: brian

.. index::
   pair: example usage; Sound
   pair: example usage; run
   pair: example usage; Network
   pair: example usage; sequence
   pair: example usage; TanCarney
   pair: example usage; ZhangSynapse
   pair: example usage; exp
   pair: example usage; set_default_samplerate
   pair: example usage; MiddleEar
   pair: example usage; tone
   pair: example usage; StateMonitor

.. _example-hears-tan_carney_2003_tan_carney_Fig11:

Example: tan_carney_Fig11 (hears/tan_carney_2003)
=================================================

Response area and phase response of a model fiber with CF=2200Hz in the 
Tan&Carney model. Reproduces Fig. 11 from:

Tan, Q., and L. H. Carney.
    "A Phenomenological Model for the Responses of Auditory-nerve Fibers.
    II. Nonlinear Tuning with a Frequency Glide".
    The Journal of the Acoustical Society of America 114 (2003): 2007.

::

    
    
    import itertools
    
    import matplotlib.pyplot as plt
    import numpy as np
    
    from brian import *
    # set_global_preferences(useweave=True)
    from brian.hears import *
    
    duration = 50*ms
    samplerate = 50*kHz
    set_default_samplerate(samplerate)
    CF = 2200
    freqs = np.arange(250.0, 3501., 50.)
    levels = [10, 30, 50, 70, 90]
    cf_level = list(itertools.product(freqs, levels))
    tones = Sound([Sound.sequence([tone(freq * Hz, duration).atlevel(level*dB).ramp(when='both',
                                                                                    duration=2.5*ms,
                                                                                    inplace=False)])
                   for freq, level in cf_level])
    
    ihc = TanCarney(MiddleEar(tones), [CF] * len(cf_level), update_interval=2)
    syn = ZhangSynapse(ihc, CF)
    s_mon = StateMonitor(syn, 's', record=True, clock=syn.clock)
    net = Network(syn, s_mon)
    net.run(duration)
    
    reshaped = s_mon.values.reshape((len(freqs), len(levels), -1))
    
    # calculate the phase with respect to the stimulus
    pi = np.pi
    min_freq, max_freq = 1100, 2900
    freq_subset = freqs[(freqs>=min_freq) & (freqs<=max_freq)]
    reshaped_subset = reshaped[(freqs>=min_freq) & (freqs<=max_freq), :, :]
    phases = np.zeros((reshaped_subset.shape[0], len(levels)))
    for f_idx, freq in enumerate(freq_subset):
        period = 1.0 / freq
        for l_idx in xrange(len(levels)):
            phase_angles = np.arange(reshaped_subset.shape[2])/samplerate % period / period * 2*pi
            temp_phases = (np.exp(1j * phase_angles) *
                           reshaped_subset[f_idx, l_idx, :])
            phases[f_idx, l_idx] = np.angle(np.sum(temp_phases))
    
    plt.subplot(2, 1, 1)
    rate = reshaped.mean(axis=2)
    plt.plot(freqs, rate)
    plt.ylabel('Spikes/sec')
    plt.legend(['%.0f dB' % level for level in levels], 0)
    plt.xlim(0, 4000)
    plt.ylim(0, 250)
    
    plt.subplot(2, 1, 2)
    relative_phases = (phases.T - phases[:, -1]).T
    relative_phases[relative_phases > pi] = relative_phases[relative_phases > pi] - 2*pi
    relative_phases[relative_phases < -pi] = relative_phases[relative_phases < -pi] + 2*pi 
    plt.plot(freq_subset, relative_phases / pi)
    plt.ylabel("Phase Re:90dB (pi radians)")
    plt.xlabel('Frequency (Hz)')
    plt.legend(['%.0f dB' % level for level in levels], 0)
    plt.xlim(0, 4000)
    plt.ylim(-0.5, 0.75)
    plt.show()