Interpolation morphing between two sound files using a stochastic kernel

importPath = "~sethares/Music/SoundFilesTest/instr/";
fileNameY = "fluteA#3";
fileNameZ = "fluteG#3";
exportPath = "~sethares/Desktop/";
chan = "L";
outFile = StringJoin[exportPath,fileNameY,"2",fileNameZ,".wav"];
$HistoryLength=0;

Read in the sound files

fullFileName = StringJoin[{importPath,fileNameY,".wav"}];
q = Import[fullFileName];
dimQ = Dimensions[q[[1,1]]];
Which[(dimQ[[1]]==2)&&(chan=="L"), ySound =q[[1,1]][[1]];,
(dimQ[[1]]==2)&&(chan=="R"), ySound =q[[1,1]][[2]];,
dimQ[[1]]==2, Total[q[[1,1]]];,True,
ySound = Flatten[q[[1,1]]];];
leny = Length[ySound];
sry = q[[1,2]];
fullFileName = StringJoin[{importPath,fileNameZ,".wav"}];
q = Import[fullFileName];
dimQ = Dimensions[q[[1,1]]];
Which[(dimQ[[1]]==2)&&(chan=="L"), zSound =q[[1,1]][[1]];,
(dimQ[[1]]==2)&&(chan=="R"), zSound =q[[1,1]][[2]];,
dimQ[[1]]==2, Total[q[[1,1]]];,True,
zSound = Flatten[q[[1,1]]];];
lenz = Length[zSound];
zSound = Norm[ySound] zSound/Norm[zSound];
srz = q[[1,2]];
sr = sry;
If[sry != srz,Print["sampling rates don't match"]];

Parameters for the analysis resynthesis

samStarty=1;samEndy = 10700;samEndy=40000;
samStartz=4000;samEndz = lenz/4; samStartz = 3 lenz/4; samEndz = lenz;
nfft = 3 1024;
morphTime =5.0; 
overlap = 4;
noverlap=Round[(overlap-1) nfft/overlap];
hop=nfft-noverlap; nfft2 = nfft/2;
ssf = N[Range[0,sr/2,sr/nfft]];
ssf2 = ssf[[2]];

Auxiliary Functions

