试验三周期信号频域分析

1、实验三周期信号频域分析一、目的(1)掌握周期信号傅立叶级数分解与合成的计算公式(2)掌握利用MATLAB实现周期信号傅立叶级数分解与综合方法(3)理解并掌握周期信号频谱特点二、周期信号傅立叶级数周期信号是定义在(*,g)区间内，按一定时间问隔(周期 T)不断重复的信号。可表小为f(t) u f (t mT)式中m为任意整数，T为周期，周期的倒数成为该信号频率。1、连续时间周期信号的分解设有周期信号f(t),周期为T,角频率。1 = 2肛=至，且满足狄里赫利条件，则该Ti周期信号可以展开成傅立叶级数，即可表示为一系列不同频率的正弦或复指数信号之和。傅立叶级数有三角形式和指数形式两种：1)三角形式

2、傅立叶级数三角形式傅立叶级数为f (t) = a。 a1 c o s h s i n it a2cos2t b? s i nt . ： a c o st bn s i n nt .0000二a。一： a” c o si( J) 一: bn si m(4)n七式中系数按如下公式求解n1 T1=-21 f (t)dtT1 2 T1=T1 f (t)cosn ；tdt T2 T1=：i f (t) sin n、tdtT12)指数形式傅立叶级数指数形式的傅立叶级数表达式为：f(t)=EFne阳 n=0,1 苫津3,n -：式中E称为傅立叶复系数，可由下式求得T1Fn = I f (t)exdtT1(3

3、-1)(3-2)(3-3)(3-4)(3-5)(3-6)傅立叶级数的指数形式和三角形式是等价的，具系数可相互转换，具体参见信号与系 统教科书。2、连续时间周期信号的傅立叶综合任何满足狄里赫利条件的周期信号，可以表示成式(3-1)或(3-5)的和式形式, 式(3-1)或(3-5)成为连续时间周期信号(CTFS)综合公式。一般说来，傅立叶级数系数有无限个非零值，即任何具有有限个间断点的周期信号 都一定有一个无限项非零系数的傅立叶级数表示。但对于数值计算来说，这是无法实现 的。在实际应用中，可以用有限项傅立叶级数求和来逼近。即：15f (t) = Fd8 =3 + an cos(nt) + bn s

4、in(n缸t)(3-7)n _Nn 1_n 1当N值取得较大时，上去就是原周期看号f(t)的一不很好近似。式(3-7)常称做f(t) 的截断傅立叶级数表示。MATLAB的符号积分函数int()可以用来求解连续时间周期信号的截断傅立叶级数 及傅立叶表示。求积函数int()的具体使用格式如下：intf=int(f,v);给出符号表达式f对指定变量v的(不带积分常数)不定积分；intf=int(f,v,a,b);给出符号表达式f对指定变量v的定积分；3、利用MATLAB实现周期信号的傅立叶级数分解与综合(1)利用MATLAB求解周期矩形脉冲傅立叶级数，并绘制出各次谐波叠加的傅立叶综 合波形图。周期矩

5、形脉冲为f(t)=京6 -nT),式中t=1,T=5。 n =二二采用三角形式傅立叶级数分解与综合形式，用式(3-2)(3-4)求出傅立叶级数分 解系数，运用MATLAB的符号运算功能，用式(3-7)实现信号的综合，谐波的阶数Nf =6。(a)实现流程利用MATLAB实现上述分析过程的流程如下：编写子函数x=time_fun_x(t),用符号表达式表示出周期信号在第一个周期内的 符号表达式，并赋值返回给符号变量 x;编写子函数y=time_fun_e(t),求出该周期信号在绘图区间内的信号样值，并赋 值给返回变量y;编写求解信号傅立叶系数及绘制合成波形图的通用CTFShchsym.m,该函数流

6、程如下： 调用函数time_fun_x(t),获取周期信号的符号表达式；求出信号的傅立叶索数；求出各次谐波；绘制各次谐波叠加波形图； 调用函数time_fun_e(t),绘制原信号波形图。(b) MATLAB算法提示及近明采用符号积分int求一个周期内时间函数的三角级数展开系数：& = Ao ,& =As, h =Bs,即计算式(3-2)(3-4)的值；用循环语句foL-end求出三角级数展开系数a, bn的数值，分别为 A_sym,B_sym;用disp()语句输出三角级数展开系数A_sym, B_sym；用傅立叶三角级数展开式(3-7)合成自续时间信号；化简表达式，据函数奇偶性可知，若 f

7、(t)为奇函数，则 a=0;若f(t)为偶函 数，M bn =0。(c)源程序编写函数文件CTFShchsym.m,这是一个计算连续时间周期信号的三角级数前6次展开系数，再用这6次谐波合成原连续时间周期信号的程序，如下所示。function A_sym,B_sym=CTFShchsym%采用符号计算柬二个周期内连续时间函数 f的三角级数展开系数，再用这些% 展开系数合成连续时间函数f.傅立叶级数16% 函数的输入输出都是数值量% Nf=6谐波的阶数% Nn输出数据的准确位数% A_sym第1元素是直流项，其后元素依次是 1,2,3次谐波cos项展开系数% B_sym第2,3,4,.元素依次是1

