# Audio time stretching, without pitch shifting

This might be a Sound Design question, or a StackOverflow question since I am attempting to do this with Java.

I would like to play back a sound at the same pitch, but stretched out in time.

My attempt has been the following strategy:

I break the sound into many granules of 2K sample points at 44100 fps (have also tried other fixed lengths). For example, if the first granule runs from 0 to 2000, the second might start at sample 250 and run to 2250. The third from 500 to 2500, etc.

The question is how to string these together in a way that will create a minimum of artifacts. I'm new to granular synthesis. It seems that a lot of writing about working with granules discusses their smooth entry and exit as a function of the envelope. But with what I want to do, I've been thinking in terms of cross-fading, and thus, coming up with an envelope where the attack is a mirror image of the decay. Is this faulty?

For cross-fades, I've tried the following functions, described here in terms of a LUT that is sized to match the number of overlap samples:

• linear:

$$\text{lut}[i] = \frac{(n-i)}{n}$$

• a cosine function:

$$\text{lut}[i] = \cos \left( \frac{\pi}{2} \times \frac{i}{n} \right)$$

• an equation that eases in and out (provided by a friend):

$$\text{lut}[i] = 3\left( \frac{i}{n} \right)^2 - 2\left( \frac{i}{n} \right)^3$$

• squares:

$$\text{lut}[i] = \frac{i \times i}{n \times n}$$

• square root:

$$\text{lut}[i] = \sqrt{ \frac{(n - i)}{n} }$$

Once the LUT is made, the decay factor for the first granule is gotten by iterating from front to back, and the attack of the second granule via iterating from back to front. I've tried different sizes of LUT for the different functions.

All of these are flawed, in terms of the resulting sound, even for longer overlaps, where the attack and decay are almost the entire granule length.

Of these, when dealing with "macro" granules in another context (lengths of 1/2 second with cross-fades of 1/5 of a second), I found that the "cos" function above did the best job of keeping the volume. But for small granules of about 0.05 of a second, for this time-stretching task, it seems the considerations must be rather different.

I know I've heard pretty nice time stretching effects before. Is there an entirely different approach that has to be used? Can an improvement be accomplished by using a different envelope shape?

Progress report, 4/11/17 In order to better understand the range of solutions offered, I am reading the wiki article Audio time-scale/pitch modification, and an article cited therein A Review of Time-Scale Modification of Musical Signals.

Progress, 4/16/17

I now have a valid implementation of an OLA algorithm using a Hamming window working. This was tricky to code and debug! A big help in understanding came by referring to concepts in the Driedger/Muller article cited above, in particular, the concepts of "analysis granule" and "synthesis granule" and their respective "hops" as a way to keep things straight.

A test that helped, in terms of debugging, was to modify my original naive attempt into something that should theoretically produce the equivalent of the OLA method: with my original approach (based on cross-fades), I used the Hamming window algorithm to control the cross-fading, and made the entire signal either fading in or out, and thus equal to a single granule in the OLA method.

