 Program Talk - Source Code Browser
```# Compute the audio novelty features of a song

import sys

import numpy as np
import scipy

from ..utils import RMS_energy

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
: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))
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
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 * wlen_ms / 2000.0, x) for x in out_peaks]
for i, p in enumerate(out_peaks):
if i == nchangepoints:
break

return [x for x in out_peaks[:nchangepoints]]

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

def smooth_hanning(x, size=11):
"""smooth a 1D array using a hanning window with requested size."""

if x.ndim != 1:
raise ValueError, "smooth_hanning only accepts 1-D arrays."
if x.size < size:
raise ValueError, "Input vector needs to be bigger than window size."
if size < 3:
return x

s = np.r_[x[size - 1:0:-1], x, x[-1:-size:-1]]
w = np.hanning(size)
y = np.convolve(w / w.sum(), s, mode='valid')
return y

def naive_peaks(vec, k=33):
"""A naive method for finding peaks of a signal.
1. Smooth vector
2. Find peaks (local maxima)
3. Find local max from original signal, pre-smoothing
4. Return (sorted, descending) peaks
"""

a = smooth_hanning(vec, k)

k2 = (k - 1) / 2

peaks = np.r_[True, a[1:] > a[:-1]] & np.r_[a[:-1] > a[1:], True]

p = np.array(np.where(peaks))
maxidx = np.zeros(np.shape(p))
maxvals = np.zeros(np.shape(p))
for i, pk in enumerate(p):
maxidx[i] = np.argmax(vec[pk - k2:pk + k2]) + pk - k2
maxvals[i] = np.max(vec[pk - k2:pk + k2])
out = np.array([maxidx, maxvals]).T
return out[(-out[:, 1]).argsort()]

if __name__ == '__main__':
novelty(sys.argv, k=int(sys.argv))
```