信号与系统实验3周期信号频域分析

一、目的

（1）掌握周期信号傅立叶级数分解与合成的计算公式

（2）掌握利用MATLAB实现周期信号傅立叶级数分解与综合方法 （3）理解并掌握周期信号频谱特点

二、周期信号傅立叶级数

周期信号是定义在(,)区间内，按一定时间间隔（周期T）不断重复的信号。可表示为

f(t)f(tmT)

式中m为任意整数，T为周期，周期的倒数成为该信号频率。

实验内容

仿照例程，实现下述周期信号的傅立叶级数分解与合成： f(t)

1

O 4 -4 -3 1

要求：

5 t （a）首先，推导出求解a0，an，bn的公式，计算出前10次系数； （b）利用MATLAB求解a0，an，bn的值，其中an，bn求解前10次系

数，并给出利用这些系数合成的信号波形。

（a）设周期信号f(t)的周期为T1，角频率12f1则该周期信号可以展开成傅立叶级数。 （1）三角形式傅立叶级数

2，且满足狄里赫利条件，T1

f(t)a0a1cos1tb1sin1ta2cos2tb2sin2t...ancosntbnsinnt...a0ancos(n1t)bnsin(n1t)n1n1

1a0T12anT12bnT1T12T12T12T12T1f(t)dtf(t)cosn1tdtf(t)sinn1tdt

2T12

（2）指数形式傅立叶级数

f(t)Fnenjn1t,n0,1,2,3,

1FnT1f(t)eT12T12jnt1dt

（b）求解a0，an，bn及合成信号波形所用程序：

function [A_sym,B_sym]=CTFShchsym

% 采用符号计算求一个周期内连续时间函数f的三角级数展开系数,再用这些 % 展开系数合成连续时间函数f.傅立叶级数 % 函数的输入输出都是数值量 % Nf=6 谐波的阶数 % Nn 输出数据的准确位数

% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 % B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % tao=1 tao/T=0.2 syms t n k x T=4; tao=T/4; a=-1.5; if nargin<4 Nf=10; end

if nargin<5 Nn=32; end

x=time_fun_x(t);

A0=int(x,t,a,T+a)/T; %求出三角函数展开系数A0

As=2/T*int(x*cos(2*pi*n*t/T),t,a,T+a); %求出三角函数展开系数As Bs=2/T*int(x*sin(2*pi*n*t/T),t,a,T+a); %求出三角函数展开系数Bs

A_sym(1)=double(vpa(A0,Nn)); %获取串数组A0所对应的

ASC2码数值数组 for k=1:Nf

A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); %获取串数组A所对应的ASC2码数值数组

B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); %获取串数组B所对应的ASC2码数值数组 end ; if nargout==0 c=A_sym;

disp(c); %输出c为三角级数展开系数:第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 d=B_sym;

disp(d); %输出d为三角级数展开系数: 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 t=-3*T:0.01:3*T;

f0=c(1); %直流 f1=c(2).*cos(2*pi*1*t/T)+d(2).*sin(2*pi*1*t/T); % 基波

f2=c(3).*cos(2*pi*2*t/T)+d(3).*sin(2*pi*2*t/T); % 2次谐波 f3=c(4).*cos(2*pi*3*t/T)+d(4).*sin(2*pi*3*t/T); % 3次谐波 f4=c(5).*cos(2*pi*4*t/T)+d(5).*sin(2*pi*4*t/T); % 4次谐波 f5=c(6).*cos(2*pi*5*t/T)+d(6).*sin(2*pi*5*t/T); % 5次谐波 f6=c(7).*cos(2*pi*6*t/T)+d(7).*sin(2*pi*6*t/T); % 6次谐波 f7=c(8).*cos(2*pi*7*t/T)+d(8).*sin(2*pi*7*t/T); % 7次谐波 f8=c(9).*cos(2*pi*8*t/T)+d(9).*sin(2*pi*8*t/T); % 8次谐波 f9=c(10).*cos(2*pi*9*t/T)+d(10).*sin(2*pi*9*t/T); % 9次谐波 f10=c(11).*cos(2*pi*10*t/T)+d(11).*sin(2*pi*10*t/T); % 10次谐波

