close all;
clear all;
clc;

N=64;
L3=10;
L2=32;
n=0:N-1;
w1=2*pi*0.2;
w2=2*pi*0.3;
ss1=(0.8.^n);
ss2=(0.7.^n);
s1=ss1.*cos(w1*n);
s2=ss2.*cos(w2*n);
x=s1+s2;
figure,
plot(s1);
hold on;
plot(s2,'r')

figure,
subplot(3,2,1)
wind=boxcar(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
title('加矩形窗的功率譜圖');
subplot(3,2,2)
wind=triang(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
title('加三角形窗的功率譜圖');
subplot(3,2,3)
wind=hamming(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
title('加漢明窗的功率譜圖');
subplot(3,2,4)
wind=hann(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
title('加漢寧窗的功率譜圖');
subplot(3,2,5)
wind=parzenwin(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
plot(w,pxy);
title('加parzen窗的功率譜圖');
wind=kaiser(L3);
[pxy,w]=pwelch(x,wind,0,256);
w=w/(2*pi);
subplot(3,2,6)
plot(w,pxy);
title('加最優窗的功率譜圖');

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

% [Bspec,waxis]=bispecd(x2,56,window,256,50);
[Bspec,waxis]=bispeci(x,10);
title('加parzen窗的雙譜圖');
figure,
[Bspec,waxis]=bispeci(x,15);
title('加parzen窗窗長加大的雙譜圖');
figure,

% [Bspec,waxis]=bispecd(x,256,window,256,50);
[Bspec,waxis]=bispeci(x,10,256,50,'unbiased',256,1);
title('加矩形窗的雙譜圖');
figure,
window=parzenwin(L2);
% [Bspec,waxis]=bispecd(x,20,window,256,50);
[Bspec,waxis]=bispeci(x,15,256,50,'unbiased',256,1);
title('加矩形窗窗長加大的雙譜圖');