一、實驗目的

1、理解傅立葉變換的MATLAB實現方法；
2、掌握系統頻率響應特性的計算方法和特性曲線的繪製方法，理解具有不同頻率響應特性的濾波器對訊號的濾波作用。
3、掌握用MATLAB語言進行系統頻響特性分析的方法。

二、實驗原理

1、傅立葉變換的MATLAB實現

訊號f(t)的傅立葉變換與反變換公式為：

傅立葉變換存在的充分條件是:

在MATLAB中有兩種實現傅立葉變換的方法，一種是利用MATLAB中的專用函式直接求解，另一種是傅立葉變換的數值計算實現法。

(1) 呼叫專用函式實現

MATLAB中實現傅立葉變換的函式為：

F=fourier( f ) 對f(t)進行傅立葉變換，其結果為F(w)
F=fourier(f,v) 對f(t)進行傅立葉變換，其結果為F(v)
F=fourier( f,u,v ) 對f(u)進行傅立葉變換，其結果為F(v)

MATLAB中實現傅立葉反變換的函式

f =ifourier( F ) 對F(w)進行傅立葉反變換，其結果為f(t)
f=ifourier(F,u) 對F(w)進行傅立葉反變換，其結果為f(u)
f=ifourier( F,v,u ) 對F(v)進行傅立葉反變換，其結果為f(u)
注意：
（1）在呼叫函式fourier( )及ifourier( )之前，要用syms命令將所有需要用到的變數（如t,u,v,w）定義成符號變數，或用符號定義符sym將其定義為符號表示式。
（2）採用fourier( )及fourier( )得到的返回函式，仍然為符號表示式。在對其作圖時要用ezplot( )函式，而不能用plot()函式。
（3）fourier( )及fourier( )函式的應用有很多侷限性，如果在返回函式中含有δ(ω)等函式，則ezplot( )函式無法作出圖。用fourier( )函式對某些訊號進行變換時，其返回函式如果包含一些不能直接表達的式子，也無法作圖。另外，在很多場合，f(t)不能表示成符號表示式，此時只能用數值計演算法進行傅氏變換。
例1：求 的傅立葉變換
matlab程式如下

clear all
syms t
F=fourier(exp(-2*abs(t)))

例2：求門函式f(t)=ε(t+1)-ε(t-1)的傅立葉變換，並畫出幅度頻譜圖。
matlab程式如下

clear all
t0=-3;t1=3;t=t0:0.02:t1;                        % 定義時間範圍
w0=-6*pi;w1=6*pi;w=w0:0.02:w1;               % 定義頻率範圍
f=sym('heaviside(t+1)-heaviside(t-1)')             % 定義符號函式f(t)
F=fourier(f)                                  % 求f(t)的傅立葉變換
F=simple(F)                                  % 化簡F(jw)的表示式
f1=subs(f,t,'t');                                % 將t陣列代入f(t)後用f1表示
fmin=min(f1)-0.2;fmax=max(f1)+0.2;             % 求f1的最大和最小值
Fv=subs(F,w,'w');                             % 將w陣列代入F(jw)後用Fv表示
F1=abs(Fv);                                  % 求F(jw)的模
P1=angle(Fv);                                % 求F(jw)的相角
subplot(3,1,1),plot(t,f1);                        % 在第一幅圖上畫f(t)
ylabel('G2(t)'); grid on;
axis([t0,t1,fmin,fmax]);
Fmin=min(F1)-0.05; Fmax=max(F1)+0.05;
subplot(3,1,2);
plot(w,F1,'color','k');                            % 在第二幅圖上畫|F(jw)|
ylabel('|F(jw)|');grid on;
axis([w0,w1,Fmin,Fmax]);
subplot(3,1,3); 
plot(w,P1*180/pi,'color','k');                      % 在第三幅圖上畫相位頻譜
ylabel('相位(度)'); grid on;

執行程式結果圖如下

(2) 傅立葉變換的數值計算實現法

