%《數字信號處理教程——MATLAB釋義與實現》第四章例4.6.3a程序hc463a
% 用IDFT由已知頻譜求連續時間傅立葉反變換IDTFT的信號
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
for i=1:2
    wc=3;Tmax=0.1*pi/wc      	  % 最大采樣周期Tmax，取算出的臨界數值的1/10
    T=input('(T<Tmax)T= ')        % 用戶選擇并輸入T
    Nmin=100*2*pi/wc/T	        % 最小樣點數目Nmin
    N=input('(N>Nmin)N=  ');    % 用戶選擇并輸入N
    dw=2*pi/(N*T); 		        % 頻譜序列長度及頻率分辨率
    %kn=floor(-(N-1)/2:-1/2); 	% 負頻率下標向量
    %kp=floor(0:(N-1)/2); 		% 正頻率下標向量
    %w=[kp,kn]*dw;     		    % 將負頻率移到正頻率的右方，形成新的頻率排序
    w=[0:N-1]*dw;
    M=floor(wc/dw);             % 有效頻率的邊界下標（整數）
    Xa=[ones(1,M+1),zeros(1,N-2*M-1),ones(1,M)];	% 按新的頻率排序輸入給定的頻譜
    xa=ifft(Xa/T);               % Xa/T是序列頻譜X(k)，其IFFT為數字序列x(n)=xa(n) 
    if i==1 plot([0:N-1]*T,xa),hold on			% 用plot語句繪圖，實現插值
    else plot([0:N-1]*T,xa,':r')
    end
    title('T= ')
end
legend('N=2048','N=4096')       % 圖例標注
set(gcf,'color','w');           % 置圖形背景色為白