Function for creating Huggins pitch stimuli.

huggins.py

# -*- coding: utf-8 -*-
"""
Created on Wed Nov 20 15:22:22 2013

@author: smathias
"""

import brian
import brian.hears as b
import numpy as np
from scipy.signal import lfilter, butter
import matplotlib.pyplot as plt

def huggins(f, dur, c=0.03):
    """
    Create a Huggins-pitch stimulus with a pitch at f in Hz and duration dur in
    seconds. If f is None, normal white noise is made instead.
    
    """
    if f:
        snd1 = b.whitenoise(dur*brian.second)
        fs = float(snd1.samplerate)
        snd1 = np.ravel(np.array(snd1))
        r = 1 - ((c * np.pi * f) / fs)
        br = np.array([(r**2), -2 * r*np.cos(2 * np.pi * f / fs), 1])
        a = br[::-1]
        snd2 = lfilter(br.T, a.T, np.array(snd1))
        snd = b.Sound((snd1,snd2))
    else:
        snd = b.whitenoise(dur*brian.second)
        snd = b.Sound((snd,snd))
    return snd

def butterworth(snd, lowcut, highcut, order=4, pad=0.5):
    nyq = 0.5 * int(snd.samplerate)
    low = lowcut / nyq
    high = highcut / nyq
    br, a = butter(order, [low, high], btype='band')
    x1 = np.ravel(np.array(snd.left))
    y1 = lfilter(br, a, x1)
    x2 = np.ravel(np.array(snd.right))
    y2 = lfilter(br, a, x2)
    return b.Sound((y1,y2))

def main():
    low, high = (200, 2000)
    F = [587.33, None, 659.255, None, 523.251,
         None, 261.626, None, 391.995, None]
    snds = [huggins(f, 0.3) for f in F]
    snd = b.sequence(snds)
    snd = butterworth(snd, low, high) # filter for aesthetic  value
    snd.level = 60 * b.dB
    snd.play(sleep=True)
    snd.save('CloseEncountHp.wav')
    plt.subplot(211)
    plt.title('Left')
    snd.left.spectrogram(high=8*brian.kHz)
    plt.xlabel('')
    plt.subplot(212)
    plt.title('Right')
    snd.right.spectrogram(high=8*brian.kHz)
    plt.show()

if __name__ == '__main__':
    main()