%《數(shù)字信號(hào)處理教程——MATLAB釋義與實(shí)現(xiàn)》第四章例4.5.2程序hc452
% 連續(xù)信號(hào)頻譜計(jì)算中有相鄰諧振峰的算例
% 電子工業(yè)出版社出版  陳懷琛編著 2004年9月
%
T0=[0.6,0.15,0.15,0.15]; % 各次計(jì)算擬采用的T,編成向量T0
N0=[256,256,256,2048];     % 各次計(jì)算擬采用的N,編成向量N0
for r=1:4  
    T=T0(r);N=N0(r); n=0:N-1;      % 根據(jù)計(jì)算順序選用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，移到對(duì)稱(chēng)位置
    k=floor(-(N-1)/2:(N-1)/2);      % 乃奎斯特頻率下標(biāo)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語(yǔ)句寫(xiě)子圖標(biāo)題
    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