close all;
N=8192;
n=0:N-1;
L2=128;
L3=64;
% u1=[-1:1/L2:1];
% dh=0.54+0.46*cos(pi*u1);
% ds=1/pi*abs(sin(pi*u1))+(1-abs(u1)).*cos(pi*u1);
% dr=ones(1,length(u1));
f1=0.12;
fai1=2*pi*rand;
f2=0.30;
fai2=2*pi*rand;
f3=f1+f2;
fai3=fai1+fai2;
f4=0.19;
fai4=2*pi*rand;
f5=0.17;
fai5=2*pi*rand;
f6=f4+f5;
fai6=fai4+fai5;
s1=cos(2*pi*f1*n);%+fai1);
s2=cos(2*pi*f2*n);%+fai2);
s3=cos(2*pi*f3*n);%+fai3);
s4=cos(2*pi*f4*n);%+fai4);
s5=cos(2*pi*f5*n);%+fai5);
s6=cos(2*pi*f6*n);%+fai6);
x=s1+s2+s3+s4+s5+s6;
figure,
plot(abs(fft(x)));
title('頻譜圖');
% y=awgn(x,-50);
figure,
% subplot(2,2,1)
window=rectwin(L2);
% [Bspec,waxis]=bispecd(x2,56,window,256,50);
[Bspec,waxis]=bispeci(x,16);
title('加parzen窗的雙譜圖');
figure,
[Bspec,waxis]=bispeci(x,32);
title('加parzen窗窗長加大的雙譜圖');
% subplot(2,2,2)
figure,
window=hamming(L2);
% [Bspec,waxis]=bispecd(x,256,window,256,50);
[Bspec,waxis]=bispeci(x,16,256,50,'unbiased',256,1);
title('加矩形窗的雙譜圖');
% subplot(2,2,3)
figure,
window=parzenwin(L2);
% [Bspec,waxis]=bispecd(x,20,window,256,50);
[Bspec,waxis]=bispeci(x,32,256,50,'unbiased',256,1);
title('加矩形窗窗長加大的雙譜圖');
figure,
subplot(3,2,1)
wind=boxcar(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加矩形窗的功率譜圖');
subplot(3,2,2)
wind=triang(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加三角形窗的功率譜圖');
subplot(3,2,3)
wind=hamming(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加漢明窗的功率譜圖');
subplot(3,2,4)
wind=hann(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加漢寧窗的功率譜圖');
subplot(3,2,5)
wind=parzenwin(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加parzen窗的功率譜圖');
wind=kaiser(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
subplot(3,2,6)
plot(w,pxy);
axis([0,0.5,0,10]);
title('加最優窗的功率譜圖');

figure,
subplot(3,2,1)
wind=boxcar(L2);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加矩形窗窗長加大的功率譜圖');
subplot(3,2,2)
wind=triang(L2);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加三角形窗窗長加大的功率譜圖');
subplot(3,2,3)
wind=hamming(L2);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加漢明窗窗長加大的功率譜圖');
subplot(3,2,4)
wind=hann(L2);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加漢寧窗窗長加大的功率譜圖');
subplot(3,2,5)
wind=parzenwin(L2);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
axis([0,0.5,0,10]);
title('加parzen窗窗長加大的功率譜圖');
wind=kaiser(L2);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
subplot(3,2,6)
plot(w,pxy);
axis([0,0.5,0,10]);
title('加最優窗窗長加大的功率譜圖');
% axis([-0.5,0.5,-0.5,0.5]);
%  bspec_abs=abs(Bspec);
% plot(bspec_abs);
%   contour(waxis,waxis,abs(Bspec),4),grid on 
% title('guji');