# Implementing a pitch shifter

I am trying to implement a pitch shifter similar to this one, in MATLAB for quick prototyping.

My goal is to do the reverse, but I assume it is similar.

The implementation details are missing. So I am trying to piecing them together.

It looks like the phase vocoder and the SOLA methods are the preferred methods. So I did some research there.

I found this one that comes with MATLAB code and a sample chapter for explaining how it works.

I understand the majority of the code, except during the phase reconstruction, there is an odd choice of picking $$\pi$$ or 0 depends on whether or not the bin index is odd.

                pizero=pi*ones(1,length(freqind));
pcent=peaksort(tk,2)-peaksort(tk,1)+1;
indpc=(2-mod(pcent,2)):2:length(freqind);
pizero(indpc)=zeros(1,length(indpc));


The textual description in the book is equally odd with a reference pointing only to a number

If these are chosen to be ... the burst generated by the IFFT will be windowed with tapered ends


Do we know why do we make that choice? What is a tapered end? Why these choices leading to it? Why is a tapered end preferable?

I also find this one, which seems to be doing the phase vocoder, but otherwise have a fairly different implementation, in particular, I can't find any peak finding work done in that code.

That makes me wonder, what is the phase vocoder, is that the first link, or the second link, or it is a class of algorithms?

When the author said pitch shifting, in reality, it seems that he is doing pitch scaling and it looks like he didn't care if some frequency content is scaled out of bound, he just ignores them.

            /* ***************** PROCESSING ******************* */
/* this does the actual pitch shifting */
memset(gSynMagn, 0, fftFrameSize*sizeof(float));
memset(gSynFreq, 0, fftFrameSize*sizeof(float));
for (k = 0; k <= fftFrameSize2; k++) {
index = k*pitchShift;
if (index <= fftFrameSize2) {
gSynMagn[index] += gAnaMagn[k];
gSynFreq[index] = gAnaFreq[k] * pitchShift;
}
}


If I apply this to the male to female transformation. We know that the male voice frequency range is 85 to 180Hz while the female voice frequency range is 165 to 255 Hz, just doing scaling is likely to generate some 300 Hz partials which are just not there for females.

I wonder if that makes sense?

• Scale the 0 to 85Hz to 0 to 165Hz.
• Scale the 85 Hz to 180Hz to 165 to 255Hz
• Scale the 180Hz to max to 255Hz to max.

So as to minimally alter the bandwidth.

I also wanted to do the formant trick as in the first link.

Formant analysis seems to be mostly done with either LPC or cepstrum. I chose LPC.

I found this which help me to estimate the formant using LPC. I can reproduce the experiment with no problem. But I wonder

How I can implement the formant shift in the phase vocoder?

My apology for such a long and convoluted question - this is so overwhelming to me, I am just getting started with audio signal processing. If I knew all this better, I can make them clear and separate questions, at this point they all look like a mess to me.

Disclaimer, this is not homework, this is just a side project learning exercise.

Here is the code I have for now. It is based on this. It generally shifts the spectrum the way I wanted, but it is nowhere close to the pitch-shifted version as in the experiment I wanted to reproduce. I suspect the high-frequency noise I observed here is the culprit.

% This implements a phase vocoder for time-stretching/compression
% Put the file name of the input file in "infile"
%     the file to be created is "outfile"
% The amount of time-stretch is determined by the ratio of the hopin
% and hopout variables. For example, hopin=242 and hopout=161.3333
% (integers are not required) increases the tempo by
% hopin/hopout = 1.5. To slow down a comparable amount,
% choose hopin = 161.3333, hopout = 242.
% 5/2005 Bill Sethares

clf
infile='male.wav';
outfile='female.wav';
time=0;                                    % total time to process
hopin=121;                                 % hop length for input
hopout=121;                                % hop length for output