Eventually, I got both methods to produce the same outcome. However, there is a strong pitch component that occurs with stretching, and correlates with the amount of stretch (the longer the stretch, the lower the tone). I will link an audio file of examples as soon as I can get to it. It may be an indication of something awry in my coding. [The artifacts are actually not so bad for a music-based example, the pitchiness seems to arise when there is some noise in the original. I haven't gotten to the bottom of this yet. But the OLA is definitely an improvement over my first attempts.]

I don't think I'll have the time to go deeper into the more advanced suggestions until after the concert, and so, I apologize for not designating an accepted answer. I want to actually try to implement the answer before selecting it.

Meanwhile, will probably be using Audacity's Paulstretch or its variable stretch tool to make multiple, fixed-length takes, rather than attempt to do the time-stretching procedurally during playback.

• you can try hanning and hamming envelopes with 50% overlap to reduce artifacts. usually for speech signals I use them. Apr 10, 2017 at 6:50
• Thanks! I'm going to have to do some research to figure out how to set them up, but that sounds like a promising approach. Apr 11, 2017 at 3:00
• if you're time-stretching, staying only in the time domain (no FFT, no phase vocoder), then you might want to consider using an autocorrelation function of some kind to get give you a better idea of what snippets of audio to be crossfading. is the audio very complex (like an orchestra or a mix of many sounds)? or, at the other extreme, is the audio normally a single note at a time? Apr 11, 2017 at 4:51
• My main interest is in working on snippets extracted from a recording of ravens in the wild. I'm hoping to play predetermined snippets with procedural, chance-derived amounts of time-stretching. It might be nice to also be able to jump into random spots on a sound cue and do this. My background is composing/performing and this is for a piece I am hoping to realize for an upcoming benefit performance. Apr 11, 2017 at 10:42
• In the past I have used different approaches, I think that I build all known types of time stretch, SOLA, PSOLA, TDHS, Phase Vocoder, WSOLA, etc I like so much of the Keith Lent, the Lent code is amazing for Monophonic signals, this works like PSOLA, If you read about PSOLA you will see that some guys try find glottal closure Instant (GCI) as epochs, I never try get CGI's, I always get the period/pitch from every frame, @robertbristow-johnson wrote a very nice paper about keith Lent algorithm Apr 11, 2017 at 18:37

The phase vocoder was first implemented by Flanagan and Golden[1] using analog filters bank, later Portnoff apply the same concept digitally using FFT.

The phase vocoder use the difference of successive phase spectra ∆φ[i], it will help find the Instantaneous Frequency used in phase synthesis.

A little math:

$$q[i] = \frac{N}{2πH} princarg \left[φ_l[i] − φ_{l−1}[i] − \frac {2πH}{N}i \right]$$

q[i] is the deviation of the partial’s Instantaneous Frequency from the bin:

Princarg maps the phase in a range of ±pi.

The Instantaneous Frequency can be determined as:

$$InstFreqs[i] = (i + q[i]) \left(\frac{Fs}{N}\right)$$

We can use matlab to show how it is done, here few lines of codes, not so hard to do, remember the standard implementation give you some artifacts like phasiness:

factor=2;

winsize=4096;
fftsize=winsize;
window=hann(winsize);
shop=winsize/4;
ahop=floor(shop/factor);

SignalLen=length(signal);
num_win = floor((SignalLen-winsize)/ahop);

OutLen=(num_win-1)*shop+winsize;
Out=zeros(OutLen,1);

TWOPI=2*pi;
CENTERFREQ = [0:fftsize-1]*TWOPI*ahop/fftsize; %Piece of equation showed in q
first=1;
PosIn=1;
PosOut=1;

for win_count=1:num_win
if first==1
framed=signal(1:winsize).*window;
X=fft(fftshift(framed),fftsize);
Mag=abs(X);
Pha=angle(X);
PhaSy=Pha;
first=0;
else
framed = signal(PosIn:PosIn+winsize-1) .* window; %framed and windowed / current position analysis
X=fft(fftshift(framed),fftsize);  %Apply FFT whith circular shift
Mag=abs(X); % Get the Magnitude
Pha=angle(X); %Get the Phase
phaseDiff = Pha - old_pha; %Difference between the current and previous phase
phaseDiff = phaseDiff - CENTERFREQ'; %Expected phase (unwrapped phase)
dphi = phaseDiff - TWOPI * round(phaseDiff /TWOPI); %principal argument, MAP phase to +/- pi
freq = (CENTERFREQ + dphi') /ahop; %true frequency
PhaSy = old_PhaSy + shop*freq'; %Phase synthesis

end

Y = Mag .* ( cos(PhaSy) + sqrt(-1) *(sin(PhaSy)) ); %Resynthesis
y_out = fftshift(real(ifft(Y,fftsize))).*window; %back to time domain

Out(PosOut:PosOut+winsize-1)= Out(PosOut:PosOut+winsize-1) + y_out; %Overlap and add

old_pha=Pha;
old_PhaSy=PhaSy;
PosIn = PosIn + ahop;
PosOut=PosOut+shop;
end

%Play
sound(Out,Fs);


The code show all the magic, it is an standard implementation of the phase vocoder, can be easy port to C/C++ or java, For better results you will need lock the phase (read Miller Puckette paper)!

Take a little time to read Portnoff

PS: Sorry my boss does not leave share TSM in time domain, but the crucial step is the Period Track (pitch track), one basic TDHS or Keith Lent code can be write in matlab with approximately 20-30 lines (Of course I'm not counting the pitch track lines)

[1] J. L. Flanagan and R. M. Golden, “Phase vocoder,” Bell Systems Technical Journal, vol. 45, pp. 1493–1509, 1966.

• Thank you for the contribution! I suspect that this was translated from another language, and could benefit from a little editing. I do not feel qualified to do so. My next step is going to be to try implementing some sort of transient detecting scheme with OLA. If the various attempts at providing correlation don't work out, at that point I'll look more deeply into writing some sort of phase vocoder and will be making further use of your answer. Apr 21, 2017 at 15:23

Use the "phase vocoder" algorithm which does this nicely with FFT's. Here are some links with further information and Matlab implementations. You need to be careful in how each FFT block is added together as the phase has a big impact, thus using the implementations already worked out below will readily get you to the results you are trying to achieve and also help you understand the fundamental details involved:

https://en.wikipedia.org/wiki/Phase_vocoder

This link contains Matlab source code for a phase vocoder implementation by Dr. Dan Ellis at Columbia University:

https://www.ee.columbia.edu/~dpwe/resources/matlab/pvoc/

Here is the top level code along with the headers describing the additional routines called. All code is at the link above:

function y = pvoc(x, r, n)
% y = pvoc(x, r, n)  Time-scale a signal to r times faster with phase vocoder
%      x is an input sound. n is the FFT size, defaults to 1024.
%      Calculate the 25%-overlapped STFT, squeeze it by a factor of r,  inverse spegram.
% 2000-12-05, 2002-02-13 [email protected].  Uses pvsample, stft, istft
% $$Header: /home/empire6/dpwe/public_html/resources/matlab/pvoc/RCS/pvoc.m,v 1.3 2011/02/08 21:08:39 dpwe Exp$$

if nargin < 3
n = 1024;
end

% With hann windowing on both input and output,
% we need 25% window overlap for smooth reconstruction
hop = n/4;
% Effect of hanns at both ends is a cumulated cos^2 window (for
% r = 1 anyway); need to scale magnitudes by 2/3 for
% identity input/output
%scf = 2/3;
% 2011-02-07: this factor is now included in istft.m
scf = 1.0;

% Calculate the basic STFT, magnitude scaled
X = scf * stft(x', n, n, hop);

% Calculate the new timebase samples
[rows, cols] = size(X);
t = 0:r:(cols-2);
% Have to stay two cols off end because (a) counting from zero, and
% (b) need col n AND col n+1 to interpolate

% Generate the new spectrogram
X2 = pvsample(X, t, hop);

% Invert to a waveform
y = istft(X2, n, n, hop)';

function D = stft(x, f, w, h, sr)
% D = stft(X, F, W, H, SR)                       Short-time Fourier transform.
%   Returns some frames of short-term Fourier transform of x.  Each column of the result
%   is one F-point fft (default 256); each successive frame is offset by H points (W/2)
%   until X is exhausted.
%   Data is hann-windowed at W pts (F),
%       or rectangular if W=0, or with W if it is a vector.
%   Without output arguments, will plot like sgram (SR will get axes right,
%       defaults to 8000).

function x = istft(d, ftsize, w, h)
% X = istft(D, F, W, H)                   Inverse short-time Fourier transform.
%   Performs overlap-add resynthesis from the short-time Fourier transform
%   data in D.  Each column of D is taken as the result of an F-point
%   fft; each successive frame was offset by H points (default
%   W/2, or F/2 if W==0). Data is hann-windowed at W pts, or
%       W = 0 gives a rectangular window (default);
%       W as a vector uses that as window.
%       This version scales the output so the loop gain is 1.0 for
%       either hann-win an-syn with 25% overlap, or hann-win on
%       analysis and rect-win (W=0) on synthesis with 50% overlap.

function c = pvsample(b, t, hop)
% c = pvsample(b, t, hop)   Interpolate an STFT array according to the 'phase vocoder'
%     b is an STFT array, of the form generated by 'specgram'.
%     t is a vector of (real) time-samples, which specifies a path through
%     the time-base defined by the columns of b.  For each value of t,
%     the spectral magnitudes in the columns of b are interpolated, and
%     the phase difference between the successive columns of b is
%     calculated; a new column is created in the output array c that
%     preserves this per-step phase advance in each bin.
%     hop is the STFT hop size, defaults to N/2, where N is the FFT size
%     and b has N/2+1 rows.  hop is needed to calculate the 'null' phase
%     advance expected in each bin.
%     Note: t is defined relative to a zero origin, so 0.1 is 90% of
%     the first column of b, plus 10% of the second.


This tutorial also looks good:

http://sethares.engr.wisc.edu/vocoders/phasevocoder.html

• I hadn't realized that this is what phase vocoder's are about. Thanks! Ultimately this seems like the way to go. Not sure where I'm going to get a Java FFT, yet. Maybe I'll try rewriting the BASIC program in Engineer and Scientist's Guide to DSP. I've been reading this book as my main source for learning more about the subject. Apr 11, 2017 at 3:16
• there are advantages and maybe a disadvantage of using the phase vocoder for time stretching. Apr 18, 2017 at 4:37

The program Rubberband is hard to beat for: quality, amount of stretchiness, number of options, and source code; along with the ability to run on a mobile phone. GitHub.

Command line options:

Rubber Band
An audio time-stretching and pitch-shifting library and utility program.

Usage: ./build/rubberband [options] <infile.wav> <outfile.wav>

You must specify at least one of the following time and pitch ratio options:

-t<X>, --time <X>       Stretch to X times original duration, or
-T<X>, --tempo <X>      Change tempo by multiple X (same as --time 1/X), or
-T<X>, --tempo <X>:<Y>  Change tempo from X to Y (same as --time X/Y), or
-D<X>, --duration <X>   Stretch or squash to make output file X seconds long

-p<X>, --pitch <X>      Raise pitch by X semitones, or
-f<X>, --frequency <X>  Change frequency by multiple X

The following options provide ways of making the time and frequency ratios
change during the audio:

-M<F>, --timemap <F>    Use file F as the source for time map

A time map (or key-frame map) file contains a series of lines, each with two
sample frame numbers separated by a single space. These are source and
target frames for fixed time points within the audio data, defining a varying
stretch factor through the audio. When supplying a time map you must specify
an overall stretch factor using -t, -T, or -D as well, to determine the
total output duration.

--pitchmap <F>   Use file F as the source for pitch map

A pitch map file contains a series of lines, each with two values: the input
sample frame number and a pitch offset in semitones, separated by a single
space. These specify a varying pitch factor through the audio. The offsets
are all relative to an initial offset specified by the pitch or frequency
option, or relative to no shift if neither was specified. Offsets are
not cumulative. This option implies realtime mode (-R) and also enables a
high-consistency pitch shifting mode, appropriate for dynamic pitch changes.
Because of the use of realtime mode, the overall duration will not be exact.

--freqmap <F>    Use file F as the source for frequency map

A frequency map file is like a pitch map, except that its second column
lists frequency multipliers rather than pitch offsets (like the difference
between pitch and frequency options above).

The following options affect the sound manipulation and quality:

-2,    --fast           Use the R2 (faster) engine

This is the default (for backward compatibility) when this tool is invoked
as "rubberband". It was the only engine available in versions prior to v3.0.

-3,    --fine           Use the R3 (finer) engine

This is the default when this tool is invoked as "rubberband-r3". It almost
always produces better results than the R2 engine, but with significantly

-F,    --formant        Enable formant preservation when pitch shifting

This option attempts to keep the formant envelope unchanged when changing
the pitch, retaining the original timbre of vocals and instruments in a
recognisable way.

-c<N>, --crisp <N>      Crispness (N = 0,1,2,3,4,5,6); default 5

This option only has an effect when using the R2 (faster) engine. See below
for details of the different levels.

The remaining options fine-tune the processing mode and stretch algorithm.
The default is to use none of these options.
The options marked (2) currently only have an effect when using the R2 engine
(see -2, -3 options above).

-R,    --realtime       Select realtime mode (implies --no-threads).
This utility does not do realtime stream processing;
the option merely selects realtime mode for the
stretcher it uses
(2)      --threads        Assume multi-CPU even if only one CPU is identified
(2)      --no-transients  Disable phase resynchronisation at transients
(2)      --bl-transients  Band-limit phase resync to extreme frequencies
(2)      --no-lamination  Disable phase lamination
(2)      --window-long    Use longer processing window (actual size may vary)
(2)      --window-short   Use shorter processing window
(2)      --smoothing      Apply window presum and time-domain smoothing
(2)      --detector-perc  Use percussive transient detector (as in pre-1.5)
(2)      --detector-soft  Use soft transient detector
--pitch-hq       In RT mode, use a slower, higher quality pitch shift
--centre-focus   Preserve focus of centre material in stereo
(at a cost in width and individual channel quality)
--ignore-clipping Ignore clipping at output; the default is to restart
with reduced gain if clipping occurs
-L,    --loose          [Accepted for compatibility but ignored; always off]
-P,    --precise        [Accepted for compatibility but ignored; always on]

-d<N>, --debug <N>      Select debug level (N = 0,1,2,3); default 0, full 3
(N.B. debug level 3 includes audible ticks in output)

The following options are for output control and administration:

-q,    --quiet          Suppress progress output
-V,    --version        Show version number and exit
-h,    --help           Show the normal help output
-H,    --full-help      Show the full help output

"Crispness" levels: (2)
-c 0   equivalent to --no-transients --no-lamination --window-long
-c 1   equivalent to --detector-soft --no-lamination --window-long (for piano)
-c 2   equivalent to --no-transients --no-lamination
-c 3   equivalent to --no-transients
-c 4   equivalent to --bl-transients
-c 5   default processing options
-c 6   equivalent to --no-lamination --window-short (may be good for drums)


The following YouTube video compares these algorithms: Elastique 3.2.3 Pro; Elastique Soloist; Elastique Efficient; Rubber Band; SoundTouch High Quality; Simple Windowed 50ms, 50% fade; CDP extension Stretch Time mode.