function[w1,w2,fn,fargu,fmod]=bispectwPARZENEWINDmod(N,M,NP,X)
[w1,w2]=meshgrid(0:pi/20:pi,0:pi/20:pi);
A=pi/M;
for i=1:M+1
K=i-1;
H1=(K/M);
if(K-M/2)
W(i)=1-6*(H1^2-H1^3);
else
W(i)=2*(1-H1)^3;
end
end
for i=1:M+1
SUMV=0;
for j=1:N-i+1
SUMV=SUMV+(X(j)-mean(X)).*(X(i+j-1)-mean(X));
end
R(i)=(SUMV)/N;
end
for j2=1:M+1
for j1=j2:M+1
sum=0;
for K=1:N-j1+1
K1=j1-1+K;
K2=j2-1+K;
sum=sum+(X(K1)-mean(X)).*(X(K2)-mean(X))'.*(X(K)-mean(X));
end
C(j1,j2)=(sum/N);
end
end
w(1)=0;
NP1=NP+1;
w(NP1)=1;
NP3=2*NP+1;
for i=2:NP3
w(i)=(i-1)./NP;
w(i);
end
for i1=1:NP1
for i2=i1:NP1
S2=R(1);s3=R(1);s4=R(1);
for i=1:M-1
w(i1)
w(i2)
S2=S2+W(i+1).*2.*(R(i+1).*cos(w(i1)*pi.*i));
s3=s3+W(i+1).*2.*(R(i+1).*cos(w(i2)*pi.*i));
s4=s4+W(i+1).*2.*(R(i+1).*cos((w(i1)+w(i2))*pi.*i));
end
fhat(i1)=(S2)./(2*pi);
fhat(i2)=(s3)./(2*pi);
fhat(i1+i2)=(s4)./(2*pi);
end
end
for i1=1:NP1
for i2=i1:NP1
ARF=0;
AIF=0;
for j2=2:M
for j1=j2+1:M+1
A1=(j1-1)*pi;
A2=(j2-1)*pi;
A12=A1-A2;
R=W(j1).*W(j2).*W(j1-j2).*C(j1,j2);
ARF=ARF+R*(cos(A1*w(i1)+A2*w(i2))+cos(A1*w(i2)+A2*w(i1))+cos(A1*w(i1)+A12*w(i2))+cos(A1*w(i2)+A12*w(i1))+cos(A12*w(i1)-A2*w(i2))+cos(A12*w(i2)-A2*w(i1)));
AIF=AIF-R*(sin(A1*w(i1)+A2*w(i2))+sin(A1*w(i2)+A2*w(i1))-sin(A1*w(i1)+A12*w(i2))-sin(A1*w(i2)+A12*w(i1))+sin(A12*w(i1)-A2*w(i2))+sin(A12*w(i2)-A2*w(i1)));
end
end
for i=2:M+1
AI=(i-1)*pi;
ARF=ARF+(C(i,1)+C(i,i)).*(W(i).^2).*(cos(AI*w(i1))+cos(AI*(w(i2)))+cos(AI*(w(i2)+w(i1))));
AIF=AIF+(C(i,1)-C(i,i)).*(W(i).^2).*(-sin(AI*w(i1))-sin(AI*(w(i2)))+sin(AI*(w(i1)+w(i2))));
end
ARF=ARF+C(1,1);
fmod(i1,i2)=(sqrt(ARF.^2+AIF.^2)).*(1/(4*pi^2));
fmod(i2,i1)=fmod(i1,i2);
fargu(i1,i2)=atan(AIF./ARF);
fn(i1,i2)=abs(fmod(i2,i1)./(sqrt(fhat(i1).*fhat(i2).*fhat(i1+i2-1))));
fn(i2,i1)=fn(i1,i2);
end
end
mesh(w1,w2,fn);
mesh(w1,w2,fmod);
end
