function [fhsmooth]=fmultitapersmooth(t,M,NP,L,X)
ja=0;r=sqrt(-1); syms lamd1; syms lamd2;
[w1,w2]=meshgrid(0:pi/NP:pi,0:pi/NP:pi);
for i=1:t
K=i;
H1=(K/t);
if (K-t/2)
h(1,i,1)=1-6.*((abs(H1)).^2-(abs(H1)).^3);
else
h(1,i,1)=2.*((1-abs(H1)).^3);
end
AI=(i)/t;
CI=AI*pi;
h(1,i,2)=sin(CI)/CI;
h(1,i,3)=.54+.46*cos(((i).*pi)/M);
end
for k=1:3
for l=1:3
for m=1:3
qq=0;
for i=1:t
qq=qq+(h(:,i,k)*h(:,i,l)*h(:,i,m));
Q(k,l,m)=qq;
end
ja=ja+(Q(k,l,m)).^2;
end
end
end
mu(1)=0;
NP1=NP+1;
mu(NP1)=1;
NP3=2*NP+1;
for I=2:NP3
mu(I)=(I-1)./NP;
mu(I);
end
for i1=1:NP1
for i2=i1:NP1
I1=(mu(i1)-lamd1)/(2*M);
I2=(mu(i2)-lamd2)/(2*M);
I3=(mu(i1)+mu(i2)-lamd1-lamd2)/(2*M);
ker1=(sin (I1)/I1).^2;
ker2=(sin (I2)/I2).^2;
ker3=(sin (I3)/I3).^2;
mker=ker1*ker2*ker3;
Kernel=mker/((2*pi).^3);
for j=1:L
result1=0;result11=0;result2=0;result22=0;result3=0;result33=0;
for k=1:3
for l=1:3
for m=1:3
tap1=0;tap2=0;tap3=0;tap11=0;tap22=0;tap33=0;
for i=1:t
tap1=tap1+sum(h(:,i,k)*X(:,i,j)*cos(lamd1*i));
tap11=tap11+sum(h(:,i,k)*X(:,i,j)*sin(lamd1*i));
tap2=tap2+sum(h(:,i,l)*X(:,i,j)*cos(lamd2*i));
tap22=tap22+sum(h(:,i,l)*X(:,i,j)*sin(lamd2*i));
tap3=tap3+sum(h(:,i,m)*X(:,i,j)*cos((lamd1+lamd2)*i));
tap33=tap33+sum(h(:,i,m)*X(:,i,j)*sin((lamd1+lamd2)*i));
end
Q(k,l,m);
result1=result1+Q(k,l,m)*tap1;
result11=result11+Q(k,l,m)*tap11;
result2=result2+Q(k,l,m)*tap2;
result22=result22+Q(k,l,m)*tap22;
result3=result3+Q(k,l,m)*tap3;
result33=result33+Q(k,l,m)*tap33;
end
end
end
sreal=result1*result2*result3+result1*result22*result33+result11*result2*result33-result11*result22*result3;
simag=result1*result2*result33-result1*result22*result3-result11*result2*result3-result11*result22*result33;
Imulti(j)=sqrt(sreal.^2+simag.^2)/(4*pi^2*ja);
mkImu = matlabFunction(Kernel.*Imulti(j));
fjjsmooth(i1,i2,j)=integral2(mkImu ,-pi,pi,-pi,pi);
end
fhsmooth(i1,i2)=sum(fjjsmooth(i1,i2,:))/(M.^2*L);
fhsmooth(i2,i1)=fhsmooth(i1,i2);
end
end
mesh(w1,w2,fhsmooth);
xlabel('w1');
ylabel('w2');
zlabel('Estimated Bispectrum modulus');
end