all2pi=2*pi*(0:100);                       % all multiples of 2 pi (used in PV-style freq search)
max_peak=50;                               % parameters for peak finding: number of peaks
eps_peak=0.005;                            % height of peaks
nfft=2^12; nfft2=nfft/2;                   % fft length
win=hanning(nfft)';                        % windows and windowing variables
% y are the samples, sr is the sampling rate.
siz=wavread(infile,'size');                % length of song in samples
stmo=min(siz); leny=max(siz);              % stereo (stmo=2) or mono (stmo=1)
if siz(2)<siz(1), y=y'; end
if time==0, time = 100000; end

% time is in unit seconds, and therefore time * sr equals to the number
% of samples for time seconds.
% in any case, we will not be dealing with 100,000 seconds anyway.
% it looks like a precaution to avoid very long loop only.
tt=zeros(stmo,ceil(hopout/hopin)*min(leny,time*sr)); % place for output

% Suppose hopin is an integer and we are advancing window by hopin
% every time, the first window is [1..nfft],
% the next window is [hopin + 1..hopin + nfft]
% and so on, so the last window should be
% [k * hopin + 1 .. k * hopin + nfft]
% and k * hopin + nfft should be closest to leny
% That is probably how lenseg is determined.
lenseg=floor((min(leny,time*sr)-nfft)/hopin) % number of nfft segments to process

% The fastest wave in discrete time takes two samples. Therefore, the
% period is 2 * sampling period. The associated frequency is sampling rate / 2
% which is also equal to sr * nfft2/nfft
% the remaining frequencies are scaled linearly.
ssf=sr*(0:nfft2)/nfft;                     % frequency vector

% ph stands for phase
outbeat=zeros(stmo,nfft); pold1=[]; pold2=[];
% hopin is the number of samples to advance, therefore dtin is the
% amount of time to advance
dtin=hopin/sr;                             % time advances dt per hop for input
dtout=hopout/sr;                           % time advances dt per hop for output
% for k=1239:1239
for k=1:lenseg-1                           % main loop - process each beat separately
% for k=1:3                           % main loop - process each beat separately
if k/1000==round(k/1000), disp(k), end % for display so I know where we are
% These are the indexes of the frame
indin=round(((k-1)*hopin+1):((k-1)*hopin+nfft));
for sk=1:stmo                          % do L/R channels separately
s=win.*y(sk,indin);                % get this frame and take FFT
ffts=fft(s);
mag=abs(ffts(1:nfft2+1));
ph=angle(ffts(1:nfft2+1));

% find peaks to define spectral mapping
% peaks is going to be a matrix with 3 columns, one row for
% each peak, the first column and the last column represents
% the bin for the peak, and the second column is index of the
% peak
peaks=findPeaks4(mag, max_peak, eps_peak, ssf);
% inds(1) is the index of the maximum peak, the mag is a
% vector, not a function
[dummy,inds]=sort(mag(peaks(:,2)));
% peaksort are the peaks, the columns are interpreted just like
% peaks, and the rows are sorted in ascending order
peaksort=peaks(inds,:);
% pc is the second column, that correspond to the indices of
% the peak, sorted in ascending order
pc=peaksort(:,2);

bestf=zeros(size(pc));
for tk=1:length(pc)                % estimate frequency using PV strategy

% tk is which peak
% pc(tk) is the index of the peak
% ph(pc(tk)) is the phase of the peak
%
% the phold term means the same thing for the last frame
% MATLAB allows adding a scalar to a vector, it means
% applying to all elements
%
% so dtheta is a set of phase differences
dtheta=(ph(pc(tk))-phold(sk,pc(tk)))+all2pi;

% In a period, the phasor runs 2pi radian
% In time dt, the phasor runs 2 pi dt/T = 2pi dt f radian
% And we know the phase differences are dtheta
% so fest are all the possible frequencies causing that
% peak.
fest=dtheta./(2*pi*dtin);      % see pvanalysis.m for same idea

% ssf(pc(tk)) correspond to the frequency of the coefficient.
% This implies we wanted the frequency that is closest
[er,indf]=min(abs(ssf(pc(tk))-fest));

% Now we have the best frequency found for the peak
bestf(tk)=fest(indf);          % find best freq estimate for each row
end

