<file_sep>from __future__ import print_function, division

import numpy as np

class Tuning():
    """
    Equal temperament tuning - allows to convert between frequency and pitch.

    - unit pitch space
      - continous, unbounded
      - 1.0 ~ one octave
    - step pitch space
      - continous, unbounded
      - N steps ~ one octave
      - unit pitch space * N
    - unit pitch class space
      - continous, bounded [0, 1.0)
      - unit pitch space % 1.0
    - step pitch class space
      - continous, bounded [0, N)
      - unit step pitch space % N
    - integer step pitch space
      - discrete, unbounded
      - floor(step pitch space)
    - integer step pitch class space
      - discrete, bounded {0, 1, .. N - 1}
      - floor(step pitch class space)
    """
    def __init__(self, base_freq=440, steps_per_octave=12, octave_ratio=2):
        self.base_freq = base_freq
        self.steps_per_octave = steps_per_octave
        self.octave_ratio = octave_ratio

    def pitch_to_freq(self, pitch):
        factor = self.pitch_to_relative_freq(pitch)
        return factor * self.base_freq

    def freq_to_pitch(self, freq):
        rel_freq = freq / self.base_freq
        if self.octave_ratio == 2:
            p = np.log2(rel_freq)
        else:
            p = np.log(rel_freq) / np.log(2)
        return p * self.steps_per_octave

    def pitch_to_relative_freq(self, pitch):
        return pow(self.octave_ratio, pitch / self.steps_per_octave)

class PitchQuantizer():
    def __init__(self, tuning, bin_division=1):
        self.tuning = tuning
        self.bin_division = bin_division

    def quantize(self, freqs):
        """
        Quantizes frequencies to nearest pitch bins (with optional division of
        bins).
        """
        return np.round(self.tuning.freq_to_pitch(freqs) * self.bin_division) / self.bin_division
<file_sep>from __future__ import print_function, division

import math
import os

import numpy as np
import scipy

from .spectrogram import db_scale, positive_freq_magnitudes, \
    select_positive_freq_fft, fftfreqs, normalized_window, scale_magnitudes
from .signal import SignalFrames
from .tuning import PitchQuantizer, Tuning
from .plots import save_raw_spectrogram_bitmap


class LinearTransform():
    def __init__(self, positive_only=True):
        # range of normalized frequencies
        self.bin_range = (0, 0.5) if positive_only else (0, 1)

    def transform_freqs(self, X_inst_freqs, sample_rate):
        output_bin_count = X_inst_freqs.shape[1]
        X_y = X_inst_freqs
        return X_y, output_bin_count, self.bin_range


class PitchTransform():
    """
    Perform the proper quantization to pitch bins according to possible
    subdivision before the actual histogram computation. Still we need to
    move the quantized pitch value a bit from the lower bin edge to ensure
    proper floating point comparison. Note that the quantizer rounds values
    from both sides towards the quantized value, while histogram2d floors the
    values to the lower bin edge. The epsilon is there to prevent log of 0
    in the pitch to frequency transformation.

    bin_range: range of pitch bins (default: A0 27.5 Hz to E10 21096.16 Hz)
    """
    def __init__(self, bin_range=(-48, 67), bin_division=1, tuning=Tuning()):
        self.tuning = tuning
        self.bin_range = bin_range
        self.bin_division = bin_division

    def transform_freqs(self, X_inst_freqs, sample_rate):
        quantization_border = 1 / (2 * self.bin_division)
        pitch_quantizer = PitchQuantizer(self.tuning, bin_division=self.bin_division)
        eps = np.finfo(np.float32).eps
        # TODO: is it possible to quantize using relative freqs to avoid
        # dependency on the fs parameter?
        X_y = pitch_quantizer.quantize(np.maximum(sample_rate * X_inst_freqs, eps) + quantization_border)
        output_bin_count = (self.bin_range[1] - self.bin_range[0]) * self.bin_division
        return X_y, output_bin_count, self.bin_range


class Spectrogram():
    """
    Represents spectrogram information of a time-domain signal which can be used
    to compute various types of reassigned spectrograms, pitchgrams, etc.
    """
    def __init__(self, signal_frames, window=scipy.hanning, positive_only=True):
        """
        :param signal_frames: signal represented as SignalFrames instance
        :param window: STFT window function - produces 1D window which will
        be normalized
        """
        self.signal_frames = signal_frames

        x_frames = signal_frames.frames
        w = normalized_window(window(signal_frames.frame_size))

        # complex spectra of windowed blocks of signal - STFT
        self.X_complex = np.fft.fft(x_frames * w)
        # linear magnitude spectrogram
        self.X_mag = abs(self.X_complex) / self.X_complex.shape[1]

        # spectra of signal shifted in time

        # This fakes looking at the previous frame shifted by one sample.
        # In order to work only with one frame of size N and not N + 1, we fill the
        # missing value with zero. This should not introduce a large error, since the
        # borders of the amplitude frame will go to zero anyway due to applying a
        # window function in the STFT tranform.
        X_prev_time = np.fft.fft(shift_right(x_frames) * w)

        # spectra shifted in frequency
        X_prev_freq = shift_right(self.X_complex)

        # cross-spectra - ie. spectra of cross-correlation between the
        # respective time-domain signals
        X_cross_time = cross_spectrum(self.X_complex, X_prev_time)
        X_cross_freq = cross_spectrum(self.X_complex, X_prev_freq)

        # instantaneous frequency estimates
        # normalized frequencies in range [0.0, 1.0] - from DC to sample rate
        self.X_inst_freqs = estimate_instant_freqs(X_cross_time)
        # instantaneous group delay estimates
        # relative coordinates within the frame with range [-0.5, 0.5] where
        # 0.0 is the frame center
        self.X_group_delays = estimate_group_delays(X_cross_freq)

        if positive_only:
            self.X_mag = positive_freq_magnitudes(self.X_mag)
            self.X_complex, self.X_inst_freqs, self.X_group_delays = [
                select_positive_freq_fft(values) for values in
                [self.X_complex, self.X_inst_freqs, self.X_group_delays]
            ]

    def reassigned(
        self,
        output_frame_size=None, transform=LinearTransform(),
        reassign_time=True, reassign_frequency=True, magnitudes='power_db'):
        """
        Reassigned spectrogram requantized both in frequency and time.

        Note it is quantized into non-overlapping output time frames which may be
        of a different size than input time frames.

        transform - transforms the frequencies
        """

        if output_frame_size is None:
            output_frame_size = self.signal_frames.hop_size

        frame_size = self.signal_frames.frame_size
        fs = self.signal_frames.sample_rate

        frame_duration = frame_size / fs
        frame_center_time = frame_duration / 2
        # group delays are in range [-0.5, 0.5] - relative coordinates within the
        # frame where 0.0 is the frame center
        input_bin_count = self.X_inst_freqs.shape[1]

        eps = np.finfo(np.float32).eps
        X_time = np.tile(self.signal_frames.start_times + frame_center_time +
            eps, (input_bin_count, 1)).T
        if reassign_time:
            X_time += self.X_group_delays * frame_duration

        if reassign_frequency:
            X_y = self.X_inst_freqs
        else:
            X_y = np.tile(fftfreqs(frame_size, fs) / fs, (self.X_inst_freqs.shape[0], 1))

        X_y, output_bin_count, bin_range = transform.transform_freqs(X_y,
            self.signal_frames.sample_rate)
        frame_duration = frame_size / fs
        end_input_time = self.signal_frames.duration
        output_frame_count = int(math.ceil((end_input_time * fs) / output_frame_size))
        print('output_frame_count', output_frame_count)
        time_range = (0, output_frame_count * output_frame_size / fs)

        output_shape = (output_frame_count, output_bin_count)
        X_spectrogram, x_edges, y_edges = np.histogram2d(
            X_time.flatten(), X_y.flatten(),
            weights=self.X_mag.flatten(),
            range=(time_range, bin_range),
            bins=output_shape)

        X_spectrogram = scale_magnitudes(X_spectrogram, magnitudes)

        return X_spectrogram

