MATLAB信號與系統分析（五）——連續時間信號的頻譜分析

一、實驗目的：

1、掌握傅立葉級數（FS），學會分析連續時間周期信號的頻譜分析及MATLAB實現；

2、掌握傅立葉變換（FT），了解傅立葉變換的性質以及MATLAB實現。

 

二、利用符號運算求傅里葉級數的系數

1、復習幾個函數：

F1=int(f,v,a,b) — 對f表達式的v變量在（a,b)區間求定積分
F2=subs(s,OLD,NEW)－用新變量NEW代替S中的指定變量OLD。
F3=vpa(x,n) ： 顯示可變精度計算；x為符號變量，n表示要精確計算的位數。

 

2、周期函數的傅里葉級數的形式

 

3、利用符號運算求傅立葉級數的系數

代碼抄襲如下：

%ex_1
%求系數
clear all; syms t x n t0;
T=10;    % 信號周期
tao_2=0.5; %脈沖寬度
Nf=7; % 分解的最高級數
Nn=6; % 有效位數

x=heaviside(t+t0)-heaviside(t-t0); % 注意：Ver 2011b
x=subs(x,t0,tao_2)  

A0=int(x,t,-tao_2,T-tao_2)/T % a0

As=int(x*2*cos(2*pi*n*t/T)/T,t,-tao_2,T-tao_2) % an

Bs=int(x*2*sin(2*pi*n*t/T)/T,t,-tao_2,T-tao_2) % bn

Fn=(As-j*Bs)/2 % fn

A(1)=double(vpa(A0,Nn)); % 直流分量
for k=1:Nf
     A(k+1)=double(vpa(subs(As,n,k),Nn));
     B(k+1)=double(vpa(subs(Bs,n,k),Nn));
end
A
B
%求各次諧波：：開始采用數值處理方法
t1=-T/2:0.01:T/2;
f0=A(1);                                             %直流分量
f1=A(2).*cos(2*pi*1*t1/T)+B(2).*sin(2*pi*1*t1/T); ;      % 基波
f2=A(3).*cos(2*pi*2*t1/T)+B(3).*sin(2*pi*2*t1/T); ;                 % 2次諧波
f3=A(4).*cos(2*pi*3*t1/T)+B(4).*sin(2*pi*3*t1/T);                 % 3次諧波
f4=A(5).*cos(2*pi*4*t1/T)+B(5).*sin(2*pi*4*t1/T); ;               % 4次諧波
f5=A(6).*cos(2*pi*5*t1/T)+B(6).*sin(2*pi*5*t1/T);                % 5次諧波
f6=A(7).*cos(2*pi*6*t1/T)+B(7).*sin(2*pi*6*t1/T);               % 6次諧波  
f7=f1+f2+f0;                      % 基波+2次諧波
f8=f7+f3+f0;                     % 基波+2次諧波+3次諧波
f9=f8+f4+f6+f0;                 % 基波+2次諧波+3次諧波+4次諧波+6次諧波 
%畫出諧波圖形
y=subs(x,t,t1);         %調用連續時間函數-周期矩形脈沖
subplot(2,2,1),plot(t1,f1),hold on;plot(t1,y,'r:');title('周期矩形波的形成—基波'),axis([-2.5,2.5,-0.5,1.1])
subplot(2,2,2),plot(t1,f7),hold on;plot(t1,y,'r:');title('周期矩形波的形成—基波+2次諧波'),axis([-2.5,2.5,-0.5,1.1])
subplot(2,2,3),plot(t1,f8),hold on;plot(t1,y,'r:');title('基波+2次諧波+3次諧波'),axis([-2.5,2.5,-0.5,1.1])
subplot(2,2,4),plot(t1,f9),hold on;plot(t1,y,'r:');title('基波+2次諧波+3次諧波+4次諧波+6次諧波'),axis([-2.5,2.5,-0.5,1.1])

%求頻譜---------------這邊給出雙邊譜的正軸部分，僅供參考
Fn=(As-j*Bs)/2
Nf=60;
A(1)=double(vpa(A0,Nn));
for k=1:Nf
     A(k+1)=double(vpa(subs(As,n,k),Nn));
     B(k+1)=double(vpa(subs(Bs,n,k),Nn));
