%《數字信號處理教程——MATLAB釋義與實現》第四章例4.5.2程序hc452
% 連續信號頻譜計算中有相鄰諧振峰的算例
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
T0=[0.6,0.15,0. % 各次計算擬采用的T,編成向量T0
N0=[256,256,256,2048];     % 各次計算擬采用的N,編成向量N0
for r=1:415,0.15];  
    T=T0(r);N=N0(r); n=0:N-1;      % 根據計算順序選用T及N
    DOmega=2*pi/(N*T);                % 頻率分辨率
    xa=exp(-0.01*n*T).*cos(2*n*T)+2*exp(-0.015*n*T).*cos(2.1*n*T); %樣本序列xa
    Xa=T*fftshift(fft(xa));Xa(1)       % 求xa的FFT，移到對稱位置
    k=floor(-(N-1)/2:(N-1)/2);      % 乃奎斯特頻率下標DOmega=2*pi/(N*T);
    set(gcf,'color','w')              % 置圖形背景色為白
    subplot(2,2,r),plot(k*DOmega,abs(Xa))    % 繪制幅頻特性
    xlabel('\Omega'),ylabel('|Xa(\Omega)|')
    switch r                                % 用switch語句寫子圖標題
    case 1, title('(a)T=0.6,N=256')
    case 2, title('(b)T=0.15,N=256')
    case 3, title('(c)T=0.15,N=256,頻率軸放大')
    case 4, title('(d)T=0.15,N=2048,頻率軸放大')
    otherwise 
    end    
    if r>2 axis([1,3,0,60]), end            % 不同r定不同比例尺
end