(* ind[y] makes a vector that is 41 elements long from y-20 to y+20. The function ind[y] doesn't allow indices less than one or greater than leny*)  
ind[y_,totalLeny_]:=Range[Max[1,y-20],Min[y+20,totalLeny]];

(* find index of minimum element of q *)
indMin[q_]:=Ordering[q][[1]];

(* find index of maximum element q *)
indMax[q_]:=Ordering[q][[Length[q]]];

(* n linearly spaced points from a to b *)
linSpace[a_,b_,n_]:=Range[a,b,(b-a)/(n-1)];

(* estimate frequency using best spline fit -- this one is faster because it varies the resolution of the estimation based on absolute frequency *)
estOneFreq[q_,indF_, ssf_]:=Module[{indRange, ssfIndVals, yyy,inc, tCheck, estFreq},
indRange = Range[indF-1,indF+1];
ssfIndVals = ssf[[indRange]];
yyy = Interpolation[Transpose[{ssfIndVals,q[[indRange]]}],InterpolationOrder -> 2];
inc = 2 1.95 ssfIndVals[[2]]/sr+0.05;
tCheck =  Range[ssfIndVals[[1]],ssfIndVals[[3]],inc];
estFreq  = tCheck[[indMax[yyy[tCheck]]]]];

fadeOut[sound_, len_]:=Module[{cont1, cont2, fcontour},
cont1 = ConstantArray[1,Length[sound]-len];
cont2 = Table[i,{i,1,0,-1./(len-1)}];
fcontour = sound *Join[cont1,cont2]
] ;

normalize[x_]:=Module[{d},d=Max[Abs[x]];If[d≠0,x/d,x]]

(* Find the maxPeak largest peaks in the vector q. All must be greater than the threshold defined by epsPeak*Max[q].*)
findPeak[q_, maxPeak_, epsPeak_]:= Module[{lenQ,thresh,allPeaks, intQwithzero, intQ, allPeaksEps, indAllPeaks},
lenQ = Length[q]-1;
thresh = epsPeak Max[q];
allPeaks = Boole[Table[(q[[i]]>q[[i-1]])&&(q[[i]]>q[[i+1]])&&(q[[i]]>thresh),{i,2,lenQ}]];
intQwithzero = Intersection[allPeaks Range[2,lenQ]];
intQ = Take[intQwithzero,-Length[intQwithzero]+1];
allPeaksEps = q[[intQ]];
If[Length[allPeaksEps]<maxPeak,indAllPeaks = Reverse[Ordering[allPeaksEps,-Length[allPeaksEps]]],indAllPeaks = Reverse[Ordering[allPeaksEps,-maxPeak]] ];
intQ[[indAllPeaks]] ];

specGramComplex[q_,n_, over_, m_,per_]:=Module[{qPart, window, sizQ, ran, fftQ},
If[ListQ[m],ran = Range[Round[m[[1]]],Round[m[[2]]]];, ran = Range[1,Max[Dimensions[q]]]];qPart = Partition[q[[ran]],n,Round[n/over]];
window = Table[1-(0.54+0.46 Cos[2 Pi i/n]),{i,1,n}];
sizQ = Dimensions[qPart];
fftQ = Table[Fourier[window qPart[[i]]],{i,1,sizQ[[1]]}]; 
Transpose[Take[fftQ,{1,sizQ[[1]]},{1,1+Round[per sizQ[[2]]]}] ] 
];

Starting and ending sounds are placed on the left and right hand boundaries

uy = specGramComplex[ySound,nfft,overlap,{samStarty,samEndy},0.5];
yAngleStart = Arg[uy[[Range[nfft2+1],1]]];
uy = Abs[uy];
fL = uy[[Range[Max[Dimensions[uy]]],Min[Dimensions[uy]]]];
uz = Abs[specGramComplex[zSound,nfft,overlap,{samStartz,samEndz},0.5]];
fR = uz[[Range[Max[Dimensions[uz]]],1]];

kernel=5;
nx=Round[morphTime overlap sr/nfft]-1;
ny = nfft2+1; 
indx =linSpace[0.01,Pi-0.01,nx];
Which[kernel==1,p[x_,y_,t_]:=2 Sin[x]^2/((4+  Sin[x]^(1/3)y^(1/3)) (Cosh[12. Abs[t-y]/y]-Cos[x] ));,kernel==2, p[x_,y_,t_]:=(Pi-x) Sin[x]^2/(2 Pi (Cosh[(t-y)]-Cos[x]));,
kernel==3, p[x_,y_,t_]:=Sin[x]^2/(Pi (Cosh[12. Abs[t-y]/y]-Cos[x]));,
kernel==4, p[x_,y_,t_]:=(Pi-x) Sin[x]^2/(2 Pi (Cosh[(t-y)/4]-Cos[x]));,
kernel==5,p[x_,y_,t_]:=Exp[ -200. (t-y)^2/(x^2 y^1.5)]  Sin[x] (Pi-x)/(x  Log[y+1]); 
indx =linSpace[0.3,Pi-0.3,nx];
]
uMorph = Table[indy = ind[y,ny]; p[x,y,indy].fL[[indy]]+p[Pi-x,y,indy].fR[[indy]],{y,1,ny},{x,indx}];

Resynthesize and recreate the starting sound, the morph, and the ending sound: read the morph matrix like a spectrogram.

uAll = Join[uy, uMorph[[All,2;;nx-1]], uz,2];
maxPeak=50;
epsPeak=0.0001;
dt=hop/sr;
phAdvance = yAngleStart;
ph=ConstantArray[0,nfft2+1];
numFrames=Min[Dimensions[uAll]];
tt=ConstantArray[0,Ceiling[numFrames hop + nfft -hop]];(* place for sound output *)
oldBestF={0};

Do[
indIn=Round[Range[(k-1) hop+1,(k-1)hop+nfft]];
mag=uAll[[All,k]];
pc=findPeak[mag[[1;;Round[0.95*nfft2]]],maxPeak,epsPeak];
bestF=ConstantArray[0,Length[pc]];Do[bestF[[tj]]=estOneFreq[mag,pc[[tj]],ssf],{tj,1,Length[pc]}];
magOut = ConstantArray[0,nfft2+1];
phOut=ph;
If[bestF≠{},
Do[
fdes=Max[bestF[[tk]],ssf2];
indF = Ceiling[(fdes+ssf2/2)/ssf2];
flo=ssf[[indF]]-ssf2/2;
fhi=ssf[[indF]]+ssf2/2;
lam=(fhi-fdes)/(fhi-flo);
indT = indMin[Abs[fdes-oldBestF]];
indS =  Ceiling[(oldBestF[[indT]]+ssf2/2)/ssf2];
freqInd=Range[indF-1,indF+1];
magOut[[freqInd]]=magOut[[freqInd]]+mag[[indF]] {lam ,1, 1-lam};
If[Abs[indF-indS]==1,
phAdvance[[indF]]=phAdvance[[indS]]+Pi+2 Pi fdes dt;,phAdvance[[indF]]=phAdvance[[indF]]+2 Pi fdes dt;];
phOut[[freqInd]]=phAdvance[[indF]]+{Pi, 0, Pi};,
{tk,1,Length[pc]}];
oldBestF=bestF;
compL=magOut Exp[Sqrt[-1] phOut];
compL[[1]]=0.;
compL[[nfft2+1]]=mag[[nfft2+1]];
compL2=Flatten[{compL, Reverse[Conjugate[compL[[Range[2,nfft2]]]]]}];outBeat=Re[InverseFourier[compL2,FourierParameters->{1,-1}]];
tt[[indIn]]=tt[[indIn]]+outBeat;],
{k,1,numFrames}];
leftSound =tt;

Repeat for the right channel

chan = "R";

Read in the sound files

fullFileName = StringJoin[{importPath,fileNameY,".wav"}];
q = Import[fullFileName];
dimQ = Dimensions[q[[1,1]]];
Which[(dimQ[[1]]==2)&&(chan=="L"), ySound =q[[1,1]][[1]];,
(dimQ[[1]]==2)&&(chan=="R"), ySound =q[[1,1]][[2]];,
dimQ[[1]]==2, Total[q[[1,1]]];,True,
ySound = Flatten[q[[1,1]]];];
leny = Length[ySound];
sry = q[[1,2]];
fullFileName = StringJoin[{importPath,fileNameZ,".wav"}];
q = Import[fullFileName];
dimQ = Dimensions[q[[1,1]]];
Which[(dimQ[[1]]==2)&&(chan=="L"), zSound =q[[1,1]][[1]];,
(dimQ[[1]]==2)&&(chan=="R"), zSound =q[[1,1]][[2]];,
dimQ[[1]]==2, Total[q[[1,1]]];,True,
zSound = Flatten[q[[1,1]]];];
lenz = Length[zSound];
zSound = Norm[ySound] zSound/Norm[zSound];
srz = q[[1,2]];
If[sry == srz,sr = sry];

Starting and ending sounds are placed on the left and right hand boundaries

uy = specGramComplex[ySound,nfft,overlap,{samStarty,samEndy},0.5];
yAngleStart = Arg[uy[[Range[nfft2+1],1]]];
uy = Abs[uy];
fL = uy[[Range[Max[Dimensions[uy]]],Min[Dimensions[uy]]]];
uz = Abs[specGramComplex[zSound,nfft,overlap,{samStartz,samEndz},0.5]];
fR = uz[[Range[Max[Dimensions[uz]]],1]];uMorph = Table[indy = ind[y,ny]; p[x,y,indy].fL[[indy]]+p[Pi-x,y,indy].fR[[indy]],{y,1,ny},{x,indx}];

Resynthesize and recreate the starting sound, the morph, and the ending sound: read the morph matrix like a spectrogram.

uAll = Join[uy, uMorph[[All,2;;nx-1]], uz,2];
maxPeak=50;
epsPeak=0.0001;
dt=hop/sr;
phAdvance = yAngleStart;
ph=ConstantArray[0,nfft2+1];
numFrames=Min[Dimensions[uAll]];
tt=ConstantArray[0,Ceiling[numFrames hop + nfft -hop]];(* place for sound output *)
oldBestF={0};

Do[
indIn=Round[Range[(k-1) hop+1,(k-1)hop+nfft]];
mag=uAll[[All,k]];
pc=findPeak[mag[[1;;Round[0.95*nfft2]]],maxPeak,epsPeak];
bestF=ConstantArray[0,Length[pc]];Do[bestF[[tj]]=estOneFreq[mag,pc[[tj]],ssf],{tj,1,Length[pc]}];
magOut = ConstantArray[0,nfft2+1];
phOut=ph;
If[bestF≠{},
Do[
fdes=Max[bestF[[tk]],ssf2];
indF = Ceiling[(fdes+ssf2/2)/ssf2];
flo=ssf[[indF]]-ssf2/2;
fhi=ssf[[indF]]+ssf2/2;
lam=(fhi-fdes)/(fhi-flo);
indT = indMin[Abs[fdes-oldBestF]];
indS =  Ceiling[(oldBestF[[indT]]+ssf2/2)/ssf2];
freqInd=Range[indF-1,indF+1];
magOut[[freqInd]]=magOut[[freqInd]]+mag[[indF]] {lam ,1, 1-lam};
If[Abs[indF-indS]==1,
phAdvance[[indF]]=phAdvance[[indS]]+Pi+2 Pi fdes dt;,phAdvance[[indF]]=phAdvance[[indF]]+2 Pi fdes dt;];
phOut[[freqInd]]=phAdvance[[indF]]+{Pi, 0, Pi};,
{tk,1,Length[pc]}];
oldBestF=bestF;
compL=magOut Exp[Sqrt[-1] phOut];
compL[[1]]=0.;
compL[[nfft2+1]]=mag[[nfft2+1]];
compL2=Flatten[{compL, Reverse[Conjugate[compL[[Range[2,nfft2]]]]]}];outBeat=Re[InverseFourier[compL2,FourierParameters->{1,-1}]];
tt[[indIn]]=tt[[indIn]]+outBeat;],
{k,1,numFrames}];
rightSound =tt;

Export morphed sound as file

ss = {normalize[leftSound],normalize[rightSound]};
w = Sound[SampledSoundList[ss,{sr,16}]];
Export[outFile, w];

ff=ListInterpolation[uAll]; {udimy,udimx}=Dimensions[uAll];
Plot3D[ff[i,j],{i,1,200},{j,1,udimx},PlotRange->All,PlotPoints->100]