數值計算方法實現連續時間訊號的傅立葉變換，實質上是藉助於MATLAB的強大數值計算功能進行的一種近似計算。其實現原理如下：
對於連續時間訊號f(t)，其傅立葉變換為：

其中τ為取樣間隔，如果f(t)是時限訊號，或者當|t|大於某個給定值時，f(t)的值已經衰減得很厲害，可以近似地看成是時限訊號，則上式中的n取值就是有限的，假定為N，有：

若對頻率變數ω進行取樣，得：

通常取：
，其中ω0是要取的頻率範圍，或訊號的頻頻寬度。採用matlab
實現上式時，其要點是要生成f(t)的N個樣本值f (nτ)的向量，以及向量
，兩向量的內積（即兩矩陣的乘積）完成上式的傅立葉變換的數值計算。

注意：

時間取樣間隔τ的確定，依據是τ 必須小於奈奎斯特（Nyquist）取樣間隔。如果f(t)不是嚴格的帶限訊號，則根據實際計算的精度要求來確定一個適當的頻率ω0為訊號的頻寬。
例3：用數值計演算法實現上面門函式f(t)=ε(t+1)-ε(t-1)的傅立葉變換，並畫出幅度頻譜圖。
分析：該訊號的頻譜F(jω)=2Sa(ω)，第一個過零點頻率為π，一般將此頻率認為是訊號的頻寬。但考慮到F(jω)的形狀（為抽樣函式），假如將精度提高到該值的50倍，即取ω0=50ωB，F0=ω0 /2π則據此確定的取樣間隔為：τ≤1/2F0 = 0.02
matlab程式如下

