%==========================================================================
%                      傅里葉變換      程序名：CTFT
%  本程序是采用數值積分的方法計算實信號的傅里葉變換。
%  對信號的要求是：
%  1、信號是時限信號，或者信號在時間t大于某一數值時，信號的值已經衰減到等于
%     或接近于零；
%  2、信號的頻譜中不能含有沖激成分，因為計算機無法表現和處理沖激信號。
%==========================================================================
%                            下段程序為初始化
%==========================================================================
clear;close all
p=0.01;
t=-5:p:5;
x=u(t+0.5)-u(t-0.5);
W1=6*pi;                           %確定信號的頻率范圍
N=512;           
k=-N:N;                              %對信號的角頻率進行采樣的點數為2N+1點
W=k*W1/N;                            %計算頻率的采樣間隔
%==========================================================================
%    下段程序計算信號的傅里葉變換，包括幅度譜和相位譜計算
%==========================================================================
X=x*exp(-j*t'*W)*p;                  %進行傅里葉數值積分運算（矩陣運算）
X1=abs(X);                           %求模——幅度譜
phi=57.3*angle(X);                     %求相位譜
%==========================================================================
%                            下段程序繪制圖形
%==========================================================================
subplot(2,1,1);plot(t,x);grid on;
xlabel('t');ylabel('The curve of x(t)');
title('The signal x(t)');
subplot(2,2,3);plot(W/(2*pi),X);grid on;
axis([-W1/(2*pi),W1/(2*pi),-0.3*max(X1),1.1*max(X1)])
title('The magnitude spectra ');
xlabel('w*2pi (radians/s)');
subplot(2,2,4);
plot(W/(2*pi),phi,'r');grid on;
xlabel('w*2pi (radians/s)');title('The phase spectra');
axis([-W1/(2*pi),W1/(2*pi),-180,180])
%===========================================================================