def cross_spectrum(spectrumA, spectrumB):
    """
    Returns a cross-spectrum, ie. spectrum of cross-correlation of two signals.
    This result does not depend on the order of the arguments.
    Since we already have the spectra of signals A and B and and want the
    spectrum of their cross-correlation, we can replace convolution in time
    domain with multiplication in frequency domain.
    """
    return spectrumA * spectrumB.conj()

def shift_right(values):
    """
    Shifts the array to the right by one place, filling the empty values with
    zeros.
    TODO: use np.roll()
    """
    # TODO: this fails for 1D input array!
    return np.hstack([np.zeros((values.shape[0], 1)), values[..., :-1]])

def arg(values):
    """
    Argument (angle) of complex numbers wrapped and scaled to [0.0, 1.0].

    input: an array of complex numbers
    output: an array of real numbers of the same shape

    np.angle() returns values in range [-np.pi, np.pi].
    """
    return np.mod(np.angle(values) / (2 * np.pi), 1.0)

def estimate_instant_freqs(crossTimeSpectrum):
    """
    Channelized instantaneous frequency - the vector of simultaneous
    instantaneous frequencies computed over a single frame of the digital
    short-time Fourier transform.

    Instantaneous frequency - derivative of phase by time.

    cif = angle(crossSpectrumTime) * sampleRate / (2 * pi)

    In this case the return value is normalized (not multiplied by sampleRate)
    to the [0.0; 1.0] interval, instead of absolute [0.0; sampleRate].
    """
    return arg(crossTimeSpectrum)

def estimate_group_delays(crossFreqSpectrum):
    "range: [-0.5, 0.5]"
    return 0.5 - arg(crossFreqSpectrum)


def process_spectrogram(filename, frame_size, hop_size, output_frame_size):
    """
    Computes three types of spectrograms (normal, frequency reassigned,
    time-frequency reassigned) from an audio file and stores and image from each
    spectrogram into PNG file.
    """
    signal_frames = SignalFrames(filename, frame_size, hop_size, mono_mix=True)

    spectrogram = Spectrogram(signal_frames)

    image_filename = os.path.basename(filename).replace('.wav', '')

    # STFT on overlapping input frames
    X_stft = db_scale(spectrogram.X_mag ** 2)
    save_raw_spectrogram_bitmap(image_filename + '_stft_frames.png', X_stft)

    linear_transform = LinearTransform(positive_only=True)

    # STFT requantized to the output frames (no reassignment)
    X_stft_requantized = spectrogram.reassigned(output_frame_size,
        linear_transform,
        reassign_time=False, reassign_frequency=False)
    save_raw_spectrogram_bitmap(image_filename + '_stft_requantized.png', X_stft_requantized)

    # STFT reassigned in time and requantized to output frames
    X_reassigned_t = spectrogram.reassigned(output_frame_size,
        linear_transform,
        reassign_time=True, reassign_frequency=False)
    save_raw_spectrogram_bitmap(image_filename + '_reassigned_t.png', X_reassigned_t)

    # STFT reassigned in frequency and requantized to output frames
    X_reassigned_f = spectrogram.reassigned(output_frame_size,
        linear_transform,
        reassign_time=False, reassign_frequency=True)
    save_raw_spectrogram_bitmap(image_filename + '_reassigned_f.png', X_reassigned_f)

    # STFT reassigned both in time and frequency and requantized to output frames
    X_reassigned_tf = spectrogram.reassigned(output_frame_size,
        linear_transform,
        reassign_time=True, reassign_frequency=True)
    save_raw_spectrogram_bitmap(image_filename + '_reassigned_tf.png', X_reassigned_tf)

    pitch_transform = PitchTransform(bin_range=(-48, 67), bin_division=1)

    # TF-reassigned pitchgram
    X_pitchgram_tf = spectrogram.reassigned(output_frame_size,
        pitch_transform,
        reassign_time=True, reassign_frequency=True)
    save_raw_spectrogram_bitmap(image_filename + '_pitchgram_tf.png', X_pitchgram_tf)

    # T-reassigned pitchgram
    X_pitchgram_t = spectrogram.reassigned(output_frame_size,
        pitch_transform,
        reassign_time=True, reassign_frequency=False)
    save_raw_spectrogram_bitmap(image_filename + '_pitchgram_t.png', X_pitchgram_t)

    # F-reassigned pitchgram
    X_pitchgram_t = spectrogram.reassigned(output_frame_size,
        pitch_transform,
        reassign_time=False, reassign_frequency=True)
    save_raw_spectrogram_bitmap(image_filename + '_pitchgram_f.png', X_pitchgram_t)

    # non-reassigned pitchgram
    X_pitchgram = spectrogram.reassigned(output_frame_size,
        pitch_transform,
        reassign_time=False, reassign_frequency=False)
    save_raw_spectrogram_bitmap(image_filename + '_pitchgram_no.png', X_pitchgram)