f11=f0+f1+f2; % 直流+基波+2次谐波

f12=f11+f3; % 直流+基波+2次谐波+3次谐波

f13=f12+f4+f5+f6; % 直流+基波+2次谐波+3次谐波+4次谐波+5次谐波+6次谐波

f14=f13+f7+f8+f9+f10; %0~10次 subplot(2,2,1)

plot(t,f0+f1),hold on

y=time_fun_e(t); %调用连续时间函数-周期矩形脉冲 plot(t,y,'r:')

title('直流+基波') axis([-8,8,-0.5,1.5]) subplot(2,2,2) plot(t,f12),hold on y=time_fun_e(t); plot(t,y,'r:')

title('1-3次谐波+直流') axis([-8,8,-0.5,1.5])

subplot(2,2,3) plot(t,f13),hold on y=time_fun_e(t); plot(t,y,'r:')

title('1-6次谐波+直流') axis([-8,8,-0.5,1.5]) subplot(2,2,4) plot(t,f14),hold on y=time_fun_e(t); plot(t,y,'r:')

title('1-10次谐波+直流') axis([-8,8,-0.5,1.5]) hold off end

function y=time_fun_e(t)

% 该函数是CTFShchsym.m的子函它由符号函数和表达式写成 a=1.5; T=4; h=1; tao=T/4;

t=-3*T:0.01:3*T;

e1=1/2+1/2.*sign(t-0.5+tao/2); e2=1/2+1/2.*sign(t-0.5-tao/2);

y=h.*(e1-e2); %连续时间函数-周期矩形脉冲

function x=time_fun_x(t)

% 该函数是CTFShchsym.m的子函数。它由符号变量和表达式写成。 h=1;

x1=sym('Heaviside(t)')*h; x=x1-sym('Heaviside(t-1)')*h;

源程序修改：

function [A_sym,B_sym]=CTFShchsym

% 采用符号计算求一个周期内连续时间函数f的三角级数展开系数,再用这些 % 展开系数合成连续时间函数f.傅立叶级数 % 函数的输入输出都是数值量 % Nf=6 谐波的阶数 % Nn 输出数据的准确位数

% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 % B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % tao=1 tao/T=0.2 syms t n k x T=5;

tao=0.2*T; a=0.5; if nargin<4 Nf=6; end

if nargin<5 Nn=32; end

x=time_fun_x(t);

A0=int(x,t,-a,T-a)/T; %求出三角函

数展开系数A0

As=2/T*int(x*cos(2*pi*n*t/T),t,-a,T-a); %求出三角函数展开系数As Bs=2/T*int(x*sin(2*pi*n*t/T),t,-a,T-a); %求出三角函数展开系数Bs

A_sym(1)=double(vpa(A0,Nn)); %获取串数组A0所对应的ASC2码数值数组 for k=1:Nf

A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); %获取串数组A所对应的ASC2码数值数组

B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); %获取串数组B所对应的ASC2码数值数组 end if nargout==0 c=A_sym;

disp(c) %输出c为三角级数展开系数:第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 d=B_sym;

disp(d) %输出d为三角级数展开系数: 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 t=-8*a:0.01:T-a;

f0=c(1); %直流 f1=c(2).*cos(2*pi*1*t/5)+d(2).*sin(2*pi*1*t/5); % 基波

f2=c(3).*cos(2*pi*2*t/5)+d(3).*sin(2*pi*2*t/5); % 2次谐波 f3=c(4).*cos(2*pi*3*t/5)+d(4).*sin(2*pi*3*t/5); % 3次谐波 f4=c(5).*cos(2*pi*4*t/5)+d(5).*sin(2*pi*4*t/5); % 4次谐波 f5=c(6).*cos(2*pi*5*t/5)+d(6).*sin(2*pi*5*t/5); % 5次谐波 f6=c(7).*cos(2*pi*6*t/5)+d(7).*sin(2*pi*6*t/5); % 6次谐波 f7=f0+f1+f2; % 直流+基波+2次谐波

f8=f7+f3; % 直流+基波+2次谐波+3次谐波

f9=f8+f4+f6; % 直流+基波+2次谐波+3次谐波+4次谐波+6次谐波 subplot(2,2,1)

