Signal Processing

from scipy import signal
import numpy as np

# Common direct imports
from scipy.signal import (
    butter, cheby1, cheby2, ellip, bessel,
    filtfilt, lfilter, sosfilt, sosfiltfilt,
    freqz, freqs, sosfreqz,
    find_peaks, peak_prominences, peak_widths,
    correlate, convolve, fftconvolve,
    spectrogram, stft, istft, welch, periodogram, lombscargle,
    resample, resample_poly, decimate,
    hilbert, chirp, gausspulse, sweep_poly,
    tf2sos, tf2zpk, zpk2sos, zpk2tf, sos2tf,
    iirfilter, firwin, firwin2, remez,
    savgol_filter, medfilt, wiener
)

from scipy.signal import butter, cheby1, cheby2, ellip, bessel

fs = 1000.0    # sample rate (Hz)
fc = 100.0     # cutoff frequency (Hz)

# Butterworth — maximally flat passband
b, a = butter(N=4, Wn=fc, btype='low', analog=False, fs=fs)
# btype: 'low', 'high', 'band', 'bandstop'

# Band-pass
b, a = butter(N=4, Wn=[50, 150], btype='band', fs=fs)

# Chebyshev Type I — equiripple passband
b, a = cheby1(N=4, rp=1, Wn=fc, btype='low', fs=fs)   # rp = max ripple (dB)

# Chebyshev Type II — equiripple stopband
b, a = cheby2(N=4, rs=40, Wn=fc, btype='low', fs=fs)  # rs = min attenuation (dB)

# Elliptic — equiripple in both bands (steepest rolloff)
b, a = ellip(N=4, rp=1, rs=60, Wn=fc, btype='low', fs=fs)

# Bessel — maximally flat group delay (linear phase)
b, a = bessel(N=4, Wn=fc, btype='low', norm='phase', fs=fs)
# norm: 'phase', 'delay', 'mag'

from scipy.signal import iirfilter, iirdesign

b, a = iirfilter(N=4, Wn=100, btype='low', ftype='butter', fs=1000, output='ba')
sos = iirfilter(N=4, Wn=100, btype='low', ftype='butter', fs=1000, output='sos')

# Design from specs (passband/stopband frequencies and attenuations)
b, a = iirdesign(wp=0.2, ws=0.4, gpass=1, gstop=40, ftype='ellip')
# wp, ws as fractions of Nyquist [0,1] when fs not given

from scipy.signal import firwin, firwin2, remez, minimum_phase

# firwin — windowed sinc design
h = firwin(numtaps=101, cutoff=100.0, window='hamming', fs=1000.0, pass_zero='low')
h = firwin(numtaps=101, cutoff=[50, 150], window='hamming', fs=1000.0, pass_zero=False)
# pass_zero: True/'low'=lowpass, False/'band'=bandpass, 'high'=highpass, 'bandstop'=bandstop

# Common windows
# 'hamming', 'hann', 'blackman', 'kaiser', ('kaiser', beta), 'flattop'

# Kaiser window — control sidelobe vs transition width
# As rule of thumb: beta~2.285*(N-1)*width/fs + 0.48*atten
h = firwin(numtaps=201, cutoff=100.0, window=('kaiser', 8.0), fs=1000.0)

# firwin2 — arbitrary frequency response
freqs_spec = [0, 100, 200, 300, 500]
gains_spec = [0, 1, 1, 0, 0]
h = firwin2(numtaps=101, freq=freqs_spec, gain=gains_spec, fs=1000.0)

# Parks-McClellan (equiripple, optimal)
h = remez(numtaps=101,
           bands=[0, 80, 100, 500],
           desired=[0, 1],            # 0 in [0,80], 1 in [100,500]
           fs=1000.0)

from scipy.signal import butter, tf2sos, tf2zpk, zpk2sos, sos2tf

# Get SOS directly (recommended)
sos = butter(N=8, Wn=100, btype='low', fs=1000, output='sos')

# Convert between formats
b, a = butter(8, 100, btype='low', fs=1000, output='ba')
sos  = tf2sos(b, a)         # ba → sos
z, p, k = tf2zpk(b, a)     # ba → zpk
sos  = zpk2sos(z, p, k)    # zpk → sos
b, a = sos2tf(sos)          # sos → ba

from scipy.signal import lfilter, filtfilt, sosfilt, sosfiltfilt

t = np.linspace(0, 1, 1000)
x = np.sin(2*np.pi*50*t) + np.random.randn(1000)*0.5

sos = butter(4, 100, btype='low', fs=1000, output='sos')

# Zero-phase filtering (forward-backward) — no phase shift, doubles effective order
y = sosfiltfilt(sos, x)                     # recommended

# Causal filtering (forward only) — introduces phase delay
y = sosfilt(sos, x)

# Initial conditions for sosfilt (useful for streaming)
zi = signal.sosfilt_zi(sos) * x[0]         # initial conditions
y, zf = sosfilt(sos, x, zi=zi)

# Legacy BA-format versions
b, a = butter(4, 100, btype='low', fs=1000, output='ba')
y = filtfilt(b, a, x)     # zero-phase
y = lfilter(b, a, x)      # causal

# Savitzky-Golay (polynomial smoothing)
from scipy.signal import savgol_filter
y_smooth = savgol_filter(x, window_length=51, polyorder=3, deriv=0)
y_deriv  = savgol_filter(x, window_length=51, polyorder=3, deriv=1, delta=1/1000)

from scipy.signal import freqz, freqs, sosfreqz

# Digital filter (ba)
w, h = freqz(b, a, worN=512)    # w in rad/sample
freq = w * fs / (2 * np.pi)     # convert to Hz
import matplotlib.pyplot as plt
plt.plot(freq, 20*np.log10(abs(h)))