8、,2,3次谐波sin项展开系数% tao=1 tao/T=0.2 syms t n k x T=5;tao=0.2*T;a=0.5;if nargin4;Nf=6;endif nargin5;Nn=32;end x=time_fun_x(t);A0=int(x,t,-a,T-a)/T;%求出三角函数展开系数A0As=int(2*x*cos(2*pi*n*t/T)/T,t,-a,T-a);%求出三角函数展开系数 AsBs=int(2*x*sin(2*pi*n*t/T)/T,t,-a,T-a);% 求出三角函数展开系数 BsA_sym(1)=double(vpa(A0,Nn);%获取用数组 A0

9、所对应的 ASC2桎值数组for k=1:NfA_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码数R数组 end if nargout=0c=A_sym;disp(c)%输出c为三角级数展开系数：第1元素是直流项，其后元素依次是1,2,a.次谐波cos项展开系数d=B_sym;disp(d) %输出d为三角级数展开系数：第2,3,4,.元素依次是1,2,3次谐波sin项展开系数t=-8*a:0.01:T-a

10、;f1=c(1)+c(2).*cos(2*pi*1*t/5)+0.*sin(2*pi*1*t/5); ;% 基波f2=c(3).*cos(2*pi*2*t/5)+0.*sin(2*pi*2*t/5); ;% 2 次t皆波f3=c(4).*cos(2*pi*3*t/5)+0.*sin(2*pi*3*t/5);% 3 次 t皆波f4=c(5).*cos(2*pi*4*t/5)+0.*sin(2*pi*4*t/5); ;% 4 次t皆波f5=c(6).*cos(2*pi*6*t/5)+0.*sin(2*pi*6*t/5);% 6 次 t皆波f6=f1+f2;% 基波+2次谐波f7=f6+f3;%基波

11、+2次谐波+3次谐波f8=f7+f4+f5;%基波+2次谐波+3次谐波+4次谐波+6次谐波subplot(2,2,1)plot(t,f1),hold ony=time_fun_e(t)%调用连续时间函数-周期矩形脉冲plot(t,y,r:)title(周期矩形波的形成一基波+直流)axis(-4,4.5,-0.5,1.5)subplot(2,2,2)plot(t,f6),hold on y=time_fun_e(t)17plot(t,y,r:)t田e(周期矩形波的形成一基波+2次谐波)axis(-4,4.5,-0.5,1,5)subplot(2,2,3)plot(t,f7),hold ony=

12、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,f8),hold ony=time_fun_e(t)plot(t,y,r:)title(基波+2次谐波+3次谐波+4次谐波+6次谐波)axis(-4,4.5,-0.5,1.5)end%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;%func

13、tion 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);%连续时间函数-周期矩形脉冲(d)程序运行结果在MATLAB命令窗口键入CTFShchsym,并按回车键，即可绘制出周期矩形波形 各次谐波合成波形，如图3-1所示。18执行如下命令就可以得到三角级数展开系数：A_sym,B_sym=CTFShchsym结杀为一A_sym =- 0.20

14、000.37420.30270.20180.09350.0000-0.0624B_sym =0000000(2)实验内容, f(t)n(a)(b)首先，推导出求解a。，an, b 利用MATLAB求解a , a ,要求：仿照例程，实现下述周期信号的傅立叶级数分解与合成：的公式，计算出前10次系数；bn的值，其中bn求解前10次系数，并给出利用这些系数合成的信号波形。三、周期信号频谱分析1、周期信号频谱对周期为Ti的信号f(t)进行傅立叶级数展开可得到f(t) =Fnen 二：二00cO二 a0 八 a” cos(nt)八 bn sin(nt)其中2 二 i = 一Ti_ 曲 1Fn = Fn

15、e 4=2 (an - jbn)(3-8)如果求出a和bn ,根据以下两式可以画出周期信号的幅度谱|Fn切和相位谱中n6。1,一一Fn =一；耳+幅频2Q = -arctan2相频an2、MATLAB 实现由上述分析可知，只要求出了周期信号的傅立叶级数幅度和相位Q ,就可以根据Fn和中n随频率8二门町的变化关系画出幅度谱和相位谱。现仍以周期矩形脉冲为f。)=京0骐-nT),(7=1,T=5)为例，来说明如何用MATLAB绘制周期信号的频谱图，并对周期信号的频谱特性进行分析。19由于绘制频谱图的前提是必须先求出周期信号的傅立叶级数系数，因此只需对上述给出的求周期信号傅立叶级数的函数 CTFShc