plot(t,f0+f1),hold on

y=time_fun_e(t); %调用连续时间函数-周期矩形脉冲 plot(t,y,'r:')

title('周期矩形波的形成—直流+基波') axis([-4,4.5,-0.5,1.5]) subplot(2,2,2) plot(t,f7),hold on y=time_fun_e(t); plot(t,y,'r:')

title('周期矩形波的形成—直流+基波+2次谐波') axis([-4,4.5,-0.5,1.5]) subplot(2,2,3) plot(t,f8),hold on y=time_fun_e(t);

plot(t,y,'r:')

title('直流+基波+2次谐波+3次谐波') axis([-4,4.5,-0.5,1.5]) subplot(2,2,4) plot(t,f9),hold on y=time_fun_e(t); plot(t,y,'r:')

title('基波+2次谐波+3次谐波+4次谐波+6次谐波') axis([-4,4.5,-0.5,1.5]) end

function y=time_fun_e(t)

% 该函数是CTFShchsym.m的子函它由符号函数和表达式写成 a=0.5; T=5; h=1;

tao=0.2*T;

t=-8*a:0.01:T-a;

e1=1/2+1/2.*sign(t+tao/2); e2=1/2+1/2.*sign(t-tao/2);

y=h.*(e1-e2); %连续时间函数-周期矩形脉冲

function x=time_fun_x(t)

% 该函数是CTFShchsym.m的子函数。它由符号变量和表达式写成。 h=1;

x1=sym('Heaviside(t+0.5)')*h; x=x1-sym('Heaviside(t-0.5)')*h;

三、周期信号频谱分析

对周期为T1 的信号f(t)进行傅立叶级数展开可得到

f(t)Fnenjn1ta0ancos(n1t)bnsin(n1t)

n1n1其中

12 T1FnFnejn1(anjbn) （3-8） 2如果求出an和bn，根据以下两式可以画出周期信号的幅度谱Fn~和相位谱n~。

Fn12anbn2 幅频 2bn 相频 anarctann实验内容

已知周期为T=4的三角波，在第一周期（-2function [A_sym,B_sym]=CTFSshbpsym(T,Nf)

% 采用符号计算求[0,T]内时间函数的三角级数展开系数。 % 函数的输入输出都是数值量 % Nn 输出数据的准确位数

% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 % B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % T T=m*tao， 信号周期 % Nf 谐波的阶数

% m (m=T/tao)周期与脉冲宽度之比,如m=4,8,16,100等 % tao 脉宽:tao=T/m syms t n y

if nargin<3

Nf=input('please Input 所需展开的最高谐波次数:Nf='); end

T=input('please Input 信号的周期T='); if nargin<5 Nn=32; end

y=time_fun_s(t);

A0=2/T*int(y,t,0,T/2);

As=2/T*int(y*cos(2*pi*n*t/T),t,0,T/2); Bs=2/T*int(y*sin(2*pi*n*t/T),t,0,T/2); A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf

A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); end

if nargout==0

An=fliplr(A_sym); %对A_sym阵左右对称交换 An(1,k+1)=A_sym(1); %A_sym的1*k阵扩展为1*(k+1)阵

An=fliplr(An); %对扩展后的S1阵左右对称交换回原位置 Bn=fliplr(B_sym); %对B_sym阵左右对称交换

Bn(1,k+1)=0; %B_sym的1*k阵扩展为1*(k+1)阵

Bn=fliplr(Bn); %对扩展后的S3阵左右对称交换回原位置

FnR=An/2-i*Bn/2; % 用三角函数展开系数A、B值合成付里叶指数系数 FnL=fliplr(FnR); N=Nf*2*pi/T; k2=-N:2*pi/T:N;

Fn=[FnL,FnR(2:end)];

%subplot(3,3,3)

%x=time_fun_e(t); % 调用连续时间函数-周期矩形脉冲

subplot(2,1,1)

stem(k2,abs(Fn)); %画出周期矩形脉冲的频谱（T=M*tao） title('连续时间函数周期三角波脉冲的双边幅度谱') axis([-80,80,0,0.12])

