%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: you should copy and paste the desired %
% transformation in the sms.m file in the line%
% TRANSFORMATIONS. It is strongly recommended %
% you read the text explanation in the chapter%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===== Filtering with arbitrary resoultion =====
Filter=[ 0  2099 2100 3000 3001 22050 ;
         0  0    1    1    0    0 ];
[syniloc,ind] = sort(iloc);
FilterEnvelope = interp1(Filter(1,:)',Filter(2,:)',syniloc/N*SR);
synival = ival(ind)+(20*log10(max(FilterEnvelope,10^-9)));
synival(ind)=synival;
syniloc(ind)=syniloc;

%=== voice to clarinet ===
syniloc = iloc;
synival = ival;
if (isHarm == 1)
for i=1:nSines
harmNum = round(((iloc(i)-1)/w1Length*SR/2)/pitchvalue);
if (mod(harmNum,2)==0)  % case of an even harmonic number
synival(i) = -100;
end
end
end

%==== Frequency Shift =====
fstretch = 300;  % frequency shift in Hz
syniloc = iloc + round(fstretch/SR*N);
syniloc = syniloc.*(syniloc<=N/2);

%===== Frequency Stretch  =====
fstretch = 1.1;
[syniloc,ind] = sort(iloc);
syniloc = syniloc.*((fstretch).^[0:nSines-1]');
syniloc = syniloc.*(syniloc<=N/2);

%==== Frequency Scale =====
fscale = 1.6;  % frequency scaling factor
syniloc = iloc * fscale;
syniloc = syniloc.*(syniloc<=N/2);

%===== Pitch transposition with timbre preservation =====
if (isHarm == 1)
	pt = 2.;                 % pitch transposition factor
	[spectralShape, shapePos] = CalculateSpectralShape(iloc, ival, MinMag, N);
[syniloc, synival] = PitchTransposition(iloc, ival, spectralShape, shapePos, pt, N);
	%--- comb filtering the residual
	CombCoef = 1;
	if (isHarm==1)
		resfft = CombFilter(resfft, N, SR/(pitchvalue*pt), CombCoef);
	end
end


%===== Pitch discretization to temperate scale =====
if (pitchvalue ~= 0)
	nst = round(12*log(pitchvalue/55)/log(2));
  discpitch = 55*((2^(1/12))^nst);           % discretized pitch
  pt = discpitch/pitchvalue ;            % pitch transposition factor
	[spectralShape, shapePos] = CalculateSpectralShape(iloc, ival, MinMag, N);
  [syniloc, synival] = PitchTransposition(iloc, ival, spectralShape, shapePos, pt, N);
  %--- comb filtering the residual
  CombCoef = 1;
  if (isHarm==1)
  	resfft = CombFilter(resfft, N, SR/(pitchvalue*pt), CombCoef);
	end
end;

%===== vibrato and tremolo =====
if (isHarm == 1)
	vtf = 5;       % vibrato-tremolo frequency in Hz
	va  = 10;     % vibrato depth in percentil
	td  = 3;       % tremolo depth in dB
	synival = ival + td*sin(2*pi*vtf*pin/SR);   % tremolo
	pt = 1 + va/200*sin(2*pi*vtf*pin/SR);    % pitch transposition factor   
	[spectralShape, shapePos] = CalculateSpectralShape(iloc, ival, MinMag, N);
	[syniloc, synival] = PitchTransposition(iloc, ival, spectralShape, shapePos, pt, N);
    %--- comb filtering the residual
    CombCoef = 1;
    resfft = CombFilter(resfft, N, SR/(pitchvalue*pt), CombCoef);
end

%===== Spectral Shape Shift (positive or negative) =====
	sss  = -200;                % spectral shape shift value in Hz
	%--- spectral shape computation
	[spectralShape, shapePos] = CalculateSpectralShape(iloc, ival, MinMag, N);
	%--- spectral shape shift
	syniloc = zeros(nSines,1);
	if shapePos > 1
		[shiftedSpectralShape,shapePos] = SpectralShapeShift(sss, iloc, ival, spectralShape, shapePos, N, SR);
	end
	syniloc = iloc;
	%--- linear interpolation of the spectral shape for the synival computation
	if shapePos > 1
		synival = interp1(shiftedSpectralShape(1,1:shapePos+1)', shiftedSpectralShape(2,1:shapePos+1)', syniloc, 'linear');
	else
		synival = ival;
	end



%===== gender change: woman to man =====	
tr='m2f'; %male to female
%tr='f2m'; %female to male
if (isHarm == 1)
	pitchmin=100;
	pitchmax=500;
	sssmax = 50;
	if (pitchvalue<pitchmin)
		sss = 0;
	elseif (pitchvalue>pitchmax)
		sss = sssmax;
	else
		sss = (pitchvalue-pitchmin)/((pitchmax-pitchmin)/sssmax);
	end
	if(tr=='f2m')
		sss=-sss;
		pt=0.5;
	else
		pt=2;
	end
   

%===== harmonizer =====	
nVoices = 2;
nSynSines = nSines*(1+nVoices);
[spectralShape, shapePos] = CalculateSpectralShape(syniloc(1:nSines), synival(1:nSines), MinMag, N);
synival(1:nSines) = synival(1:nSines) - 100;
pt = [1.3 1.5]; % pitch transposition factor
ac = [-1 -2]; % amplitude change factor in dB

for i=1:nVoices
   [tmpsyniloc, tmpsynival] = ...
      PitchTransposition(syniloc(1:nSines), synival(1:nSines), spectralShape, shapePos, pt(i), N);
   tmpsynival = tmpsynival + ac(i);
   syniloc(nSines*i+1:nSines*(i+1)) = tmpsyniloc;
   synival(nSines*i+1:nSines*(i+1)) = tmpsynival;
   if (pin > 0)
      distminindex(nSines*i+1:nSines*(i+1)) = distminindex(1:nSines)+nSines*i;
   end
end


%===== hoarseness =====	
rgain = 2; % gain factor applied to the residual


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: for doing the morph you actually need %
% to modify the sms.m code, please read the   %
% text.                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%=== Morphing ===
%--- sorting the frequencies in bins; iloc1, iloc2 
[syniloc1, ind1] = sort(iloc1);
synival1   = ival1(ind1);
syniphase1 = iphase1(ind1);
distminindex1 = distminindex1(ind1);
[syniloc2, ind2] = sort(iloc2);
synival2   = ival2(ind2);
syniphase2 = iphase2(ind2);
distminindex2 = distminindex2(ind2);
%--- interpolation -----
alpha = 0.5;      % interpolation factor    
syniloc     = alpha*syniloc1 + (1-alpha)*syniloc2;
synival     = alpha*synival1 + (1-alpha)*synival2;
%--- pitch computation
isHarmsyn = isHarm1*isHarm2;
if (isHarmsyn ==1)
	npitchpeaks = min(50,nPeaks);
	[pitchvalue,pitcherror] = TWM(syniloc(1:npitchpeaks),synival(1:npitchpeaks),N,SR);
else
	pitchvalue = 0;
	pitcherror = 0;
end
if (pin==0)	%--- for the first frame
	nNewPeaks = nSines;
else
	%--- creation of new born tracks
	for i=1:nSines
		if (previoussyniloc(i)==0)
			[previoussyniloc(i), previoussynival(i)] = CreateNewTrack ...
				(syniloc, synival, previoussyniloc, previoussynival, nSines, MinMag);
			nNewPeaks = nNewPeaks - 1;
		end
	end
	%--- peak tracking of the peaks of the synthetized signal
	[syniloc, synival, syniphase, previoussyniloc, previoussynival, distminindex] = ...
	peakTrackSimple (nSines, nPeaks, N, SR, pitchvalue, syniloc, synival, ...
		syniphase, isHarmsyn, previoussyniloc, previoussynival);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: these are auxiliary functions that are%
% used in some of the transformations, read   %
% the text in order to decide the ones you    %
% in your particular case.                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===== Spectral Shape Computation =====
[spectralShape, shapePos] = CalculateSpectralShape(iloc, ival, MinMag, N);
	%--- spectral shape shift
	syniloc = zeros(nSines,1);
	if shapePos > 1
		[shiftedSpectralShape,shapePos] = SpectralShapeShift(sss, iloc, ival, spectralShape, shapePos, N, SR);
	end
	syniloc = iloc;
	%--- linear interpolation of the spectral shape for the synival computation
	if shapePos > 1
		synival = interp1(shiftedSpectralShape(1,1:shapePos+1)', shiftedSpectralShape(2,1:shapePos+1)', syniloc, 'linear');
	else
		synival = ival;
	end
	%--- pitch transposition
	pt = 0.5;                    
	[syniloc, synival] = PitchTransposition(iloc, ival, spectralShape, shapePos, pt, N);
	%--- comb filtering the residual
	CombCoef = 1;
	if (isHarm==1)
		resfft = CombFilter(resfft, N, SR/(pitchvalue*pt), CombCoef);
	end
end   

%===== Spectral Shape Shift =====
function [shiftedSpectralShape,shapePos] = SpectralShapeShift(sss, iloc, ival, spectralShape, shapePos, N, SR)
shiftedSpectralShape = spectralShape;
sssn = round (sss*N/SR);   % spectral shape shift in number of bins
if sssn > 0
   shiftedSpectralShape(1,2:shapePos) = min(N/2,spectralShape(1,2:shapePos) + sssn);
   for i=shapePos:-1:1
      if shiftedSpectralShape(1,i) < N/2
         shapePos = i;
         break;
      end;
   end;
else
   shiftedSpectralShape(1,2:shapePos) = max(1,spectralShape(1,2:shapePos) + sssn);
   for i=1:shapePos
      if shiftedSpectralShape(1,i) > 1
         shiftedSpectralShape(1,2:2+shapePos+1-i) = shiftedSpectralShape(1,i:shapePos+1);
         shapePos = shapePos-(i-2);
         break;
      end;
   end;
end;


%===== Comb Filter =====
function combFT = CombFilter(FT, N, delay, ampl)
%
%===> Comb filter in the frequential domain
% data:
%       combFT:  FT of the signal comb filtered
%       FT:      FT of the signal to filtered
%       N:       size of the FT
%       delay:   delay to apply, in samples
%       ampl:    amplitude of the multiplying coefficient (in [0,1])

coef = ampl * exp(-2*j*pi*delay*(0:N-1)/N)';
combFT = FT .* (1 + coef + coef.^2);

%===== Pitch Transposition =====
function [syniloc, synival] = PitchTransposition(iloc, ival, spectralShape, shapePos, pt, N)
syniloc = iloc.*pt;
syniloc = syniloc.*(syniloc<=N/2);
%--- linear interpolation of the spectral shape for the synival computation
if shapePos > 1
   synival = interp1(spectralShape(1,:)',spectralShape(2,:)',syniloc);
else
   synival = ival;
end