%《數字信號處理教程——MATLAB釋義與實現》第四章例4.3.4a程序hc434a
% 計算無限長序列頻譜時截取長度的自動選擇程序
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
a=1; b=100; beta=1;             % 給定初始數據
while b>beta                    % 判斷是否應結束循環運算
    N1=2^a;n1=0:N1-1;           % 確定數據長度N1
    x1=0.7.^n1; X1=fft(x1);     % 求長度N1的序列x及其FFT X1   
    N2=2*N1; n2=0:N2-1;         % 數據長度加倍為N2=2*N1
    x2=0.7.^n2;X2=fft(x2);      % 求長度N2的序列x2及其FFT X2
    k1p=0:N1/2;k2p=2*k1p;       % 確定兩序列對應點的下標
    d=max(abs(X1(k1p+1)-X2(k2p+1))); % 求對應點上FFT的誤差
    Xm=max(abs(X1(k1p+1)));     % 求X1幅特性的最大值
    b=d/Xm*100;                 % 求相對誤差的百分數
    a=a+1;                      % 序列加長一倍
N2,b                          % 結束循環后顯示達到要求的長度N2和相對誤差b
k=floor(-N2/2+0.5:N2/2-0.5);dw=2*pi/N2;     % 繪圖用的乃奎斯特頻率范圍
subplot(1,2,1) 
plot(k*dw,abs(fftshift(X2))),title('(a)振幅頻譜'),hold on % 畫幅特性
xlabel('頻率'),ylabel('幅度')
set(gcf,'color','w')              % 置圖形背景色為白
subplot(1,2,2) 
plot(k*dw,angle(fftshift(X2))),title('(b)相位頻譜'),hold on % 畫相特性
xlabel('頻率'),ylabel('相位')
line([-5,5],[0,0])
end