16、hsym.m进行适当修改，即可编写出绘制周 期信号频谱的通用函数。需注意的是，由于周期信号的频谱是离散的，故在绘制频谱时，采用的是stem命令而不是plot命令，下面给出实现上述通用程序 CTFStpshsym.m的过程。(1)实现流程此处采用三角形式傅立叶级数分解，求出分解系数an和bn,再利用式(3-8)求出傅立叶复指数系数Fn,画出Fn的幅度谱Fn 8。谐波的阶数Nf可任意指定，此处指定Nf =60。(a)实现流程如下：编写子函数y=time_fun_s(t),用符号表达式表示出周期信号在第一个周期内的符号表达式，并赋值返回给符号变量 V；编写子函数x=time_fun_e(t),求出该

17、周期信号在绘图区间内的信号样值，并赋 值给返回变量x;编写求解信号傅立叶复指数系数 Fn的通用函数，该函数流程如下： 调用函数time_fun_s(t),获取周期信号的符号表达式；求出信号的三帚形氐的傅立叶系数an和bn;求出信号的傅立复指数系数Fn；绘制Fn的幅度频谱图； 调用函数time_fun_e(t)绘制信号波形。(b)算法提示采用符号积分int求0,T内时间信号的三角级数展开系数：4 = A0 , a = As,bn = Bs;用循环语句foL-end求出三角级数展开系数an和bn的数值，并赋值给变量A_sym(k+1), B_sym(k+1)；心角级数展开家数a和bn得到复指数展开

18、系数Fn： Fn可以根据式(38)求出, 需要注意的是an、bn和Fn的自变量取值情况， 耳、R的变量h的取值范围为 n=0,1,2,3，N，而匕的变量n的取值范围为n = 0,1,2,3,，土N ,为了从出和 bn得到E,需要用到MATLAB的反折函数fliplr来实现频谱的反折；例如已知0, 1, 2, 3, 4。要组成员序列反折后再-4, -3, -2, -1, 0, 1, 2, 3, 4)。程序一个序列a=5 4 3 2 1,采样位置为: 叠加原序列的一个新序列(采样位置： 如下：n=0:4;a=5 4 3 2 1;subplot(2,1,1);stem(n,a);b=fliplr(a

19、);k=-4:4;c=b,a(2:end);subplot(2,1,2);stem(k,c);结果如右图所示。20(c)源程序代码编写函数文件CTFStpshsym.m,如下所示：CTFStpshsym.mfunction 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谐波的阶数%Nn输出数

20、据的准确位数% m (m=T/tao)周期与脉冲宽度之比，如m=4,8,16,100等%tao脉宽:tao=T/msyms t n yNf=input(please Input所需展开的最高谐波次数：Nf=);T=input(please Input 信号的周期 T=);Nn=32;y=time_fun_s(t);A0=2*int(y,t,0,T)/T;As=int(2*y*cos(2*pi*n*t/T)/T,t,0,T);Bs=int(2*y*sin(2*pi*n*t/T)/T,t,0,T);A_sym(1)=double(vpa(A0,Nn); for k=1:NfA_sym(k+1)=d

21、ouble(vpa(subs(As,n,k),Nn);B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn); endS1=fliplr(A_sym);%对A_sym阵左右对称交换S1(1,k+1)=A_sym(1); %A_sym 的 1*k 阵扩展为 1*(k+1)阵S2=fliplr(1/2*S1);%对扩展后的S1阵左右对称交换回原位置S3=fliplr(1/2*B_sym); %对 B_sym 阵左右对称交换S3(1,k+1)=0;%B_sym 的 1*k 阵扩展为 1*(k+1)阵S4=fliplr(S3);%3扩展后的S3阵左右对称交换回原位置S5=S2-i

22、*S4;%用三角函数展开系数 A、B值合成付里叶指数系数S6=fliplr(S5);N=Nf*2*pi/T;k2=-N:2*pi/T:N;S7=S6,S5(2:end);subplot(3,3,3)x=time_fun_e1(t) %调用连续时间函数-周期矩形脉冲stem(k2,abs(S7); %画出周期矩形脉冲的频谱(T=M*tao)title(连续时间函数周期矩形脉冲的双边幅度谱)axis(-80,80,0,0.12)line(-80,80,0,0)line(0,0,0,0.12)21%function y=time_fun_s(t,T)%该函数是CTFSshbpsym.m的子函数。它由

23、符号变量和表达式写成。syms a a1T=input(please Input 信号的周期 T=);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);%function x=time_fun_e1(t)%该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。% t是时间数组% T是周期 duty=tao/T=0.2T=5;t=-2*T:0.01:2*T;tao=T/5;x=rectpuls(t,1); %产生一个宽度tao=1的矩

24、形脉冲plot(t,x)hold onx=rectpuls(t-5,1); %产生一个宽度tao=1的矩形脉，中心位置在t=5处 plot(t,x) hold onx=rectpuls(t+5,1); %产生一个宽度tao=1的矩形脉，中心位置在t=-5处 plot(t,x)title(周期为T=5,脉宽tao=1的矩形脉冲) axis(-10,10,0,1.2)(d)运行结果及分析程序运行结果调用CTFStpshsym.m函数文件，即可绘出周期矩形脉冲波形信号的双边幅度频谱。指令如下：在 MATLAB 命令窗口键入 CTFStpshsym,回车，命令窗口将会出现：please Input所需展开的最高