# 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. – Arpit Jain Apr 10 '17 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. – Phil Freihofner Apr 11 '17 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? – robert bristow-johnson Apr 11 '17 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. – Phil Freihofner Apr 11 '17 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 – ederwander Apr 11 '17 at 18:37

The phase vocoder was first implemented by Flanagan and Golden using analogic 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 synthesys.

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)

 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. – Phil Freihofner Apr 21 '17 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:

http://www.ee.columbia.edu/~dpwe/LabROSA/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 dpwe@ee.columbia.edu.  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. – Phil Freihofner Apr 11 '17 at 3:16
• there are advantages and maybe a disadvantage of using the phase vocoder for time stretching. – robert bristow-johnson Apr 18 '17 at 4:37