| """resample-poly: polyphase rational resampling for audio. |
| |
| Each filter phase stores only its own support rather than being padded to a |
| common width, so an output sample multiplies through its taps and nothing else. |
| The filter matches torchaudio's construction, so results agree to floating-point |
| rounding and this is a drop-in replacement. |
| |
| High-level API: |
| Resampler(orig_freq, new_freq, ...) plan holding the polyphase bank |
| resample(x, orig_freq, new_freq) one-shot convenience |
| """ |
| import math |
| from typing import Optional |
|
|
| import torch |
|
|
| from ._ops import ops |
|
|
|
|
| def sinc_kernel(orig_freq: int, new_freq: int, lowpass_filter_width: int = 6, |
| rolloff: float = 0.99, method: str = "sinc_interp_hann", |
| beta: Optional[float] = None): |
| """Windowed-sinc polyphase bank, matching torchaudio's construction. |
| |
| Returns (kernels [L, W], width), with frequencies already reduced by their |
| greatest common divisor. |
| """ |
| if orig_freq <= 0 or new_freq <= 0: |
| raise ValueError("frequencies must be positive") |
| if lowpass_filter_width <= 0: |
| raise ValueError("lowpass_filter_width must be positive") |
| g = math.gcd(int(orig_freq), int(new_freq)) |
| M = int(orig_freq) // g |
| L = int(new_freq) // g |
|
|
| base_freq = min(M, L) * rolloff |
| width = math.ceil(lowpass_filter_width * M / base_freq) |
|
|
| idx = torch.arange(-width, width + M, dtype=torch.float64) / M |
| t = (torch.arange(0, -L, -1, dtype=torch.float64)[:, None] / L) + idx[None, :] |
| t = t * base_freq |
| t = t.clamp_(-lowpass_filter_width, lowpass_filter_width) |
|
|
| if method == "sinc_interp_hann": |
| window = torch.cos(t * math.pi / lowpass_filter_width / 2) ** 2 |
| elif method == "sinc_interp_kaiser": |
| b = 14.769656459379492 if beta is None else float(beta) |
| window = (torch.i0(b * torch.sqrt( |
| 1 - (t / lowpass_filter_width) ** 2).clamp(min=0)) / torch.i0( |
| torch.tensor(b, dtype=torch.float64))) |
| else: |
| raise ValueError(f"unknown method {method!r}") |
|
|
| tp = t * math.pi |
| k = torch.where(tp == 0, torch.ones_like(tp), torch.sin(tp) / tp) |
| k = k * window * (base_freq / M) |
| return k.to(torch.float32), width, L, M |
|
|
|
|
| class Resampler: |
| """Cached polyphase plan for a fixed rate pair.""" |
|
|
| def __init__(self, orig_freq: int, new_freq: int, |
| lowpass_filter_width: int = 6, rolloff: float = 0.99, |
| method: str = "sinc_interp_hann", beta: Optional[float] = None): |
| kern, width, L, M = sinc_kernel(orig_freq, new_freq, |
| lowpass_filter_width, rolloff, |
| method, beta) |
| self.orig_freq, self.new_freq = int(orig_freq), int(new_freq) |
| self.L, self.M, self.width = L, M, width |
| self.kernel = kern |
|
|
| |
| |
| |
| starts, lens, offs, taps = [], [], [], [] |
| cursor = 0 |
| for i in range(L): |
| row = kern[i] |
| nz = torch.nonzero(row, as_tuple=False).flatten() |
| if nz.numel() == 0: |
| starts.append(-width); lens.append(0); offs.append(cursor) |
| continue |
| j0, j1 = int(nz[0]), int(nz[-1]) + 1 |
| starts.append(j0 - width) |
| lens.append(j1 - j0) |
| offs.append(cursor) |
| taps.append(row[j0:j1].contiguous()) |
| cursor += j1 - j0 |
| self.start = torch.tensor(starts, dtype=torch.int32) |
| self.len = torch.tensor(lens, dtype=torch.int32) |
| self.off = torch.tensor(offs, dtype=torch.int64) |
| self.taps = (torch.cat(taps) if taps else |
| torch.zeros(0, dtype=torch.float32)).contiguous() |
|
|
| |
| |
| |
| |
| |
| |
| |
| self.taps_per_output = float(self.taps.numel()) / float(L) |
| self.fused = self.taps_per_output >= 16.0 |
|
|
| @property |
| def density(self) -> float: |
| """Stored taps over the padded-kernel size the conv1d form would use.""" |
| return float(self.taps.numel()) / float(self.kernel.numel()) |
|
|
| def _bank(self, x: torch.Tensor) -> torch.Tensor: |
| """Padded-bank formulation, for rate pairs whose phases are too short |
| for the fused path to pay for its reduction.""" |
| B, T = x.shape |
| pad = torch.nn.functional.pad(x, (self.width, self.width + self.M)) |
| r = torch.nn.functional.conv1d(pad[:, None], self.kernel[:, None, :], |
| stride=self.M) |
| return r.transpose(1, 2).reshape(B, -1)[:, :self.out_len(T)] |
|
|
| def out_len(self, n: int) -> int: |
| return int(math.ceil(self.L * n / self.M)) |
|
|
| def __call__(self, x: torch.Tensor) -> torch.Tensor: |
| """x [T] or [B, T] -> resampled, matching torchaudio's output length.""" |
| squeeze = x.dim() == 1 |
| if squeeze: |
| x = x.unsqueeze(0) |
| if x.dim() != 2: |
| raise ValueError("x must be [T] or [B, T]") |
| x = x.to(torch.float32).contiguous() |
| B, T = x.shape |
| if not self.fused: |
| y = self._bank(x) |
| return y.squeeze(0) if squeeze else y |
| out = torch.empty(B, self.out_len(T), dtype=torch.float32) |
| ops.rp_resample(out, x, self.taps, self.start, self.len, self.off, |
| self.L, self.M) |
| return out.squeeze(0) if squeeze else out |
|
|
|
|
| def resample(x: torch.Tensor, orig_freq: int, new_freq: int, |
| **kw) -> torch.Tensor: |
| """One-shot resample. Prefer Resampler when converting repeatedly.""" |
| return Resampler(orig_freq, new_freq, **kw)(x) |
|
|
|
|
| __all__ = ["Resampler", "resample", "sinc_kernel"] |
|
|
