%對高斯白噪聲n(t)中的信號x(t)=cos(2*pi*60*t)+cos(2*pi*63*t)+0.1*n(t),分別采用周期圖法、Welch法、Burg法、
%Music法進行譜估計。其中fs=1000Hz,nfft=1024和2048，加窗類型自選。要求（1)畫出x(t)的波形圖（2）畫出各種方法的譜估計
%圖（3）對譜估計結果進行分析。
fs = 1000;   
t = 0:1/fs:1; 
s = cos(2*pi*60*t) + cos(2*pi*63*t);
w = 0.1*randn(size(t)); 
x = s + w;  
figure(1); plot(t,x); title('含噪聲信號')
xlabel('t/s');
nfft1=1024; nfft2=2048;
figure(2); subplot(2,1,1); hp1=spectrum.periodogram('hamming'); psd(hp1,x,'nfft',nfft1,'fs',fs); title('周期圖譜估計1024點');
           subplot(2,1,2); hp2=spectrum.periodogram('hamming'); psd(hp2,x,'nfft',nfft2,'fs',fs); title('周期圖譜估計2048點');
figure(3); subplot(2,1,1); hw1 = spectrum.welch('hamming',64,50); psd(hw1,x,'nfft',nfft1,'fs',fs); title('welch譜估計1024點');
           subplot(2,1,2); hw2 = spectrum.welch('hamming',64,50); psd(hw2,x,'nfft',nfft2,'fs',fs); title('welch譜估計2048點');
figure(4); subplot(2,1,1); pburg(x,50,nfft1,fs); title('burg譜估計1024點');
           subplot(2,1,2); pburg(x,50,nfft2,fs); title('burg譜估計2048點');
figure(5); subplot(2,1,1); pmusic(x,2,nfft1,fs); title('music譜估計1024點');
           subplot(2,1,2); pmusic(x,2,nfft2,fs); title('music譜估計2048點');