#!/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). 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 def phase_vocoder_stretch(signal: np.ndarray, stretch: float, frame_size: int = FRAME_SIZE, analysis_hop: int = ANALYSIS_HOP) -> 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. """ 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) 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) # --- 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 # 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 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()