Source code for niacin.timeseries.freq

#!/usr/bin/env python
# -*- coding: utf-8

"""Frequency domain transforms
"""

import numpy as np
from numpy import random
from scipy.fft import irfft, rfft


[docs]def add_discrete_phase_shifts(x: np.ndarray, p: float=0.01, m: float=0.1): r"""Shift each frequency component with probability p by distance m*len(x). Compute real FFT on sequence. Then, for each frequency entry, with probability p, swap it with an entry that is ±m*len(x)//2 steps away. Swapped entries are tagged so that they are not moved twice. The 0th component is ignored. Args: x: sequence p: per-frequency probability of shifting m: magnitude of shift Returns: enriched sequence Examples: >>> x = np.sin(np.linspace(0, 6*np.pi, 100)) | +1 -\\ -\\ --\ | / \\ / \\ / \ | / \ / \ / \ | / \ / \ / \ |/ \ / \ / \ |/ \ / \ / \ +--------\--------/--------\---------/-------\---------+-- +0 \ / \ / \ +99 | \ / \ / \ / | \ / \ / \ / | \ / \ / \ / | \ / \\ / \\ / -1 --/ \-/ \-/ | >>> ts.add_discrete_phase_shifts(x, 1.0, 0.1) | +1.03587-| \ -| \ \| \ | |/| | | /| || /| / | | | /| ||| | | || || | | | | | | | | || | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | +---|--|--|---|--|---|--|---|--|--|---|--|---|--|--|---+-- +0 | | | | | | | | | | | | | | +99 | | | | | | | | | | | | | | | | * | | | | | | | | | | | | | | | | | | | | | | | | | | | | || | | || | \| | | \ | | | \| || \| |/ -1.01053 / /| / | -| | -/ | """ fx = rfft(x) step = round(len(fx) * m) ps = random.binomial(1, p, len(fx)) signs = random.binomial(1, 0.5, len(fx)) * 2 - 1 # skip 0 Hz component for i in range(1, len(fx)): if ps[i] == 1: j = min(max(i + signs[i] * step, 0), len(fx)-1) if ps[j] == 1: fx[i], fx[j] = fx[j], fx[i] # don't shift the same component more than once ps[i], ps[j] = -1, -1 return irfft(fx)
[docs]def add_random_frequency_noise(x: np.ndarray, p: float=0.01, m: float=0.1): r"""Add gaussian noise to each frequency with probability p and magnitude N(0,1) * m * max(fft(x)) Args: x: sequence p: per-frequency probability of noise m: magnitude of noise Returns: enriched sequence Examples: >>> x = np.sin(np.linspace(0, 6*np.pi, 100)) | +1 -\\ -\\ --\ | / \\ / \\ / \ | / \ / \ / \ | / \ / \ / \ |/ \ / \ / \ |/ \ / \ / \ +--------\--------/--------\---------/-------\---------+-- +0 \ / \ / \ +99 | \ / \ / \ / | \ / \ / \ / | \ / \ / \ / | \ / \\ / \\ / -1 --/ \-/ \-/ | >>> ts.add_random_frequency_noise(x, 0.5, 0.1) | +1.59613 \| | -\/\ | || |\| / \ -||| \/||\/| |/| | |/ | \ |/ -/-\ | \ | | | ||\| | | | +-------|-|-|--|-|-\\|--\|||-|-------|-------|-\------|+-- +0 ||| | | |/ /| | | | | | /| +99* | ||| | || / | -/| / | | | |\ |\| |\\-\| | | | | / \| //|| | || \\| | | || / || | | | | -1.90176 | | """ fx = rfft(x) ps = random.binomial(1, p, len(fx)) noise = random.randn(2*len(fx)).view(np.complex128) * m * fx.max() fx[1:] = fx[1:] + ps[1:] * noise[1:] return irfft(fx)
[docs]def add_high_frequency_noise(x: np.ndarray, p: float=0.01, m: float=0.1): r"""Add gaussian noise to single highest frequency component with probability p and magnitude N(0,1) * m * max(fft(x)) Args: x: sequence p: per-sequence probability of noise m: magnitude of noise Returns: enriched sequence Examples: >>> x = np.sin(np.linspace(0, 6*np.pi, 100)) | +1 -\\ -\\ --\ | / \\ / \\ / \ | / \ / \ / \ | / \ / \ / \ |/ \ / \ / \ |/ \ / \ / \ +--------\--------/--------\---------/-------\---------+-- +0 \ / \ / \ +99 | \ / \ / \ / | \ / \ / \ / | \ / \ / \ / | \ / \\ / \\ / -1 --/ \-/ \-/ | >>> ts.add_high_frequency_noise(x, 1.0, 1.0) | +2.11065 |||| || || | |||||| |||||| ||||||| ||||||||| |||||||| ||||||||| ||||||||| | |||||||||| ||||||||||| * |||||||||||| |||||||||||| ||||||||||||| | ||||||||||||| ||||||||||||||| | ||||||||||||||| | || +||||||||||||||||||||||||||||||||||||||||||||||||||||||+-- +0| | ||||||||||||||| | ||||||||||||||| |||||||||+99 || ||||||||||||| |||||||||||| ||||||||||| | ||||||||||| ||||||||| | |||||||||| | ||||||||| |||||||| |||||||| | ||||||| |||||| |||||| -2.11065 || || |||| |||| | """ fx = rfft(x) if random.binomial(1, p): noise = random.randn(2).view(np.complex128) * m * fx.max() fx[-1] = fx[-1] + noise return irfft(fx)
[docs]def remove_random_frequency(x: np.ndarray, p: float=0.01, m=None): r"""Remove each frequency component with probability p. Args: x: sequence p: per-frequency probability of removal m: ignored Returns: enriched sequence Examples: >>> x = np.sin(np.linspace(0, 6*np.pi, 100)) | +1 -\\ -\\ --\ | / \\ / \\ / \ | / \ / \ / \ | / \ / \ / \ |/ \ / \ / \ |/ \ / \ / \ +--------\--------/--------\---------/-------\---------+-- +0 \ / \ / \ +99 | \ / \ / \ / | \ / \ / \ / | \ / \ / \ / | \ / \\ / \\ / -1 --/ \-/ \-/ | >>> ts.reverse(x, 1.0) | +1 | | | | | +------------------------------------------------------+-- +0 +99 | | | | -1 """ fx = rfft(x) ps = random.binomial(1, 1-p, len(fx)) fx = fx * ps return irfft(fx)