% exa060204_BT.m,  for example 6.2.4 
% to test BT.m  ;

clear all;
Fs=1000;% 采樣頻率
% 產生含有噪聲的序列
var=sqrt(1/exp(1.0));
n=0:1/Fs:1;
N=length(n);
e=var*randn(1,N);
w0=100*pi;
w1=50*pi;
xn=exp(j*w0*n-j*pi)+exp(j*w1*n-j*0.7*pi)+e;

% 繪制信號波形
subplot(311)
plot(n,abs(xn))
xlabel('n')
title('x(n)=exp(j*w0*n-j*pi)+exp(j*w1*n-j*0.7*pi)+e(n)')

% 計算序列的自相關函數
m=-500:500
[r,lag]=xcorr(xn,500,'biased')
subplot(312)
hndl=stem(m,r);
set(hndl,'Marker','.')
set(hndl,'MarkerSize',2);
ylabel('自相關函數R(m)')

% 利用間接法計算功率譜
k=0:1000;
w=(pi/500)*k;
M=k/500;
X=r*(exp(-j*pi/500).^(m'*k));
magX=abs(X);
subplot(313)
plot(M,10*log10(magX));
xlabel('功率譜的BT法估計')