#!/usr/bin/env python3 """ Phase vocoder: time-stretch an audio signal without changing its pitch. Usage: python3 phase_vocoder.py input.wav 1.5 output.wav The third argument is the stretch factor: 2.0 makes the audio twice as long (half speed), 0.5 makes it half as long (double speed), 1.0 leaves the duration unchanged. --- The idea behind a phase vocoder ------------------------------------ Simply resampling audio to change its duration also changes its pitch (think of slowing down a vinyl record). A phase vocoder decouples the two by working in the time-frequency domain: 1. Analysis: split the signal into overlapping windowed frames and take the FFT of each frame (this is the Short-Time Fourier Transform, STFT). Each frame gives us, per frequency bin, a magnitude (loudness) and a phase (timing offset within the cycle). 2. Phase propagation: the frames are re-spaced in time -- read at the original hop size, written at a scaled hop size -- to stretch or compress the signal. Magnitudes are just carried over to the new frame spacing. Phases need more care: each bin's phase should keep advancing at that bin's true instantaneous frequency, not simply be copied. So for every frame we: a. compute how much the phase actually advanced since the last frame, compared to how much a bin at exactly its nominal frequency would have advanced (the "phase deviation"); b. unwrap that deviation into the range [-pi, pi] to find the true instantaneous frequency of the bin; c. accumulate ("propagate") the output phase using that true frequency, scaled by the new hop size. This keeps sinusoids within each bin phase-coherent across frames even though the frames are now spaced differently in time -- which is exactly what preserves pitch while changing duration. 3. Synthesis: turn each frame's (magnitude, propagated phase) back into a time-domain signal with an inverse FFT, and overlap-add the frames at the new hop size to reconstruct the stretched waveform (this is the classic overlap-add, OLA, method). Two refinements on top of that basic recipe noticeably improve audio quality: 4. Phase locking: propagating every bin independently lets bins that belong to the same partial drift out of relative phase with each other, which is heard as a smeared, reverberant "phasiness". Each frame we find the local magnitude peaks (one per partial) and re-lock every other bin's phase to keep the same offset from its nearest peak that it had in the analysis frame, instead of letting it wander off on its own propagated trajectory. 5. Transient detection: a sudden onset (drum hit, pluck, consonant) is naturally phase-coherent across all bins in that instant, but normal propagation forces its phase to follow on smoothly from the previous frame, smearing the attack. We watch for a sharp rise in spectral energy ("spectral flux") and, when one is detected, reset the output phase straight to the frame's own analysis phase instead of propagating, keeping the attack sharp. Only numpy (for FFT/array math) and scipy (for the window function and WAV file I/O) are used -- there is no specialized "phase vocoder" or "time stretch" library involved, so every DSP step above is visible in the code below. """ import sys import numpy as np from scipy.io import wavfile from scipy.signal.windows import hann # STFT parameters. A larger frame gives better frequency resolution but # worse time resolution; 2048 samples (~46 ms at 44.1 kHz) is a common # middle ground for music. The hop size (frame advance) is a quarter of # the frame size, giving 75% overlap between consecutive frames -- a # standard choice that gives the phase-propagation step enough # resolution to track instantaneous frequency accurately. FRAME_SIZE = 2048 ANALYSIS_HOP = FRAME_SIZE // 4 # A frame is flagged as a transient when its spectral flux (the total # positive rise in per-bin magnitude versus the previous frame) exceeds # this fraction of the previous frame's total magnitude. Lower catches # softer onsets but risks false positives on vibrato/tremolo; higher # only catches hard attacks. TRANSIENT_FLUX_THRESHOLD = 1.5 def _peak_region_map(magnitude: np.ndarray) -> np.ndarray | None: """Map each bin to the bin index of its nearest local magnitude peak. Used for identity phase locking: every bin is "owned" by the peak (partial) closest to it, with the boundary between two peaks' regions falling at their midpoint. Returns None if the spectrum has no local maxima (e.g. it's completely flat or silent), in which case locking is skipped for that frame. """ peaks = np.where((magnitude[1:-1] >= magnitude[:-2]) & (magnitude[1:-1] >= magnitude[2:]))[0] + 1 if len(peaks) == 0: return None boundaries = (peaks[:-1] + peaks[1:]) / 2.0 region = np.searchsorted(boundaries, np.arange(len(magnitude)), side="right") return peaks[region] def phase_vocoder_stretch(signal: np.ndarray, stretch: float, frame_size: int = FRAME_SIZE, analysis_hop: int = ANALYSIS_HOP, phase_locking: bool = True, transient_detection: bool = True) -> np.ndarray: """Time-stretch `signal` by `stretch` using the phase vocoder algorithm. `signal` is a 1-D float array. Returns a 1-D float array whose length is approximately len(signal) * stretch. `phase_locking` and `transient_detection` are on by default; turn them off to hear the effect of the plain, unlocked phase vocoder (useful for comparing quality, or as a teaching baseline). """ synthesis_hop = int(round(analysis_hop * stretch)) window = hann(frame_size, sym=False) # Pad the input so the last frame still fits, and so overlap-add at # the start/end doesn't taper off the very first/last samples. padded = np.concatenate([ np.zeros(frame_size), signal, np.zeros(frame_size), ]) num_frames = (len(padded) - frame_size) // analysis_hop # Bin frequencies, in radians/sample, that a phase-coherent sinusoid # sitting exactly at the centre of each FFT bin would advance by # every analysis hop. This is the reference we compare the *actual* # measured phase advance against, in order to find the true # instantaneous frequency of whatever is really in that bin. bin_freqs = 2 * np.pi * np.arange(frame_size // 2 + 1) / frame_size expected_phase_advance = bin_freqs * analysis_hop output_len = frame_size + num_frames * synthesis_hop output = np.zeros(output_len) # Normalizes the overlap-add so overlapping windows sum back to a # flat gain instead of amplitude ripple. window_sum = np.zeros(output_len) previous_input_phase = np.zeros(frame_size // 2 + 1) previous_magnitude = np.zeros(frame_size // 2 + 1) accumulated_output_phase = np.zeros(frame_size // 2 + 1) for i in range(num_frames): analysis_pos = i * analysis_hop frame = padded[analysis_pos: analysis_pos + frame_size] * window spectrum = np.fft.rfft(frame) magnitude = np.abs(spectrum) phase = np.angle(spectrum) # --- Transient detection --- # A transient shows up as a sudden rise in energy across many # bins at once ("spectral flux"). Propagating phase as usual # would smear that sharp attack over several frames, since it # forces this frame's phase to follow on smoothly from the # last one. If a transient is detected we skip propagation # entirely for this frame instead. spectral_flux = np.sum(np.maximum(magnitude - previous_magnitude, 0.0)) is_transient = (transient_detection and spectral_flux > TRANSIENT_FLUX_THRESHOLD * (np.sum(previous_magnitude) + 1e-8)) previous_magnitude = magnitude # --- Phase propagation --- # How far did the phase actually move since the previous frame, # versus how far a bin at its nominal frequency would have # moved? That difference is the phase deviation caused by the # real signal not sitting exactly on the bin's centre frequency. phase_deviation = phase - previous_input_phase - expected_phase_advance # Wrap the deviation into [-pi, pi] so it represents the extra # phase rotation within a single bin, not a multiple of 2*pi. phase_deviation = phase_deviation - 2 * np.pi * np.round(phase_deviation / (2 * np.pi)) # True instantaneous frequency of this bin, in radians/sample. true_freq = bin_freqs + phase_deviation / analysis_hop previous_input_phase = phase if is_transient: # The analysis phase is already coherent across bins at this # instant, so use it directly as the output phase rather # than propagating from the previous frame. accumulated_output_phase = phase.copy() else: # Advance the *output* phase by the true frequency, scaled # to the synthesis hop instead of the analysis hop. This is # what actually changes the timing without changing the # pitch: the phase keeps accumulating at the signal's real # frequency, just sampled at the new frame spacing. accumulated_output_phase = accumulated_output_phase + true_freq * synthesis_hop # --- Phase locking --- # Re-lock every bin's phase to the nearest spectral peak's # just-propagated phase, preserving the offset it had from # that peak in the analysis frame. This keeps bins that make # up the same partial phase-coherent with each other, which # is what removes the smeared "phasiness" of a plain phase # vocoder (Laroche & Dolson, 1999). if phase_locking: peak_of_bin = _peak_region_map(magnitude) if peak_of_bin is not None: accumulated_output_phase = (accumulated_output_phase[peak_of_bin] + (phase - phase[peak_of_bin])) new_spectrum = magnitude * np.exp(1j * accumulated_output_phase) new_frame = np.fft.irfft(new_spectrum, n=frame_size) * window synthesis_pos = i * synthesis_hop output[synthesis_pos: synthesis_pos + frame_size] += new_frame window_sum[synthesis_pos: synthesis_pos + frame_size] += window ** 2 # Avoid division by zero in the padded tail where windows don't overlap. nonzero = window_sum > 1e-8 output[nonzero] /= window_sum[nonzero] # Strip the padding added before the loop to line up with the # unpadded, stretched signal. return output[frame_size: frame_size + int(round(len(signal) * stretch))] def main(): if len(sys.argv) != 4: print(f"Usage: {sys.argv[0]} input.wav stretch_factor output.wav") sys.exit(1) input_path, stretch_str, output_path = sys.argv[1], sys.argv[2], sys.argv[3] stretch = float(stretch_str) if stretch <= 0: print("stretch_factor must be positive") sys.exit(1) sample_rate, samples = wavfile.read(input_path) original_dtype = samples.dtype # Normalize integer PCM to float in [-1, 1] so the DSP math (window # multiplication, FFT, phase math) isn't done on raw integer codes. if np.issubdtype(original_dtype, np.integer): max_value = np.iinfo(original_dtype).max samples_float = samples.astype(np.float64) / max_value else: samples_float = samples.astype(np.float64) max_value = 1.0 is_stereo = samples_float.ndim == 2 if is_stereo: channels = [samples_float[:, ch] for ch in range(samples_float.shape[1])] else: channels = [samples_float] stretched_channels = [phase_vocoder_stretch(ch, stretch) for ch in channels] if is_stereo: stretched = np.stack(stretched_channels, axis=1) else: stretched = stretched_channels[0] # Convert back to the original PCM dtype, clipping to avoid wraparound # on any small overshoot introduced by the overlap-add reconstruction. stretched = np.clip(stretched, -1.0, 1.0) if np.issubdtype(original_dtype, np.integer): out_samples = (stretched * max_value).astype(original_dtype) else: out_samples = stretched.astype(original_dtype) wavfile.write(output_path, sample_rate, out_samples) print(f"Wrote {output_path}: {sample_rate} Hz, " f"{len(out_samples) / sample_rate:.2f} s " f"(stretch factor {stretch})") if __name__ == "__main__": main()