def reassigned_spectrogram(signal_frames, output_frame_size=None, magnitudes='power_db',
    reassign_time=True, reassign_frequency=True):
    """
    From frames of audio signal it computes the frequency reassigned spectrogram
    requantized back to the original linear bins.

    Only the real half of spectrum is given.
    """
    return Spectrogram(signal_frames).reassigned(
        output_frame_size, LinearTransform(),
        reassign_time, reassign_frequency, magnitudes=magnitudes)

# [-48,67) -> [~27.5, 21096.2) Hz
def pitchgram(signal_frames, output_frame_size=None, bin_range=(-48, 67), bin_division=1, magnitudes='power_db'):
    """
    From frames of audio signal it computes the frequency reassigned spectrogram
    requantized to pitch bins (pitchgram).
    """
    return Spectrogram(signal_frames).reassigned(
        output_frame_size, PitchTransform(bin_range, bin_division), magnitudes=magnitudes)

if __name__ == '__main__':
    import sys
    process_spectrogram(filename=sys.argv[1], frame_size=4096, hop_size=1024, output_frame_size=1024)
<file_sep>from __future__ import print_function, division

import numpy as np
import math
import scipy

from .signal import mean_power
from .signal import SignalFrames

def spectrogram(filename, frame_size=2048, hop_size=512, magnitudes='power_db'):
    """
    Computes an STFT magnitude power spectrogram from an audio file.
    Returns: spectrogram, audio_samples, frame_times
    """
    signal_frames = SignalFrames(filename, frame_size, hop_size, mono_mix=True)
    x = signal_frames.frames
    times = signal_frames.start_times
    w = create_window(frame_size)
    X = stft_spectrogram(x, w, magnitudes)
    return X, x, times

def stft_spectrogram(x, w, magnitudes):
    """
    Computes an STFT magnitude power spectrogram from an array of samples
    already cut to frames.
    Input:
    - x - time-domain samples - array of shape (frames, frame_size)
    - w - window - array of shape (frame_size)
    - magnitudes - indicates whether to scale the
    Output: spectrogram
    """
    X = magnitude_spectrum(x * w) ** 2
    if magnitudes:
        X = db_scale(X)
    return X

def magnitude_spectrum(x):
    '''
    Magnitude spectrum scaled so that each bin corresponds to the original sine
    amplitude. Only the real part of the spectrum is returned.
    x - 1D sampled signal (possibly already windowed)

    For signal in range [-1., 1.] the output range is [0., 1.].

    The energy is not preserved, it's scaled down
    (energy_out = energy_in / (N//2)).
    '''
    X = np.fft.fft(x)
    Xr = positive_freq_magnitudes(X)
    N = Xr.shape[-1]
    return abs(Xr) / N