end
Fn1(1) = A(1);
Fn1(2:Nf+1)=(A(2:Nf+1)-j*B(2:Nf+1))/2;
absFn1 = abs(Fn1);
angleFn1 = angle(Fn1);
figure;
subplot(211);stem([0:Nf],absFn1,'r.');title('雙邊幅度譜的正半軸，注意與p200圖比較幅度值')
subplot(212);stem([0:Nf],angleFn1,'r.');title('雙邊相位譜的正半軸')

注意：課件里的基波定義有誤，不應該含有直流分量

 

 

% 例9_1
% 觀察周期方波信號的分解與合成
% m：傅里葉級數展開的項數
clc;clear all;close all;

display('Please input the value of m (傅里葉級數展開的項數)'); % 在命令窗口顯示提示信息
m = input('m = ');                                             % 鍵盤輸入傅里葉級數展開的項數
t = -2*pi:0.01:2*pi;                                            % 時域波形的時間范圍-2π~2π，采樣間隔0.01
n = round(length(t)/4);                                        % 根據周期方波信號的周期，計算1/2周期的數據點數
f = [ones(n,1);-1*ones(n,1);ones(n,1);-1*ones(n+1,1)];         %構造周期方波信號
y = zeros(m+1,max(size(t)));
y(m+1,:) = f'; 
figure(1);
plot(t/pi,y(m+1,:));                                           %繪制方波信號
grid;                                                          %在圖形中加入柵格
axis([-2 2 -1.5 1.5]);                                         %指定圖形顯示的橫坐標范圍和縱坐標范圍
title('周期方波');                                             %給顯示的圖形加上標題
xlabel('單位pi','Fontsize', 8);                                %顯示橫坐標單位
x = zeros(size(t));
kk = '1';
for k=1:2:2*m-1                                                %循環顯示諧波疊加圖形
    pause;
    x = x+sin(k*t)/k;
    y((k+1)/2,:) = 4/pi*x;                                     %計算各次諧波疊加和
    plot(t/pi,y(m+1,:));
    hold on;
    plot(t/pi,y((k+1)/2,:));                                   %繪制諧波疊加信號
    hold off;
    grid;
    axis([-2 2 -1.5 1.5]);
    title(strcat('第',kk,'次諧波疊加'));
    xlabel('單位pi','Fontsize', 8);
    kk = strcat(kk,'、',num2str(k+2));
end
pause;
plot(t/pi,y(1:m+1,:));
grid;
axis([-2 2 -1.5 1.5]);
title('各次諧波疊加波形');
xlabel('單位pi','Fontsize', 8);
% End

 

4、周期信號的頻譜分析：

%ex_3求頻譜
clear all;syms t x n t0;
T=5;tao=0.5;Nf=60;Nn=6;
x=heaviside(t+t0)-heaviside(t-t0)
x=subs(x,t0,tao)
A0=int(x,t,-tao,T-tao)/T
As=int(x*2*cos(2*pi*n*t/T)/T,t,-tao,T-tao)
Bs=int(x*2*sin(2*pi*n*t/T)/T,t,-tao,T-tao)
Fn=(As-j*Bs)/2
A(1)=double(vpa(A0,Nn));
for k=1:Nf
     A(k+1)=double(vpa(subs(As,n,k),Nn));
     B(k+1)=double(vpa(subs(Bs,n,k),Nn));
end
t1=-2.5:0.003:2.5;
y=subs(x,t,t1);
subplot(3,1,1),plot(t1,y),title('矩形脈沖'),axis([-2.5,2.5,-0.1,1.1])
%單邊譜
Fs(1)=A(1);
Fs(2:Nf+1)=abs(A(2:Nf+1)-j.*B(2:Nf+1));
Ns=0:Nf;
subplot(3,1,2),stem(Ns,Fs),title('連續時間函數的單邊譜')
%雙邊譜
N=Nf*2*pi/T;
K=-N:2*pi/T:N;
Fs(2:Nf+1)=abs(A(2:Nf+1)-j.*B(2:Nf+1))/2;
Fd=[fliplr(Fs),Fs(2:Nf+1)];
subplot(3,1,3),stem(K,Fd),title('連續時間函數的雙邊譜')

 

5、利用MATLAB實現典型周期信號的頻譜

（1）周期方波脈沖頻譜的Matlab實現

% 周期方波信號頻譜分析
function  CTFS_SQ
%  繪制並觀察周期方波信號頻譜   
%  Nf：傅里葉級數展開的項數
%  an：各次諧波余弦項系數
display('Please input the value of Nf ');
Nf = input('Nf = ');
an = zeros(Nf+1,1);
cn(1) = 0;
for i = 1:Nf 
    an(i+1) = 2/(i*pi)*sin(i*pi/2);                     %計算系數an
    cn(i+1) = abs(2/(i*pi)*sin(i*pi/2));                %計算幅度譜