# Digital filter (sos)
w, h = sosfreqz(sos, worN=512, fs=fs)   # fs given → w in Hz

# Analog filter
w, h = freqs(b, a, worN=np.logspace(0, 4, 200))

# Group delay
w, gd = signal.group_delay((b, a))

from scipy.signal import welch, periodogram, spectrogram, stft, istft

fs = 1000.0
x = np.random.randn(10000)

# Welch's power spectral density (averaged)
f, Pxx = welch(x, fs=fs, nperseg=256, noverlap=128, window='hann',
                scaling='density')   # 'density' (V²/Hz) or 'spectrum' (V²)

# Periodogram (single FFT, noisy)
f, Pxx = periodogram(x, fs=fs, window='hann', scaling='density')

# Spectrogram (time-frequency)
f, t_spec, Sxx = spectrogram(x, fs=fs, nperseg=256, noverlap=200,
                               window='hann', scaling='density')
# Sxx.shape = (freq_bins, time_bins)

# STFT / ISTFT
f, t_stft, Zxx = stft(x, fs=fs, nperseg=256, noverlap=128, window='hann')
x_rec, t_rec = istft(Zxx, fs=fs, nperseg=256, noverlap=128, window='hann')

# Lomb-Scargle (non-uniform time sampling)
from scipy.signal import lombscargle
t_uneven = np.sort(np.random.rand(200)) * 10
x_uneven = np.sin(2*np.pi*2*t_uneven)
f_ls = np.linspace(0.1, 10, 500)
omega = 2 * np.pi * f_ls
pgram = lombscargle(t_uneven, x_uneven, omega, normalize=True)

from scipy.signal import find_peaks, peak_prominences, peak_widths

x = np.sin(np.linspace(0, 10, 500)) + np.random.randn(500)*0.1

# Basic peak finding
peaks, props = find_peaks(x)

# With constraints
peaks, props = find_peaks(x,
    height=0.5,              # minimum peak height (scalar or [min, max])
    threshold=0.1,           # min vertical distance to neighbors
    distance=20,             # min samples between peaks
    prominence=0.3,          # min prominence
    width=5,                 # min peak width
    wlen=100,                # window for prominence calculation
    rel_height=0.5)          # fractional height for width measurement

# Returned properties dict (only requested keys present)
print(props.keys())  # 'peak_heights', 'prominences', 'widths', etc.

# Separate computations
prominences, l_bases, r_bases = peak_prominences(x, peaks)
widths, width_heights, l_ips, r_ips = peak_widths(x, peaks, rel_height=0.5)

# Find valleys (negate signal)
valleys, _ = find_peaks(-x, prominence=0.3)

from scipy.signal import convolve, correlate, fftconvolve, oaconvolve

a = np.array([1, 2, 3, 4, 5], dtype=float)
b = np.array([1, 0, -1], dtype=float)

# 1-D convolution
c = convolve(a, b, mode='full')     # len = len(a)+len(b)-1
c = convolve(a, b, mode='same')     # len = max(len(a), len(b))
c = convolve(a, b, mode='valid')    # len = max(len(a), len(b)) - min(...) + 1

# Cross-correlation (no flip of b)
r = correlate(a, b, mode='full')

# FFT-based (faster for large arrays)
c = fftconvolve(a, b, mode='same')

# Overlap-add (best for very different-length arrays)
c = oaconvolve(a, b, mode='same')

# 2-D
from scipy.signal import convolve2d, correlate2d
img = np.random.rand(100, 100)
kernel = np.ones((5, 5)) / 25
smoothed = convolve2d(img, kernel, mode='same', boundary='symm')

from scipy.signal import resample, resample_poly, decimate

x = np.sin(np.linspace(0, 2*np.pi, 1000))

# Resample to N samples (FFT-based)
x_resampled = resample(x, 250)    # downsample 4x

# Resample by rational factor (polyphase, anti-aliased)
x_up = resample_poly(x, up=3, down=2)          # 3/2 rate change
x_up = resample_poly(x, 3, 2, window=('kaiser', 5.0))

# Decimate (integer downsampling, applies anti-alias filter)
x_down = decimate(x, q=4, ftype='iir', zero_phase=True)  # q must be <= 13 for iir
x_down = decimate(x, q=10, ftype='fir')

from scipy.signal import chirp, gausspulse, sweep_poly, unit_impulse

t = np.linspace(0, 1, 1000)

# Linear chirp (frequency sweep)
x = chirp(t, f0=10, f1=100, t1=1, method='linear')
x = chirp(t, f0=1, f1=100, t1=1, method='logarithmic')
x = chirp(t, f0=1, f1=100, t1=1, method='quadratic')

# Gaussian pulse
x, envelope = gausspulse(t - 0.5, fc=50, bw=0.5, retenv=True)

# Unit impulse
imp = unit_impulse(100, idx='mid')      # impulse at center
imp = unit_impulse(100, idx=10)         # impulse at index 10

# Square wave
x = signal.square(2 * np.pi * 5 * t)          # duty cycle = 0.5
x = signal.square(2 * np.pi * 5 * t, duty=0.3)

# Sawtooth
x = signal.sawtooth(2 * np.pi * 5 * t)       # rising
x = signal.sawtooth(2 * np.pi * 5 * t, width=0)  # falling

from scipy.signal import hilbert

x = np.sin(2*np.pi*10*t)

analytic = hilbert(x)            # complex analytic signal
envelope = np.abs(analytic)      # instantaneous amplitude
inst_phase = np.unwrap(np.angle(analytic))  # instantaneous phase
inst_freq = np.diff(inst_phase) / (2*np.pi*dt)  # instantaneous frequency