def select_positive_freq_fft(X):
    """
    Select the positive frequency part of the spectrum in a spectrogram.
    """
    N = X.shape[1]
    return X[:, :N//2]

# TODO: we should probably multiply the whole result by 2, to conserve energy
def positive_freq_magnitudes(X):
    """
    Select magnitudes from positive-frequency half of the spectrum in a
    spectrogram. The DC term shared for positive and negative frequencies is
    halved.

    Note this is not a complete information to reconstruct the full spectrum,
    since we throw away the bin at the negative Nyquist frequency (index N/2+1).
    """
    X_pos = select_positive_freq_fft(X).copy()
    X_pos[:, 0] *= 0.5
    return X_pos

def create_window(size):
    """
    A normalized Hanning window of given size. Useful for analyzing sinusoidal
    signals.
    """
    return normalized_window(scipy.hanning(size))

def normalized_window(w):
    """
    Normalizes an FFT window so that it has energy equal to its length, and mean
    power equal to 1.0.
    """
    return w / mean_power(w)

def db_scale(magnitude_spectrum, normalized=False):
    """
    Transform linear magnitude to dbFS (full-scale) [-120, 0] (for input range
    [0.0, 1.0]) which can be optionally normalized to [0.0, 1.0].
    """
    scaled = 20 * np.log10(np.maximum(1e-6, magnitude_spectrum))

    # map from raw dB [-120.0, 0] to [0.0, 1.0]
    if normalized:
        scaled = (scaled / 120) + 1
    return scaled

def scale_magnitudes(X_mag, transform):
    if transform == 'linear':
        return X_mag
    elif transform == 'power':
        return X_mag ** 2
    elif transform == 'power_db':
        return db_scale(X_mag ** 2)
    elif transform == 'power_db_normalized':
        return db_scale(X_mag ** 2, normalized=True)
    else:
        raise ValueError('Unknown magnitude scaling transform ' + transform)

# -- extras --

def energy_weighted_spectrum(x):
    N = x.shape[-1]
    X = np.fft.fft(x)
    # np.allclose(energy(abs(X) / math.sqrt(N)), energy(x))
    # np.allclose(energy(abs(X[:N//2]) / math.sqrt(N//2)), energy(x))
    return abs(X) / math.sqrt(N)

def fftfreqs(frame_size, fs):
    """
    Positive FFT frequencies from DC (incl.) until Nyquist (excl.).
    The size of half of the FTT size.
    """
    return np.fft.fftfreq(frame_size, 1/fs)[:frame_size // 2]

def inverse_spectrum(spectrum, window):
    '''
    inverse_spectrum(np.fft.fft(x * window), window) == x
    '''
    return np.real(np.fft.ifft(spectrum)) / window
<file_sep>from __future__ import print_function, division

import numpy as np
import os

from tfr import SignalFrames, Tuning, pitchgram, reassigned_spectrogram
from tfr.reassignment import shift_right, arg

DATA_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')

def test_shift_right():
    assert np.allclose(shift_right(np.array([[1, 2, 3]])), np.array([0, 1, 2]))

def test_pitchgram_on_single_tone_should_have_peak_at_that_tone():
    pitch = 12 + 7 # G5
    f = Tuning().pitch_to_freq(pitch)
    fs = 44100
    x = sine(sample_time(0, 1, fs=fs), freq=f)
    frame_size = 4096
    hop_size = 2048
    output_frame_size = hop_size
    signal_frames = SignalFrames(x, frame_size, hop_size, sample_rate=fs, mono_mix=True)
    bin_range = [-48, 67]
    x_pitchgram = pitchgram(signal_frames,
        output_frame_size, magnitudes='power_db', bin_range=bin_range, bin_division=1)

    max_bin_expected = pitch - bin_range[0]
    max_bin_actual = x_pitchgram.mean(axis=0).argmax()

    assert x_pitchgram.shape == (22, 115), x_pitchgram.shape
    assert max_bin_actual == max_bin_expected

def test_arg():
    values = np.array([-5.-1.j, -1.-5.j,  2.-2.j,  3.+4.j,  2.+0.j,  2.-5.j, -3.-3.j,
       -3.+1.j, -2.-5.j,  0.+2.j])
    args = arg(values)
    expected_args=np.array([0.53141648, 0.71858352, 0.875 , 0.14758362, 0.,
        0.81055947, 0.625 , 0.44879181, 0.68944053, 0.25])
    assert np.allclose(args, expected_args)

def test_reassigned_spectrogram_values_should_be_in_proper_range():
    frame_size = 4096
    hop_size = frame_size
    output_frame_size = 1024
    audio_file = os.path.join(DATA_DIR, 'she_brings_to_me.wav')
    signal_frames = SignalFrames(audio_file, frame_size, hop_size, mono_mix=True)
    X_r = reassigned_spectrogram(signal_frames, output_frame_size, magnitudes='power_db')
    assert np.all(X_r >= -120), 'min value: %f should be >= -120' % X_r.min()
    assert np.all(X_r <= 0), 'max value: %f should be <= 0' % X_r.max()

def test_reassigned_pitchgram_values_should_be_in_proper_range():
    frame_size = 4096
    hop_size = frame_size
    output_frame_size = 1024
    audio_file = os.path.join(DATA_DIR, 'she_brings_to_me.wav')
    signal_frames = SignalFrames(audio_file, frame_size, hop_size, mono_mix=True)
    X_r = pitchgram(signal_frames, output_frame_size, magnitudes='power_db')
    assert np.all(X_r >= -120), 'min value: %f should be >= -120' % X_r.min()
    assert np.all(X_r <= 0), 'max value: %f should be <= 0' % X_r.max()

# --- helper functions ---

def sample_time(since, until, fs=44100.):
    '''
    Generates time sample in given interval [since; until]
    with given sampling rate (fs).
    '''
    return np.arange(since, until, 1. / fs)

def sine(t, freq=1., amplitude=1., phase=0.):
    '''
    Samples the sine function given the time samples t,
    frequency (Hz), amplitude and phase [0; 2 * np.pi).
    '''
    return amplitude * np.sin(2 * np.pi * freq * t + phase)
<file_sep>from __future__ import absolute_import

from .reassignment import Spectrogram, LinearTransform, PitchTransform, \
    reassigned_spectrogram, pitchgram
from .signal import SignalFrames
from .sklearn import PitchgramTransformer
from .tuning import Tuning, PitchQuantizer
<file_sep>from __future__ import print_function, division

import numpy as np

from tfr import SignalFrames

def test_split_to_frames():
    signal_frames = SignalFrames(np.arange(23), frame_size=8, hop_size=6,
        sample_rate=44100)
    assert np.allclose(np.array([
        [ 0,  1,  2,  3,  4,  5,  6,  7],
        [ 6,  7,  8,  9, 10, 11, 12, 13],
        [12, 13, 14, 15, 16, 17, 18, 19],
        [18, 19, 20, 21, 22,  0,  0,  0],
    ]), signal_frames.frames)
<file_sep>import numpy as np
import scipy.misc
import sys


def save_raw_spectrogram_bitmap(file_name, spectrogram):
    # input:
    # rows = frequency bins (low to high)
    # columns = time
    # output:
    # rows = frequency bins (bottom to top)
    # columns = time (left to right)
    scipy.misc.imsave(file_name, spectrogram.T[::-1])

def spectrogram_to_image(npz_file):
  save_raw_spectrogram_bitmap(npz_file + '.png', np.load(npz_file)['arr_0'])

if __name__ == '__main__':
    spectrogram_to_image(sys.argv[1])
<file_sep>from __future__ import absolute_import

from sklearn.base import BaseEstimator, TransformerMixin

from .signal import SignalFrames
from .reassignment import pitchgram


class PitchgramTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, sample_rate=44100, frame_size=4096, hop_size=2048,
        output_frame_size=None,
        bin_range=[-48, 67], bin_division=1):
        self.sample_rate = sample_rate
        self.frame_size = frame_size
        self.hop_size = hop_size
        # if no output frame size is specified the input hop size is the default
        self.output_frame_size = output_frame_size if output_frame_size is not None else hop_size
        self.bin_range = bin_range
        self.bin_division = bin_division

    def transform(self, X, **transform_params):
        """
        Transforms audio clip X into a normalized pitchgram.
        Input: X - mono audio clip - numpy array of shape (samples,)
        Output: X_pitchgram - numpy array of shape (frames, bins)
        """
        signal_frames = SignalFrames(X, self.frame_size, self.hop_size,
            self.sample_rate, mono_mix=True)
        X_pitchgram = pitchgram(
            signal_frames,
            self.output_frame_size,
            magnitudes='power_db_normalized',
            bin_range=self.bin_range,
            bin_division=self.bin_division)
        return X_pitchgram

    def fit(self, X, y=None, **fit_params):
        return self
<file_sep>for spec_type in stft reassigned pitchgram; do
  python -m tfr.spectrogram_features \
    she_brings_to_me.wav she_brings_to_me_$spec_type.npz \
    -t $spec_type -b 4096 -p 2048 --output-frame-size=2048
done
<file_sep>from __future__ import print_function, division

import math
import numpy as np
import soundfile as sf

class SignalFrames():
    """
    Represents frames of time-domain signal of regular size with possible
    overlap plus it's metadata.

    The signal can be read from an numpy array a file via the soundfile library.
    The input array can be of shape `(samples,)` or `(samples, channels)`. By
    default the signal is mixed to mono to shape `(samples,)`. This can be
    disabled by specifying `mono_mix=False`.

    It is split into frames of `frame_size`. In case `hop_size < frame_size` the
    frames are overlapping. When the last frame is not fully covered by the
    signal it's padded with zeros.

    When reading signal from a file the sample rate can usually determined
    automatically, otherwise you shoud provide `sample_rate`.

    Attributes:
    - `frames` - signal split to frame, shape `(frames, frame_size)`
    - `frame_size`
    - `hop_size`
    - `length` - length of the source signal (in samples)
    - `duration` - duration of the source signal (in seconds)
    - `start_times` - array of start times of each frame (in seconds)

    Example usage:

    ```
    signal_frames = SignalFrames('audio.flac', frame_size=4096, hop_size=1024)
    spectrogram = np.fft.fft(signal_frames.frames * window)
    ```

    :param source: source of the time-domain signal - numpy array, file name,
    file-like object
    :param frame_size: size of each frame (in samples)
    :param hop_size: hop between frame starts (in samples)
    :param sample_rate: sample rate (required when source is an array)
    :mono_mix: indicates that multi-channel signal should be mixed to mono
    (mean of all channels)
    """
    def __init__(self, source, frame_size=4096, hop_size=2048, sample_rate=None,
        mono_mix=True):
        if type(source) == np.ndarray:
            signal = source
            self.sample_rate = sample_rate
        else:
            signal, self.sample_rate = sf.read(source)

        if mono_mix:
            signal = self._to_mono(signal)

        self.frames = self._split_to_frames(signal, frame_size, hop_size)
        self.frame_size = frame_size
        self.hop_size = hop_size
        self.length = len(signal)
        self.duration = self.length / self.sample_rate
        self.start_times = np.arange(0, self.duration, self.hop_size / self.sample_rate)

    def _split_to_frames(self, x, frame_size, hop_size):
        """
        Splits the input audio signal to frame of given size (in samples).
        Start position of each frame is determined by given hop size.
        The last frame is right zero-padded if needed.
        input:
        x - array-like representing the audio signal
        """
        if hop_size is None:
            hop_size = frame_size
        frame_count = int(math.ceil(len(x) / hop_size))
        def pad(x, size, value=0):
            padding_size = size - len(x)
            if padding_size:
                x = np.pad(x, (0, padding_size), 'constant', constant_values=(0, 0))
            return x
        frames = np.vstack(
            pad(x[start:start + frame_size], frame_size) \
            for start in np.arange(0, hop_size * frame_count, hop_size))
        return frames

    def _to_mono(self, samples):
        if samples.ndim == 1:
            return samples
        else:
            return samples.mean(axis=-1)

def mean_power(x_frames):
    return np.sqrt(np.mean(x_frames**2, axis=-1))

def power(x_frames):
    return np.sqrt(np.sum(x_frames**2, axis=-1))

def mean_energy(x_frames):
    """
    Example usage:

    import matplotlib.pyplot as plt
    import soundfile as sf
    from analysis import read_frames

    def analyze_mean_energy(file, frame_size=1024):
        frames, t, fs = read_frames(x, frame_size)
        y = mean_energy(frames)
        plt.semilogy(t, y)
        plt.ylim(0, 1)
    """
    return np.mean(x_frames**2, axis=-1)

def energy(x_frames):
    return np.sum(x_frames**2, axis=-1)
<file_sep>"""
The goal is to transform an audio signal into an STFT spectrogram in a form
suitable as features for machine learning.
"""

from __future__ import print_function, division

import numpy as np
import os

from .signal import SignalFrames
from .reassignment import reassigned_spectrogram, pitchgram


def spectrogram_features(file_name, frame_size, output_frame_size, hop_size, spectrogram_type, magnitudes='power_db_normalized'):
    signal_frames = SignalFrames(file_name, frame_size, hop_size, mono_mix=True)

    if spectrogram_type == 'stft':
        X = reassigned_spectrogram(signal_frames, output_frame_size,
            magnitudes=magnitudes, reassign_time=False, reassign_frequency=False)
    elif spectrogram_type == 'reassigned':
        X = reassigned_spectrogram(signal_frames, output_frame_size,
            magnitudes=magnitudes)
    elif spectrogram_type == 'pitchgram':
        X = pitchgram(signal_frames, output_frame_size, magnitudes=magnitudes)
    else:
        raise ValueError('unknown spectrogram type: %s' % spectrogram_type)

    return X

def spectrogram_features_to_file(input_filename, output_filename, frame_size, output_frame_size, hop_size, spectrogram_type, magnitudes='power_db'):
    X = spectrogram_features(input_filename, frame_size, output_frame_size, hop_size, spectrogram_type, magnitudes)
    np.savez_compressed(output_filename, X)
    # scipy.misc.imsave(output_filename.replace('.npz', '.png'), X.T[::-1])

def default_output_filename(input_file_name, type):
    return os.path.basename(input_file_name).replace('.wav', '_power_spectrogram_%s.npz' % type)

def parse_args():
    import argparse

    parser = argparse.ArgumentParser(description='Extracts STFT magnitude spectrogram features.')
    parser.add_argument('input_file', metavar='INPUT', help='input file in WAV format')
    parser.add_argument('output_file', metavar='OUTPUT', nargs='?', help='output file in NumPy npz format')
    parser.add_argument('-b', '--frame-size', type=int, default=2048, help='STFT frame size')
    parser.add_argument('-p', '--hop-size', type=int, default=512, help='STFT hop size')
    parser.add_argument('-o', '--output-frame-size', type=int, default=512, help='output frame size')
    parser.add_argument('-t', '--type', default='stft', help='plain "stft", "reassigned" spectrogram or "pitchgram"')
    parser.add_argument('-m', '--magnitudes', default='power_db_normalized',
        choices=['linear', 'power', 'power_db', 'power_db_normalized'])

    return parser.parse_args()

def main():
    args = parse_args()

    output = args.output_file if args.output_file else default_output_filename(args.input_file, args.type)

    spectrogram_features_to_file(args.input_file, output, args.frame_size,
        args.output_frame_size, args.hop_size, args.type, args.magnitudes)

if __name__ == '__main__':
    main()
<file_sep># tfr - time-frequency reassignment in Python

[![PyPI version](https://img.shields.io/pypi/v/tfr.svg)](https://pypi.python.org/pypi/tfr)
![Supported Python versions](https://img.shields.io/pypi/pyversions/tfr.svg)
![License](https://img.shields.io/pypi/l/tfr.svg)

Spectral audio feature extraction using [time-frequency reassignment](https://en.wikipedia.org/wiki/Reassignment_method).

<img src="examples/multicomponent_animation.gif" alt="reassigned spectrogram illustration" width="516" height="436">

Besides normal spectrograms it allows to compute reassigned spectrograms, transform them (eg. to log-frequency scale) and requantize them (eg. to musical pitch bins). This is useful to obtain good features for audio analysis or machine learning on audio data.

A reassigned spectrogram often provides more precise localization of energy in the time-frequency plane than a plain spectrogram. Roughly said in the reassignment method we use the phase (which is normally discarded) and move the samples on the time-frequency plane to a more suitable place computed from derivatives of the phase.

This library supports reassignment in both frequency and time (both are optional). As well it does requantization from the input overlapping grid to an non-overlapping output grid.

It is a good building block to compute [chromagram features](https://en.wikipedia.org/wiki/Chroma_feature) (aka pitch class profiles) where pitch is transformed into pitch class by ignoring the octave. See also [harmonic pitch class profiles](https://en.wikipedia.org/wiki/Harmonic_pitch_class_profiles).

## Installation

```
pip install tfr
```

Or for development (all code changes will be available):

```
git clone https://github.com/bzamecnik/tfr.git
pip install -e tfr
```

## Usage

### Split audio signal to frames

You can read time-domain signal from an audio file (using the `soundfile` library) and split it into frames for spectral processing.

```
import tfr
signal_frames = tfr.SignalFrames('audio.flac')
```

`SignalFrames` instance contains the signal split into frames and some metadata useful for further processing.

The signal values are normalized to [0.0, 1.0] and the channels are converted to mono.

It is possible to provide the signal a numpy array as well.

```
import tfr
x = np.sin(2 * np.pi * 10 * np.linspace(0, 1, 1000))
signal_frames = tfr.SignalFrames(x)
```

### Minimal example - pitchgram from audio file

```
import tfr
x_pitchgram = tfr.pitchgram(tfr.SignalFrames('audio.flac'))
```

From audio frames it computes a reassigned pitchgram of shape `(frame_count, bin_count)` with values being log-magnitudes in dBFS `[-120.0, 0.0]`. Sensible parameters are used by default, but you can change them if you wish.

### Reassigned spectrogram

Like normal one but sharper and requantized.

```
import tfr
x_spectrogram = tfr.reassigned_spectrogram(tfr.SignalFrames('audio.flac'))
```

### Signal frames with specific parameters

- `frame_size` - affects the FFT size - trade-off between frequency and time resolution, good to use powers of two, eg. 4096
- `hop_size` - affects the overlap between frames since a window edges fall to zero, eg. half of frame_size (2048)

```
import tfr
signal_frames = tfr.SignalFrames('audio.flac', frame_size=1024, hop_size=256)
```

### General spectrogram API

The `pitchgram` and `reassigned_spectrogram` functions are just syntax sugar for the `Spectrogram` class. You can use it directly to gain more control.

General usage:

```
x_spectrogram = tfr.Spectrogram(signal_frames).reassigned()
```

From one Spectrogram instance you can efficiently compute reassigned spectrograms with various parameters.

```
s = tfr.Spectrogram(signal_frames)
x_spectrogram_tf = s.reassigned(output_frame_size=4096)
x_spectrogram_f = s.reassigned(output_frame_size=512)
```

Different window function (by default we use Hann window):

```
import scipy
x_spectrogram = tfr.Spectrogram(signal_frames, window=scipy.blackman).reassigned()
```

Different output frame size (by default we make it the same as input hop size):

```
x_spectrogram = tfr.Spectrogram(signal_frames).reassigned(output_frame_size=512)
```

Disable reassignment of time and frequency separately:

```
s = tfr.Spectrogram(signal_frames)
x_spectrogram = s.reassigned(reassign_time=False, reassign_frequency=False)
x_spectrogram_t = s.reassigned(reassign_frequency=False)
x_spectrogram_f = s.reassigned(reassign_time=False)
x_spectrogram_tf = s.reassigned()
```

Disable decibel transform of output values:

```
x_spectrogram = tfr.Spectrogram(signal_frames).reassigned(magnitudes='power')
```

Magnitudes in the spectrogram can be transformed at the end in multiple ways given by the `magnitudes` parameter:

- `linear` - energy spectrum
- `power` - power spectrum
- `power_db` - power spectrum in decibels, range: [-120, 0]
- `power_db_normalized` - power spectrum in decibels normalized to range: [0, 1]
  - this is useful as a feature

Use some specific transformation of the output values. `LinearTransform` (default) is just for normal spectrogram, `PitchTransform` is for pitchgram. Or you can write your own.

```
x_spectrogram = tfr.Spectrogram(signal_frames).reassigned(transform=LinearTransform())
```

```
x_pitchgram = tfr.Spectrogram(signal_frames).reassigned(transform=PitchTransform())
```

```
class LogTransform():
  def __init__(self, bin_count=100)
    self.bin_count = bin_count

  def transform_freqs(self, X_inst_freqs, sample_rate):
      X_y = np.log10(np.maximum(sample_rate * X_inst_freqs, eps))
      bin_range = (0, np.log10(sample_rate))
      return X_y, self.bin_count, bin_range

x_log_spectrogram = tfr.Spectrogram(signal_frames).reassigned(transform=LogTransform())
```

### Pitchgram parameters

In pitchgram the frequencies are transformed into pitches in some tuning and then quantized to bins. You can specify the tuning range of pitch bins and their subdivision.

- `tuning` - instance of `Tuning` class, transforms between pitch and frequency
- `bin_range` is in pitches where 0 = 440 Hz (A4), 12 is A5, -12 is A3, etc.
- `bin_division` - bins per each pitch

### Extract features via CLI

```
# basic STFT spectrogram
python -m tfr.spectrogram_features audio.flac spectrogram.npz
# reassigned STFT spectrogram
python -m tfr.spectrogram_features audio.flac -t reassigned reassigned_spectrogram.npz
# reassigned pitchgram
python -m tfr.spectrogram_features audio.flac -t pitchgram pitchgram.npz
```

Look for other options:

```
python -m tfr.spectrogram_features --help
```

### scikit-learn transformer

In order to extract pitchgram features within a sklearn pipeline, we can use `PitchgramTransformer`:

```
import soundfile as sf
x, fs = sf.read('audio.flac')

from tfr.signal import to_mono
from tfr.sklearn import PitchgramTransformer
ct = PitchgramTransformer(sample_rate=fs)
x_pitchgram = ct.transform(x)

# output:
#  - shape: (frame_count, bin_count)
#   - values in dBFB normalized to [0.0, 1.0]
```

## Status

Currently it's alpha. I'm happy to extract it from some other project into a separate repo and package it. However, the API must be completely redone to be more practical and obvious.

## About

- Author: <NAME> ([@bzamecnik](http://twitter.com/bzamecnik))
- License: MIT

### Support the project

Need some consulting or coding work regarding audio processing, machine learning or big data? Drop me a message via [email](mailto:<EMAIL>?subject=Work+inquiry+-+based+on+tfr) or [LinkedIn](https://www.linkedin.com/in/bohumirzamecnik). Or just say hello :).

## Literature

- [A Unified Theory of Time-Frequency Reassignment](https://arxiv.org/abs/0903.3080) - <NAME>, <NAME>, Digital Signal Processing 30 September 2005
- [Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications](http://acousticslab.org/learnmoresra/files/fulopfitz2006jasa119.pdf) - <NAME>, <NAME>, Journal of Acoustical Society of America, Jan 2006
- [Time Frequency Reassignment: A Review and Analysis](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.4.1053&rep=rep1&type=pdf) - <NAME>, <NAME>,
Technical Report, Cambridge University Engineering Dept.
- [Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method](http://perso.ens-lyon.fr/patrick.flandrin/IEEE_SP1995.pdf) - <NAME>, <NAME>, IEEE Transactions on Signal Processing, vol. 43, no. 5, May 1995
- [Time–frequency reassignment: from principles to algorithms](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.331.5416&rep=rep1&type=pdf) - <NAME>, <NAME>, <NAME>, CRC Press 2003
- [Time-frequency toolbox for Matlab, user’s guide and reference guide](http://iut-saint-nazaire.univ-nantes.fr/~auger/tftb.html) - F.Auger, P.Flandrin, P.Goncalves, O.Lemoine
<file_sep>"""
An example showing analysis of multicomponent sound with and without time or
frequency reassignment. With both TF-reassignment the resulting spectrogram is
much sharper altough it still cannot resolve inseparable places.
"""

import numpy as np
import soundfile as sf
import subprocess
import tfr

def sample_time(since, until, fs=44100.):
    """
    Generates time sample in given interval [since; until]
    with given sampling rate (fs).
    """
    return np.arange(since, until, 1. / fs)

def freq_mod_sine(t, carrier_freq, mod_freq, mod_amp):
    return np.sin(2 * np.pi * carrier_freq * t
        + mod_amp * np.sin(2 * np.pi * mod_freq * t))

def sinusoid(t, freq):
    """
    t - array of time samples
    freq - frequency - constant or array of same size as t
    """
    return np.sin(2 * np.pi * freq * t)

def linear_chirp(t, start_freq, end_freq):
    slope = (end_freq - start_freq) / (t[-1] - t[0])
    return np.sin(2 * np.pi * (start_freq * t + 0.5 * slope * t**2))

def generate_example_sound(fs=44100, duration=3, carrier_freq=2000, mod_freq=1, mod_amp=1000, click_freq=10):
    t = sample_time(0, duration, fs)

    component_count = 4
    amplitude = 1 / component_count

    # FM component
    x = amplitude * freq_mod_sine(t, carrier_freq, mod_freq, mod_amp)
    # constant tone component
    x += amplitude * sinusoid(t, carrier_freq)
    # linear chirps
    x += amplitude * linear_chirp(t, carrier_freq-mod_amp, carrier_freq+mod_amp)
    x += amplitude * linear_chirp(t, carrier_freq+mod_amp, carrier_freq-mod_amp)
    # sawtooth pulse-train component
    for i in range(duration*click_freq):
        idx = int(i * fs / click_freq)
        x[idx:idx+10] = 1
    sf.write('multicomponent.wav', x, fs)
    return x, fs

def compute_example_spectrograms():
    x, fs = generate_example_sound()
    signal_frames = tfr.SignalFrames(x, sample_rate=fs, frame_size=4096, hop_size=512)

    def band_pass(X, y_range):
        return X[:, slice(*y_range)]

    y_range = (84, 287) # 900-3100 Hz

    # plain non-reassigned spectrogram
    X_stft = band_pass(
        tfr.reassigned_spectrogram(
            signal_frames,
            reassign_time=False,
            reassign_frequency=False),
        y_range)
    tfr.plots.save_raw_spectrogram_bitmap('multicomponent_stft.png', X_stft)

    # time-reassigned spectrogram
    X_t = band_pass(
        tfr.reassigned_spectrogram(signal_frames, reassign_frequency=False),
        y_range)
    tfr.plots.save_raw_spectrogram_bitmap('multicomponent_t.png', X_t)

    # frequency-reassigned spectrogram
    X_f = band_pass(
        tfr.reassigned_spectrogram(signal_frames, reassign_time=False),
        y_range)
    tfr.plots.save_raw_spectrogram_bitmap('multicomponent_f.png', X_f)

    # time-frequency reassigned spectrogram
    X_tf = band_pass(tfr.reassigned_spectrogram(signal_frames), y_range)
    tfr.plots.save_raw_spectrogram_bitmap('multicomponent_tf.png', X_tf)


def make_animation():

    def add_label_and_resize(source_file, dest_file, label):
        subprocess.call([
            'convert',
            source_file,
            # '-filter', 'point', '-resize', '200%',
            # '-pointsize', '30',
            'label:%s' % label, '+swap', '-gravity', 'Center', '-append',
            dest_file])

    for name, label in [('stft', 'no'), ('f', 'frequency'), ('t', 'time'), ('tf', 'time-frequency')]:
        add_label_and_resize(
            'multicomponent_%s.png' % name,
            'multicomponent_%s_label.png' % name,
            '%s reassignment' % label)

    subprocess.call([
        'convert',
        '-delay', '100',
        '-loop', '0',
        'multicomponent_stft_label.png',
        'multicomponent_f_label.png',
        'multicomponent_tf_label.png',
        'multicomponent_stft_label.png',
        'multicomponent_t_label.png',
        'multicomponent_tf_label.png',
        'multicomponent_animation.gif'
    ])

if __name__ == '__main__':
    compute_example_spectrograms()
    make_animation()
<file_sep>from __future__ import print_function, division

from tfr import Tuning


def test_pitch_to_relative_freq():
    tuning_step1 = Tuning(steps_per_octave=1)
    tuning_step12 = Tuning(steps_per_octave=12)
    assert 1. == tuning_step1.pitch_to_relative_freq(0.)
    assert 2. == tuning_step1.pitch_to_relative_freq(1.)
    assert 4. == tuning_step1.pitch_to_relative_freq(2.)
    assert 0.5 == tuning_step1.pitch_to_relative_freq(-1.)

    assert 1. == tuning_step12.pitch_to_relative_freq(0.)
    assert 2. == tuning_step12.pitch_to_relative_freq(12.)
    assert 4. == tuning_step12.pitch_to_relative_freq(24.)
    assert 0.5 == tuning_step12.pitch_to_relative_freq(-12.)

def test_pitch_to_freq():
    tuning = Tuning()
    assert 440. == tuning.pitch_to_freq(0.)
    assert 880. == tuning.pitch_to_freq(12.)
    assert 1760. == tuning.pitch_to_freq(24.)
    assert 220. == tuning.pitch_to_freq(-12.)

    assert abs(466.1637615180899 - tuning.pitch_to_freq(1.)) < 1e-10
    assert abs(415.3046975799451 - tuning.pitch_to_freq(-1.)) < 1e-10
    assert abs(1318.5102276514797 - tuning.pitch_to_freq(12 + 7)) < 1e-10

# TODO: test:
# - freq_to_pitch()
# - PitchQuantizer
# - various configurations of Tuning
<file_sep>"""
A regression test for computing three kinds of spectrogram.

Just to ensure we didn't break anything.
"""

from __future__ import print_function, division

import numpy as np
import os
from tfr.spectrogram_features import spectrogram_features

DATA_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')

def test_spectrograms():
    for spectrogram_type in ['stft', 'reassigned', 'pitchgram']:
        yield assert_spectrogram_is_ok, spectrogram_type

def assert_spectrogram_is_ok(spectrogram_type):
    audio_file = os.path.join(DATA_DIR, 'she_brings_to_me.wav')
    X = spectrogram_features(audio_file, frame_size=4096, hop_size=2048,
        output_frame_size=2048,
        spectrogram_type=spectrogram_type, magnitudes='power_db_normalized')
    npz_file = os.path.join(DATA_DIR, 'she_brings_to_me_%s.npz' % spectrogram_type)
    X_expected = np.load(npz_file)['arr_0']
    assert X.shape == X_expected.shape, (X.shape, X_expected.shape)
    print('spectrogram [%s]: max abs error' % spectrogram_type, abs(X - X_expected).max())
    assert np.allclose(X, X_expected)
<file_sep>from __future__ import print_function, division

import numpy as np
import os
from tfr import SignalFrames
from tfr.signal import energy, mean_power
from tfr.spectrogram import create_window, stft_spectrogram

DATA_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')

def test_window_should_be_normalized():
    def assert_ok(size):
        w = create_window(size)
        assert np.allclose(energy(w), len(w))
        assert np.allclose(mean_power(w), 1.0)

    for size in [16, 100, 512, 777, 4096]:
        yield assert_ok, size

def test_spectrogram_db_magnituds_should_be_in_proper_range():
    frame_size = 4096
    hop_size = 4096
    audio_file = os.path.join(DATA_DIR, 'she_brings_to_me.wav')
    signal_frames = SignalFrames(audio_file, frame_size, hop_size, mono_mix=True)
    w = create_window(frame_size)
    X = stft_spectrogram(signal_frames.frames, w, magnitudes='power_db')
    assert np.all(X >= -120), 'min value: %f should be >= -120' % X.min()
    assert np.all(X <= 0), 'max value: %f should be <= 0' % X.max()