% generate output mag and phase

magout=mag; phout=ph;
for tk=1:length(pc)
% We have only one best frequency for a peak
fin=bestf(tk);
% Here are all the input spectral indexes
inbinfrom = peaksort(tk,1);
inbinto = peaksort(tk,3);
freqindin=(inbinfrom:inbinto);

fdes = map_frequency(fin, sr);
inbinfromfreq = inbinfrom * sr/nfft;
inbintofreq = inbinto * sr/nfft;

outbinfromfreq = map_frequency(inbinfromfreq, sr);
outbintofreq = map_frequency(inbintofreq, sr);

outbinfrom = round(outbinfromfreq * nfft / sr);
outbinto = round(outbintofreq * nfft / sr);

freqindout=(outbinfrom:outbinto);

% specify magnitude and phase of each partial

% The output magnitude is interpolated
inputx = (0:length(freqindin)-1)/(length(freqindin));
inputy = mag(freqindin);
outputx = (0:length(freqindout)-1)/(length(freqindout));
outputy = interp1(inputx, inputy, outputx,'linear','extrap');
magout(freqindout)= outputy;

% We know the real frequency, so that phase should increase by
% 2 * pi * fdes * dtout
% We are estimating only a single phadvance for the whole bin

% These weird expression zero outs some entries in pizero. This
% is a mechanism that is not totally understood around the idea
% of tapered end as in the notes
pizero=pi*ones(1,length(freqindout));
pcent=peaksort(tk,2)-peaksort(tk,1)+1;
indpc=(2-mod(pcent,2)):2:length(freqindout);
pizero(indpc)=zeros(1,length(indpc));

% Now we can set the phase for the bin
end

% reconstruct time signal (stretched or compressed)
compl=magout.*exp(sqrt(-1)*phout);
compl(nfft2+1)=ffts(nfft2+1);
compl=[compl,fliplr(conj(compl(2:(nfft2))))];
wave=real(ifft(compl));
outbeat(sk,:)=wave;
phold(sk,:)=ph;

end % end stereo
indout=round(((k-1)*hopout+1):((k-1)*hopout+nfft));
tt(:,indout)=tt(:,indout)+outbeat;
end

tt=0.8*tt/max(max(abs(tt)));
[rtt,ctt]=size(tt); if rtt==2, tt=tt'; end
wavwrite(tt,sr,16,outfile);
fclose('all');


The findPeaks4 function is given, all I did is write some comments to document my understanding.

function  peaks = findPeaks4( Amp, MAX_PEAK, EPS_PEAK, SSF )

