Return points of high "novelty" in a song
    (e.g., significant musical transitions)

    :param song: Song to analyze
    :type song: :py:class:`radiotool.composer.Song`
    :param k: Width of comparison kernel (larger kernel finds coarser differences in music)
    :type k: int
    :param wlen_ms: Analysis window length in milliseconds
    :type wlen_ms: int
    :param start: Where to start analysis within the song (in seconds)
    :type start: float
    :param duration: How long of a chunk of the song to analyze (None analyzes the entire song after start)
    :type duration: float
    :param nchangepoints: How many novel change points to return
    :type nchangepoints: int
    :param feature: Music feature to use for novelty analysis
    :type feature: "rms" or "mfcc" (will support "chroma" eventually)
    :returns: List of change points (in seconds)
    :rtype: list of floats

def novelty(song, k=64, wlen_ms=100, start=0, duration=None, nchangepoints=5, feature="rms"):
    """Return points of high "novelty" in a song
    (e.g., significant musical transitions)

    :param song: Song to analyze
    :type song: :py:class:`radiotool.composer.Song`
    :param k: Width of comparison kernel (larger kernel finds coarser differences in music)
    :type k: int
    :param wlen_ms: Analysis window length in milliseconds
    :type wlen_ms: int
    :param start: Where to start analysis within the song (in seconds)
    :type start: float
    :param duration: How long of a chunk of the song to analyze (None analyzes the entire song after start)
    :type duration: float
    :param nchangepoints: How many novel change points to return
    :type nchangepoints: int
    :param feature: Music feature to use for novelty analysis
    :type feature: "rms" or "mfcc" (will support "chroma" eventually)
    :returns: List of change points (in seconds)
    :rtype: list of floats
    """

    if feature != "rms" and feature != "mfcc":
        raise ValueError, "novelty currently only supports 'rms' and 'mfcc' features"

    if feature == "rms":
        frames = song.all_as_mono()

        wlen_samples = int(wlen_ms * song.samplerate / 1000)

        if duration is None:
            frames = frames[start * song.samplerate:]
        else:
            frames = frames[start * song.samplerate:(start + duration) *
                 song.samplerate]
                 
        # Compute energies
        hamming = np.hamming(wlen_samples)
        nwindows = int(2 * song.duration / wlen_samples - 1)
        energies = np.empty(nwindows)
        for i in range(nwindows):
            energies[i] = RMS_energy(
                hamming *
                frames[i * wlen_samples / 2:
                       i * wlen_samples / 2 + wlen_samples]
            )

        energies_list = [[x] for x in energies]
    elif feature == "mfcc":
        analysis = song.analysis
        energies_list = np.array(analysis["timbres"])
    
    # Compute similarities
    S_matrix = 1 - scipy.spatial.distance.squareform(
        scipy.spatial.distance.pdist(energies_list, 'euclidean'))
                    
    # smooth the C matrix with a gaussian taper
    C_matrix = np.kron(np.eye(2), np.ones((k,k))) -\
               np.kron([[0, 1], [1, 0]], np.ones((k,k)))
    g = scipy.signal.gaussian(2*k, k)
    C_matrix = np.multiply(C_matrix, np.multiply.outer(g.T, g))
    
    # Created checkerboard kernel
    
    N_vec = np.zeros(np.shape(S_matrix)[0])
    for i in xrange(k, len(N_vec) - k):
        S_part = S_matrix[i - k:i + k, i - k:i + k]
        N_vec[i] = np.sum(np.multiply(S_part, C_matrix))
        
    # Computed checkerboard response

    peaks = naive_peaks(N_vec, k=k / 2 + 1)
    out_peaks = []
    
    if feature == "rms":
        # ensure that the points we return are more exciting
        # after the change point than before the change point
        for p in peaks:
            frame = p[0]
            if frame > k:
                left_frames = frames[int((frame - k) * wlen_samples / 2):
                                     int(frame * wlen_samples / 2)]
                right_frames = frames[int(frame * wlen_samples / 2):
                                      int((frame + k) * wlen_samples / 2)]
                if RMS_energy(left_frames) <\
                   RMS_energy(right_frames):
                   out_peaks.append(p)

        out_peaks = [(x[0] * wlen_ms / 2000.0, x[1]) for x in out_peaks]
        for i, p in enumerate(out_peaks):
            if i == nchangepoints:
                break
        
        return [x[0] for x in out_peaks[:nchangepoints]]

    elif feature == "mfcc":
        beats = analysis["beats"]
        return [beats[int(b[0])] for b in peaks[:nchangepoints]]