%《數字信號處理教程——MATLAB釋義與實現》第九章例9.5.3程序hc953
% 窗函數對頻譜分析的影響
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
N=128;
fr=[0.25,0.28];
R0=[128,256];
for i=1:2                               % 計算兩組參數下的結果
    R=R0(i);
    n=0:N-1;
    x=0.5*sin(2*pi*n*fr(1))+sin(2*pi*n*fr(2));          % 列出信號序列
    X=fft(x,R);k=0:R-1;                                 % 求離散頻譜
    subplot(3,2,2*i-1),stem(k,abs(X),'.');grid on       % 畫頻譜圖
    title(['(a) N=',num2str(N),',R=',num2str(R),',矩形窗DFT'])
    ylabel('|X(k)|')
    x1=x.*hamming(N)';                                  % 對信號序列加哈明窗
    X1=fft(x1,R);k=0:R-1;                               % 求離散頻譜
    subplot(3,2,2*i),stem(k,abs(X1),'.');grid on        % 畫頻譜圖
    title(['(b) N=',num2str(N),',R=',num2str(R),',哈明窗DFT'])
    ylabel('|X(k)|')
end
subplot(3,2,5),plot(k/R,abs(X)),grid on                 % 畫X的連續頻譜
title(['(e) N=',num2str(N),',R=',num2str(R),',矩形窗DTFT'])
xlabel('f=k/R'),ylabel('|X(f)|')
subplot(3,2,6),plot(k/R,abs(X1)),grid on                % 畫X1的連續頻譜
title(['(f) N=',num2str(N),',R=',num2str(R),',哈明窗DTFT'])
xlabel('f=k/R'),ylabel('|X(f)|')
set(gcf,'color','w')                                % 置圖形背景色為白
subplot(3,2,1),axis([29,39,0,70])                   % 設定比例尺
subplot(3,2,2),axis([29,39,0,40])
subplot(3,2,3),axis([58,78,0,70])
subplot(3,2,4),axis([58,78,0,40])
subplot(3,2,5),axis([0.24,0.3,0,60])
subplot(3,2,6),axis([0.24,0.3,0,40])