line([-80,80],[0,0],'color','r') line([0,0],[0,0.12],'color','r') end

function x=time_fun_e(t)

% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 % t 是时间数组 % T 是周期 duty=tao/T=0.2 T=5;

t=-2*T:0.01:2*T; tao=T/5;

x=rectpuls(t,tao); %产生一个宽度tao=1的矩形脉冲 subplot(2,2,2) plot(t,x) hold on

x=rectpuls(t-5,tao); %产生一个宽度tao=1的矩形脉,中心位置在t=5处 plot(t,x) hold on

x=rectpuls(t+5,tao); %产生一个宽度tao=1的矩形脉,中心位置在t=-5处 plot(t,x)

title('周期为T=5，脉宽tao=1的矩形脉冲') axis([-10,10,0,1.2])

function y=time_fun_s(t) syms t y=1-abs(t);

x1=sym('Heaviside(t+2)'); x=x1-sym('Heaviside(t-2)'); y=y*x;

ezplot(t,y,[-10,10]) grid

源程序修改：

function [A_sym,B_sym]=CTFSshbpsym(T,Nf)

% 采用符号计算求[0,T]内时间函数的三角级数展开系数。 % 函数的输入输出都是数值量 % Nn 输出数据的准确位数

% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 % B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % T T=m*tao， 信号周期 % Nf 谐波的阶数

% m (m=T/tao)周期与脉冲宽度之比,如m=4,8,16,100等 % tao 脉宽:tao=T/m syms t n y if nargin<3

Nf=input('please Input 所需展开的最高谐波次数:Nf='); end

T=input('please Input 信号的周期T='); if nargin<5 Nn=32; end

y=time_fun_s(t); A0=2/T*int(y,t,0,T);

As=2/T*int(y*cos(2*pi*n*t/T),t,0,T); Bs=2/T*int(y*sin(2*pi*n*t/T),t,0,T); A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf

A_sym(k+1)=double(vpa(subs(As,n,k),Nn));

B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); end

if nargout==0

An=fliplr(A_sym); %对A_sym阵左右对称交换 An(1,k+1)=A_sym(1); %A_sym的1*k阵扩展为1*(k+1)阵

An=fliplr(An); %对扩展后的S1阵左右对称交换回原位置 Bn=fliplr(B_sym); %对B_sym阵左右对称交换 Bn(1,k+1)=0; %B_sym的1*k阵扩展为1*(k+1)阵

Bn=fliplr(Bn); %对扩展后的S3阵左右对称交换回原位置

FnR=An/2-i*Bn/2; % 用三角函数展开系数A、B值合成付里叶指数系数 FnL=fliplr(FnR); N=Nf*2*pi/T; k2=-N:2*pi/T:N;

Fn=[FnL,FnR(2:end)];

subplot(3,3,3)

x=time_fun_e(t); % 调用连续时间函数-周期矩形脉冲

subplot(2,1,1)

stem(k2,abs(Fn)); %画出周期矩形脉冲的频谱（T=M*tao） title('连续时间函数周期矩形脉冲的双边幅度谱') axis([-80,80,0,0.12]) line([-80,80],[0,0]) line([0,0],[0,0.12]) end

function x=time_fun_e(t)

% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 % t 是时间数组

% T 是周期 duty=tao/T=0.2 T=5;

t=-2*T:0.01:2*T; tao=T/5;

x=rectpuls(t,tao); %产生一个宽度tao=1的矩形脉冲 subplot(2,2,2) plot(t,x) hold on

x=rectpuls(t-5,tao); %产生一个宽度tao=1的矩形脉,中心位置在t=5处 plot(t,x) hold on

x=rectpuls(t+5,tao); %产生一个宽度tao=1的矩形脉,中心位置在t=-5处 plot(t,x)

title('周期为T=5，脉宽tao=1的矩形脉冲') axis([-10,10,0,1.2])

function y=time_fun_s(t)

% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 syms a a1

T=input('please Input 信号的周期T='); M=input('周期与脉冲宽度之比M='); A=1; tao=T/M; a=tao/2;

y1=sym('Heaviside(t+a1)')*A; y=y1-sym('Heaviside(t-a1)')*A; y=subs(y,a1,a); y=simple(y);