end
t = -5:0.01:5;
x = square(pi*(t+0.5));                                 %用square函數求方波信號
subplot(211);
plot(t,x);                                              %繪制方波信號
axis([-5 5 -1.5 1.5]);
title('周期方波信號','Fontsize',8);
xlabel('t   (單位:s)', 'Fontsize',8);
subplot(212);
k = 0:Nf;
stem(k,cn);
hold on;
plot(k,cn);
title('幅度頻譜','Fontsize',8);
xlabel('諧波次數', 'Fontsize',8);
% End

 

（2）周期鋸齒波脈沖頻譜的Matlab實現

% 周期鋸齒脈沖信號頻譜分析
function  CTFS_JC
% 繪制並觀察周期鋸齒脈沖信號頻譜特性
%     Nf：諧波的次數
% bn：第1,2,3,...次諧波正弦項系數
display('Please input the value of Nf ');
Nf = input('Nf = ');
bn = zeros(Nf+1,1);
cn = zeros(Nf+1,1);
bn(1) = 0;
for i = 1:Nf
    bn(i+1) = (-1)^(i+1)*1/(i*pi);                          % 計算系數bn
    cn(i+1) = abs(bn(i+1));                                 %計算幅度頻譜
end
t = -5:0.01:5;
x = sawtooth(pi*(t+1));                                     %用sawtooth函數構造周期鋸齒脈沖信號
subplot(211);
plot(t,x);
axis([-5 5 -1.5 1.5]);
title('周期鋸齒脈沖信號','Fontsize',8);
xlabel('t   (單位:s)', 'Fontsize',8);
subplot(212);
k = 0:Nf;
stem(k,cn);
hold on;
plot(k,cn);
title('幅度頻譜','Fontsize',8);
xlabel('諧波次數', 'Fontsize',8);
% End

 

（3）周期三角波脈沖頻譜的Matlab實現

% 周期三角脈沖信號頻譜分析
function  CTFS_SJ
%  繪制周期三角脈沖信號頻譜
%     Nf：諧波的次數
%  an：第1,2,3,...次諧波余弦項系數
display('Please input the value of Nf ');
Nf = input('Nf = ');
an = zeros(Nf+1,1);
cn = zeros(Nf+1,1);
an(1) = 1/2;
an(1) = an(1);
for i = 1:Nf
    an(i+1) = 4*sin(i*pi/2)*sin(i*pi/2)/(i*i*pi*pi);
    cn(i+1) = abs(an(i+1));
end
t = -5:0.01:5;
x = (sawtooth(pi*(t+1),0.5)+1)/2;                        %構造三角脈沖信號
subplot(211);
plot(t,x);
axis([-5 5 -1.5 1.5]);
title('周期三角脈沖信號','Fontsize',8);
xlabel('t   (單位:s)', 'Fontsize',8);
subplot(212);
k = 0:Nf;
stem(k,cn);
hold on;
plot(k,cn);
title('幅度頻譜','Fontsize',8);
xlabel('諧波次數', 'Fontsize',8);
% End

 

三、非周期函數的傅立葉變換

1、利用符號函數求傅立葉變換

傅立葉變換：F=fourier(f);        F,f應為符號表達式
反傅立葉變換：f=ifourier(F);

 

2、連續時間信號傅立葉變換的數值計算

代碼實現

%ex_4
clear all;
R=0.01;
t=-2:R:2;
f=Heaviside(t+1)-Heaviside(t-1);
W1=2*pi*5;
N=1000;k=-N:N;W=k*W1/N;

F=f*exp(-j*t'*W)*R;
F_r=real(F);
figure(1)
subplot(2,1,1);plot(t,f);
xlabel('t');ylabel('f(t)');title('f(t)=u(t+1)-u(t-1)');
subplot(2,1,2);plot(W,F_r);
xlabel('w');ylabel('F(w)');title('f(t)的付氏變換F(w)');
F_A=abs(F);%幅頻特性
F_P=angle(F);%相頻特性
figure(2)
subplot(2,1,1),plot(W,F_A),xlabel('w');ylabel('abs(F(w))');title('f(t)幅頻特性)');
subplot(2,1,2),plot(W,F_P),xlabel('w');ylabel('angle(F(w))');title('f(t)相頻特性)');