clear all
R=0.02; %取樣間隔τ=0.02 
t = -2:R:2; % t為從-2到2，間隔為0.02的行向量,有201個樣本點
ft=[zeros(1,50),ones(1,101),zeros(1,50)]; % 產生f(t)的樣值矩陣（即f(t)的樣本值組成的行向量）
W1=10*pi; %取要計算的頻率範圍
M=500; k=-M:M; w=k*W1/M; %頻域取樣數為M, w為頻率正半軸的取樣點
Fw=f t*exp(-j*t'*w)*R; %求傅氏變換F(jw) 
FRw=abs(Fw); %取振幅
subplot(2,1,1) ; plot(t,ft) ;grid;  %畫出原始函式f(t)的波形
xlabel('t') ; ylabel('f(t)'); title('f(t)=u(t+1)-u(t-1)'); 
subplot(2,1,2); plot(w,FRw) ;grid on; %畫出振幅頻譜的波形
xlabel ('w') ; ylabel ('F(w)'); 

執行程式結果如下

2、用matlab計算連續時間LTI系統的頻率響應

頻率特性是指系統在正弦訊號激勵下的穩態響應隨頻率變化的情況，包括響應的幅度隨頻率的變化情況和響應的相位隨頻率的變化情況兩個方面。

f (t)、y(t)分別為系統的激勵和響應，h(t)是系統的單位衝激響應，三者的關係如下：
時域：
頻域：
或：

H(jω)為系統的頻域數學模型，它實際上就是系統的單位衝激響應h(t)的傅立葉變換。即：

其中，׀H(jω)׀為幅度頻率相應，反映訊號經過系統之後，訊號各頻率分量的幅度變化情況，為相位特性，反映訊號經過系統後，訊號各頻率分量在相位變換情況。
當系統的激勵，則其響應為

若輸入訊號為正弦訊號，即f(t) = sin(0t)，則系統響應為

系統對某一頻率分量的影響體現為，訊號的幅度被׀H(jω)׀加權，訊號的相位被移相。
關鍵是確定系統的頻率響應。
matlab裡面求系統頻率響應的函式為freqs（），該函式可求出系統頻率響應的數值解，並可繪出系統的幅頻及相頻響應曲線。
當系統的頻率響應H(jω)是jω的有理多項式時，

matlab裡面求系統頻率響應的函式為freqs，該函式可直接計算系統的頻率響應的數值解，並可繪出系統的幅頻及相頻響應曲線。其呼叫格式如下
（1）h=freqs(b,a,w) 該呼叫格式中，b和a分別對應於（7-7）式的向量[b1,b2,…,bm]和[a1,a2,…,an]，w為形如w1:p:w2的冒號運算定義的系統頻率響應的頻率範圍，w1為頻率起始值，w2為頻率終止值，p 為頻率取樣間隔。向量h為返回在向量w所定義的頻率點上，系統頻率響應的樣值。
（2）[h,w]=freqs(b,a)
該呼叫格式將計算預設頻率範圍內200個頻率點的系統頻率響應的樣值，並賦值給返回變數h，200個頻率點記錄在w中。
（3）[h,w]=freqs(b,a,n)
該呼叫格式將計算預設頻率範圍內n個頻率點上系統頻率響應的樣值，並賦值給返回變數h，n個頻率點記錄在w中。
（4）freqs(b,a)
該格式並不返回系統頻率響應的樣值，而是以對數座標的方式繪出系統的幅頻響應和相頻響應曲線。
舉例說明如下：
a=[1 2 1]; b=[0 1]; h=freqs(b,a,0:0.5:2*pi)；%計算0～2π頻率範圍內以間隔0.5取樣的系統頻率響應的樣值；
例4：某連續時間LTI系統的微分方程如下

編寫MATLAB程式繪製系統的幅度響應、相位響應、頻率響應的實部和頻率響應的虛部。
matlab程式如下

clear  all
b=[1];      
a=[1 3 2];   
[H,w]=freqs(b,a);  
Hm=abs(H);          
phai=angle(H);    
Hr=real(H);         
Hi = imag(H);  
subplot(2,2,1)     
plot(w,Hm), grid on,  title('Magnitude response'), 
xlabel('Frequency in rad/sec')
subplot(2,2,3)
plot(w,phai); grid on; title('Phase response'); xlabel('Frequency in rad/sec')
subplot(2,2,2)
plot(w,Hr); grid on; title('Real part of frequency response'),  
xlabel('Frequency in rad/sec')
subplot(2,2,4)
plot(w,Hi); grid on; title('Imaginary part of frequency response'),  
xlabel('Frequency in rad/sec')

執行結果如下：

例5:已知一RC電路如圖所示，系統的輸入電壓為f(t)，輸出訊號為電阻兩端的電壓y(t)，當RC=0.04，f(t)=cos5t+cos100t, 試求該系統的零狀態響應y(t)

由圖知其頻率響應為

餘弦訊號cos(ω0t)通過LTI系統的零狀態響應為

matlab程式如下

clear all
RC=0.04;
t=linspace(-2,2,1024);
w1=5;w2=100;
H1=j*w1/(j*w1+1/RC);
H2=j*w2/(j*w2+1/RC);
f=cos(5*t)+cos(100*t);
y=abs(H1)*cos(w1*t+angle(H1))+ abs(H2)*cos(w2*t+angle(H2));
subplot(2,1,1);
plot(t,f); ylabel('f(t)'); xlabel('Time(s)');
subplot(2,1,2);
plot(t,y); ylabel('y(t)'); xlabel('Time(s)');

三、實驗內容

1.已知三個LTI系統的微分方程如下

系統1：
系統2：

系統3：

用matlab繪製幅度響應、相位響應、頻率響應的實部和頻率響應的虛部曲線圖。
並分析系統具有什麼濾波特性？

系統1：

程式如下：

clear  all
b=[1];      
a=[1 1 25];  
[H,w]=freqs(b,a);  
Hm=abs(H);          
phai=angle(H);    
Hr=real(H);         
Hi = imag(H);  
subplot(2,2,1)     
plot(w,Hm), grid on,  title('·ù¶ÈÏìÓ¦'), 
xlabel('Frequency in rad/sec')
 
subplot(2,2,3)
plot(w,phai); grid on; title('ÏàλÏìÓ¦'); xlabel('Frequency in rad/sec')
 
subplot(2,2,2)
plot(w,Hr); grid on; title('頻率響應的實部'),  
xlabel('Frequency in rad/sec')
 
subplot(2,2,4)
plot(w,Hi); grid on; title('頻率響應的虛部'),  
xlabel('Frequency in rad/sec')

執行結果如下：

系統2：

程式如下：

clear  all
b=[1 1];      
a=[1 1];  
[H,w]=freqs(b,a);  
Hm=abs(H);          
phai=angle(H);    
Hr=real(H);         
Hi = imag(H);  
subplot(2,2,1)     
plot(w,Hm), grid on,  title('幅度響應'), 
xlabel('Frequency in rad/sec')
 
subplot(2,2,3)
plot(w,phai); grid on; title('相位響應'); xlabel('Frequency in rad/sec')
 
subplot(2,2,2)
plot(w,Hr); grid on; title('頻率響應的實部'),  
xlabel('Frequency in rad/sec')
 
subplot(2,2,4)
plot(w,Hi); grid on; title('頻率響應的虛部'),  
xlabel('Frequency in rad/sec')

執行結果如下：

系統3：

程式如下：**

clear  all
b=[262];      
a=[10 48 148 306 401 262];  
[H,w]=freqs(b,a);  
Hm=abs(H);          
phai=angle(H);    
Hr=real(H);         
Hi = imag(H);  
subplot(2,2,1)     
plot(w,Hm), grid on,  title('幅度響應'), 
xlabel('Frequency in rad/sec')
 
subplot(2,2,3)
plot(w,phai); grid on; title('相位響應'); xlabel('Frequency in rad/sec')
 
subplot(2,2,2)
plot(w,Hr); grid on; title('頻率響應的實部'),  
xlabel('Frequency in rad/sec')
 
subplot(2,2,4)
plot(w,Hi); grid on; title('頻率響應的虛部'),  
xlabel('Frequency in rad/sec')

執行結果如下：

2. 分別用傅立葉變換的數值計演算法和呼叫專用函式編寫MATLAB程式，計算

f(t)= 3e-2tε(t)的傅立葉變換。
程式如下：

clear all
syms t
F=fourier(3*exp(-2*t)*heaviside(t))

執行結果如下：

F =
 
3/(2 + w*1i)

3.已知一RC電路如圖所示 系統的輸入電壓為f(t)，輸出訊號為電阻兩端的電壓y(t).當R1=10kΩ , R2=30kΩ C=100µF，f(t)=cos50t+1, 試求該系統的零狀態響應y(t)

程式如下：

clear all
RC=1;
t=linspace(-2,2,1024);
w1=50;
H1=j*w1/(j*w1+1/RC);
 
f=cos(50*t)+1;
y=abs(H1)*cos(w1*t+angle(H1))
subplot(2,1,1);
plot(t,f); ylabel('f(t)'); xlabel('Time(s)');
subplot(2,1,2);
plot(t,y); ylabel('y(t)'); xlabel('Time(s)');

執行結果如下：

四、 實驗報告要求

1.簡述實驗目的和實驗原理；

2.寫出其對應的matlab程式；

3.計算相應的衝激響應、零狀態響應及卷積積分的理論值，並與實驗結果進行比較。

4.上機除錯程式的方法及實驗中的心得體會。

本次上級實驗，我掌握了連續時間訊號與系統的頻域分析方法，從頻域的角度對訊號與系統的特性進行分析；掌握了連續時間訊號傅立葉變換與傅立葉逆變換的實現方法；掌握了連續時間傅立葉變換的特點及應用；掌握了連續時間傅立葉變換的數值計算方法及繪製訊號頻譜的方法。熟悉了
（1）.syms命令:定義符號變數
呼叫格式:
symsx，定義x為符號變數，可以直接使用。
symsxy,定義x和y為符號變數，可以直接使用。
（2）.ezplot函式:實現一元函式的繪圖。相比plot函式要制定自變數範圍，ezplot 無需
資料準備，可以直接畫圖，尤其適用於符號函式。