Different scalings of y work better at different frequencies

This sets up morphs at y values around 8, 16, 32 and 64. The best scaling in the Cosh term at thees values is y, y/2, y/4, and y/8. Demonstrated here with pictures.

p[x_,y_]:=Sin[x]^2/(Cosh[y]-Cos[x]);
Print[p[x,y]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-16]+DiracDelta[t-32]+DiracDelta[t-64];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-18]+DiracDelta[t-36]+DiracDelta[t-72];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,t-y] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,t-y] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,0,80},PlotRange->All,MaxRecursion->5, PlotPoints->{20,100},ImageSize->800]

above saved as harm1.png

p[x_,y_]:=Sin[x]^2/(Cosh[y/2]-Cos[x]);
Print[p[x,y]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-15]+DiracDelta[t-30]+DiracDelta[t-58];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-17]+DiracDelta[t-34]+DiracDelta[t-66];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,t-y] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,t-y] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,0,80},PlotRange->All,MaxRecursion->5, PlotPoints->{20,100}]

p[x_,y_]:=Sin[x]^2/(Cosh[y/4]-Cos[x]);
Print[p[x,y]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-15]+DiracDelta[t-30]+DiracDelta[t-58];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-17]+DiracDelta[t-34]+DiracDelta[t-66];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,t-y] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,t-y] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,0,80},PlotRange->All,MaxRecursion->5, PlotPoints->{20,100}]

p[x_,y_]:=Sin[x]^2/(Cosh[y/6]-Cos[x]);
Print[p[x,y]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-15]+DiracDelta[t-30]+DiracDelta[t-58];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-17]+DiracDelta[t-34]+DiracDelta[t-66];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,t-y] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,t-y] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,0,80},PlotRange->All,MaxRecursion->5, PlotPoints->{20,100}]

p[x_,y_]:=Sin[x]^2/(Cosh[y/8]-Cos[x]);
Print[p[x,y]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-15]+DiracDelta[t-30]+DiracDelta[t-58];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-17]+DiracDelta[t-34]+DiracDelta[t-66];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,t-y] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,t-y] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,0,80},PlotRange->All,MaxRecursion->5, PlotPoints->{20,100}]

Can't really have this depend on y ad t separately in this form. Here is the same thing rewritten as a function of three variables p[x,y,t]

p[x_,y_,t_]:=Sin[x]^2/(Cosh[(y-t)/3]-Cos[x]);
Print[p[x,y,t]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-15]+DiracDelta[t-30]+DiracDelta[t-58];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-17]+DiracDelta[t-34]+DiracDelta[t-66];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,y,t] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,y,t] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,1,80},PlotRange->All,MaxRecursion->5, PlotPoints->{20,200}]

Plot[Tooltip[{Cosh[y],Cosh[y/2],Cosh[y/4],Cosh[y/8]}],{y,-10,10}]

p[x_,y_,t_]:=Sin[x]^2/(Cosh[8 (t-y)/y]-Cos[x]);
Print[p[x,y,t]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-15]+DiracDelta[t-30]+DiracDelta[t-58];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-17]+DiracDelta[t-34]+DiracDelta[t-66];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,y,t] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,y,t] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,1,100},PlotRange->All,MaxRecursion->5, PlotPoints->{20,200}]

p[x_,y_,t_]:=Sin[x]^2/(Cosh[8 (t-y)/y]-Cos[x]);
Print[p[x,y,t]];
sL[t_]:=DiracDelta[t-8]+DiracDelta[t-16]+DiracDelta[t-32]+DiracDelta[t-64];
sR[t_]:=DiracDelta[t-9]+DiracDelta[t-18]+DiracDelta[t-36]+DiracDelta[t-72];
sxy[x_,y_]=1/(2 Pi) Integrate[p[x,y,t] sL[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0]+1/(2 Pi) Integrate[p[Pi-x,y,t] sR[t],{t,-Infinity,Infinity},Assumptions->Pi>x>0];
Plot3D[sxy[x,y],{x,0,Pi},{y,0,80},PlotRange->All,MaxRecursion->5,PlotPoints->{20,200},ImageSize->800]

saved as harm2.png