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));
phout(freqind)=phadvance(sk,peaksort(tk,2))+pizero;
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.
Updates:
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.
[y,sr]=wavread(infile); % read song file
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
phold=zeros(stmo,nfft2+1); phadvance=zeros(stmo,nfft2+1);
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
phadvance(sk,peaksort(tk,2))=phadvance(sk,peaksort(tk,2))+2*pi*fdes*dtout;
% 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
phout(freqindout)=phadvance(sk,peaksort(tk,2))+pizero;
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