%   This version modified from findPeaks.m by P. Moller-Nielson
%   28.3.03, pm-n. ( see http://www.daimi.au.dk/~pmn/sound/ )

SPECTRUM_SIZE=length(Amp);
% The very first entry of the spectrum is the DC value, we do not consider
% that. Therefore only [3..SPECTRUM_SIZE-1] can be a local maximum
% The first two comparison make sure it is a local maximum, and the last
% multiplication extract the value of that local maximum.
%
% The resulting vector index 1 correspond to index 3 of the Amp vector (*)
peakAmp = ( Amp(3:SPECTRUM_SIZE-1) > Amp(2:SPECTRUM_SIZE-2) ) .* ...
( Amp(3:SPECTRUM_SIZE-1) > Amp(4:SPECTRUM_SIZE) ) .* ...
Amp(3:SPECTRUM_SIZE-1);
peakPos = zeros( MAX_PEAK, 1);

maxAmp = max( peakAmp );
nPeaks = 0;
for p = 1 : MAX_PEAK
% The maximum value is m, the index of that maximum value is b
[m, b] = max( peakAmp );
if m <= ( EPS_PEAK * maxAmp )
break;
end;
% The + 2 is due to (*)
peakPos(p) = b+2;
% By setting peakAmp(b) to 0, the max in the next iteration will find the
% next max
peakAmp(b) = 0;
nPeaks = p;
end;

peakPos = sort( peakPos );
peaks = zeros( nPeaks, 3 );

% b stands for bin. WTF! Saving just 2 characters to let me guess?
last_b = 1;
for p = 1 : nPeaks
b = peakPos(MAX_PEAK-nPeaks+p);
first_b = last_b+1;
if p == nPeaks
last_b = SPECTRUM_SIZE;
else
next_b = peakPos(MAX_PEAK-nPeaks+p+1);
% Between the peaks, find the minimum
% note that rel_min is the index starting from b
% rel_min = 1 meaning Amp(b) is minimal
[dummy, rel_min] = min( Amp(b:next_b));
% Therefore last_b is the index of the minimum
last_b = b+rel_min-1;
end;
peaks(p,1) = first_b;
peaks(p,2) = b;
peaks(p,3) = last_b;
end;


The map_frequency function is what you would expect

function fdes = map_frequency(fin, sr)
fmax = sr/2;
if fin < 85
fdes = fin / 85 * 165;
else
if fin < 180
fdes = (fin - 85)/(180-85) * (255 - 165) + 165;
else
fdes = (fin - 180)/(fmax-180) * (fmax - 255) + 255;
end
end


There's a lot that you want in this question. It's a lot to expect someone to read. My understanding is that you wanted someone to explain what a phase vocoder is which is a lot to do.

Using the SOLA technique, you will need a really good pitch detector. But this technique breaks the input vocal waveform into little wavelets or grains. Then launches those grains (and overlaps them) at a rate of the intended output pitch. If the grains are not stretched or squished in time, the formants are not moved, even if the pitch is. If the grains are compressed or squished in time, the formants are moved up. If the grains are stretched in time, the formants are moved down.

If you were to do this with a phase vocoder instead of SOLA, then you might have to do some LPC analysis to get the formant filters and to flatten the spectrum before pitch shifting and then to apply the spectrum on the shifted output afterward.

These are really big topics. Not easily treated completely on an SE site.

• Thanks a lot for chiming in. Are you suggesting this is too hard for a beginner? What would you suggest is more suitable for a beginner to get started? Jan 3, 2021 at 6:40
• i don't know exactly what a "beginner" can do and cannot do. i might have some MATLAB code doing a phase vocoder to time-stretch or time-compress audio. pitch shifting would require resampling. Jan 3, 2021 at 6:58
• Yeah, I guess that is hard to tell. If you could share your code, or point me to some references that you think is good, that would be nice too. Jan 3, 2021 at 7:31
• Can you please take a look at this too? dsp.stackexchange.com/questions/72383/… Jan 3, 2021 at 7:33

I'm going to address only the phase vocoder /pitch shifter part of your question.

That makes me wonder, what is the phase vocoder, is that the first link, or the second link, or it is a class of algorithms?

The pitch shifting code you linked to is actually a simple implementation of the phase vocoder. The phase vocoder is a process/tool where you first convert a time domain signal into the frequency domain (usually via a short time fourier transform) and then use either the magnitude (or phase) information (or both) to do some modification of the original signal (whether it's pitch shifting, time stretching, frequency boosting/attenuating etc).

When the author said pitch shifting, in reality, it seems that he is doing pitch scaling and it looks like he didn't care if some frequency content is scaled out of bound, he just ignores them.

Yes, that pitch shifting algorithm discards any frequency content above the Nyquist frequency (if the pitch shifting factor > 1). You would first need to resample the signal to a higher sampling rate to keep the entire spectrum. That said, the human voice has a relatively narrow bandwidth ( a traditional telephone line has a 4kHz bandwidth if I remember correctly) so one could argue that implementing a pitch shifter even without the resampling is sufficient to get decent results (at least for smaller pitch shifting factors, say 1.0-2.0).

The problem is, however, that pitch shifting the entire spectrum isn't going to make a male voice sound like that of a female. You need to shift the formants to match that of a female instead but even that alone may (probably) still not be enough. A woman's voice can often have different qualities that are absent in men's voice (e.g breathiness) so you will probably need to emulate this as well to get good results.