MATLAB计算方法和技巧8-1方波

#x(t)x(t)=(8.1.1)(8.1.2)方波的傅立叶分析［问题］周期性方波如B8.1a图所示，对方波分解为傅立叶级数并画出方波合成曲线。［数学模型］方法一：用解析式。一个复杂的振动是由一系列不同频率的简谐振动合成的，把一个复杂振动分解为不同频率的简谐振动的方法称为傅立叶分析。周期性方波可表示为A,(0<t<T/2)0,(T/2<t<T)根据傅立叶定理，周期为T的函数x(t)可以表示为余弦函数和正弦函数之和，设最小圆频率为①，则其他圆频率是最小频率的整数倍，即x(t)=厶+(acosnet+bsinnet)2nnn=122a=cosnetdt，b=sinnetdtnTnT00级数中的第一项a0/2是x(t)的时间平均值。周期与圆频率之间的关系为t=2n®3称为基频，n®称为n次倍频或n次谐频。对于n次谐频，振幅为(8.1.3)初相为A=Ja2+b2nnn(8.1.4)(8.1.5a)-b申=arctannnanA和y称为振动的频谱，实际应用中只讨论振幅的频谱。nn方波的系数分别为(8.1.5b)2T/2a=JAdt=A0T0(8.1.6a)2T/22A1a=JAcosnetdt=sinnetnTTne0T/2=0(8.1.6b)2T/22A-1b=JAsinnetdt=cosnetnTTne0其中偶数项为零，因此系数为T/2A=[1-(-1)n]nn02Ab=2n-1(2n一1)n(8.1.6c)方波的展开式为/、A2Av®sin[(2n一1対t]x(t)=+乙2n2n-1n=1方波的频谱就是b。n[算法]取周期T为时间单位，取A为方波的单位，则得(8.1.7)2申sin[(2n一1)2nt*]x*=+—乙n2n一1n=1其中t*=t/T，x*=x/A。[程序]zqy8_1_1fourier.m如下。%方波的傅立叶分析(用解析式)clear%清除变量a0=0.5;tm=2;dt=0.001;t=0:dt:tm;%常数项%周期个数也是最大时间(1)%时间增量%时间向量(周期的倍数)%t0=[0,0,0.5,0.5,1];x0=[0,1,1,0,0];figureplot([t0,t0+1],[x0,x0],'LineWidth',2)gridonaxis([0,tm,-0.65,1.2])fs=16;xlabel('\itt\rm/\itT','FontSize',fs)ylabel('\itx\rm/\itA','FontSize',fs)%矩形波自变量%矩形波函数%创建图形窗口%画矩形波(2)%加网格%曲线范围%字体大小%x标签%y标签%holdonplot([0,tm],[a0,a0])%保持图像%画水平线(4)title('方波的傅立叶分析','FontSize',fs)%标题text(0,a0,num2str(a0),'FontSize',fs)n=1;x=a0;An=[];while1%显示平均值文本%求和下标%常量%频谱变量置空%无限循环h=plot(t,zeros(size(t)),'r','EraseMode','Xor');%画线并取句柄(5)%bn=2/pi/(2*n-1);An=[Anbn];xn=bn*sin((2*n-1)*2*pi*t);%%求正弦项的系数%连接频谱%求分量x=x+xn;set(h,'YData',x)%累加分量(9)%设合成线(10)plot(t,xn,'Color',rand(3,1))%画分支线yn=input('还要吗(y/n)?','s');%键盘输入iflower(yn)=='n',break,end%判断退出(11)n=n+1;%波的个数加1end%结束循环An(An<=5e-5)=nan;%振幅太小则取为非数(12)figure%创建图形窗口stem(1:2:2*n-1,An)%画杆图gridon%加网格xlabel('\it\omega','FontSize',fs)%x标签ylabel('\itA_n/A','FontSize',fs)%y标签title('方波的频谱图','FontSize',fs)%标题text(1:2:2*n-1,An,num2str(An',4))%显示频谱程序执行后首先画出脉冲方波和一个基频波以及基频波与平均值的和的曲线；在命令窗口不断回车，就不断画出谐频分波和分波的叠加曲线。最后按n回车，结束叠加过程。[图示]⑴如P8.1.1a图所示，当谐波很多时，它们的叠加曲线就越接近脉冲波。频谱图如P8.1.1b图所示，频率越大，频谱值越小。当频谱值比较大时，对合成波的影响也比较大。P8.1.1a图P8.1.1aP8.1.1a图方法二：用符号积分。如B8.1b图所示，单位阶梯函数又称为海维赛得函数，其定义为[0,(t<a)u(t—a)=](8.1.8)[l,(a<t)利用单位阶梯函数，方波可表示为x(t)=u(t)-u(t—).5)(8.1.9)利用欧拉公式cosx=(eix+e-ix)/2，sinx=(eix-e-ix)/2i，傅立叶级数可化为x(t)_a十牙exp(in®t)+exp(-in®t)十匕exp(in®t)-exp(-inwt).

兀⑴_~2+"n2+n2i)n_1a_—0■-22nn+a_—0■-22nn+2区(a-bi)exp(血)+2区(a+bi)exp(-i讪)nnn_1n_1利用系数关系a=—fx(t)cos(-net)dt=a-nTn0b=—fx(t)sin(-net)dt=-b-nTtOaB8.1b图可得x(利用系数关系a=—fx(t)cos(-net)dt=a-nTn0b=—fx(t)sin(-net)dt=-b-nTtOaB8.1b图可得x(t)=—艺(a2n=—g-bi)exp(inet)=艺cexp(in®t)nnn=—g(8.1.10)其中复系数为c=—(a-bi)=n2nnT01fx(t)(cosnet-isinn®t]dtc=—fx(t)exp(-inet)dtnT0(8.1.11)—、,:c*-n-n2nnnc=—(a-bi)=—(a+bi)=-n2c*表示c的复共轭。由此可知：nna0=2c0

a=2Re(c)

nn

b=-2Im(c)

nn(8.1.12a)(8.1.12b)(8.1.12c)n次谐波的振幅为A=\：a2+b2=21cInnn(8.1.13)设x=cnnexp(inat)，贝I」=[Tfx(T)exp(-ine)dT]exp(inet)=命fx(T)exp[ine(t—t)]dT00x-n0-nnn［算法］取周期T为时间单位，取A为方波的单位，复系数可表示为c=fx*(t*)exp(-i2nt*)dt*n0其中x*(t*)屮0<t*<1/2)(8.—.——*)(8.1.1*)10,(1/2<t*<1)

(8.—.——*)(8.1.1*)傅立叶级数可表示为x*(t*)=£cexp(in2n*)(8.1.10*)nn=—g利用海维赛德指令heaviside指令形成被积分函数，利用符号积分指令int求积分。替换整数即可形成系数，从而求得各项之值。［程序］zqy8_l_2fourier.m如下。%方波的傅立叶分析(用符号积分)clear%清除变量70symsstsn%定义符号变量sx=heaviside(st)-heaviside(st-0.5);a0=int(sx,st,0,1);%用单位阶梯函数形成矩形波%求常数项a0=double(a0);%常数项化为数值cn=int(sx*exp(-i*2*pi*sn*st),st,0,1);%求复数系数tm=2;%周期个数也是最大时间dt=0.001;%时间增量t=0:dt:tm;%t0=[0,0,0.5,0.5,1];%时间向量(周期的倍数)%矩形波自变量x0=[0,1,1,0,0];%矩形波函数figure%创建图形窗口plot([t0,t0+1],[x0,x0],'LineWidth',2)%画矩形波gridon%加网格axis([0,tm,-0.65,1.2])%曲线范围fs=16;%字体大小xlabel('\itt\rm/\itT','FontSize',fs)%x标签ylabel('\itx\rm/\itA','FontSize',fs)%y标签title('方波的傅立叶分析','FontSize',fs)%标题holdon%保持图像plot([0,tm],[a0,a0])%画水平线(4)h=plot(t,zeros(size(t))，‘r','EraseMode','Xor');%画线并取句柄(5)text(0,a0,num2str(a0),'FontSize',fs)%显示平均值文本n=1;%求和下标x=a0;%常量An=[];%频谱变量置空while1%c=double(subs(cn,sn,n));%无限循环%虚指数项系数化为数值An=[An,abs(2*c)];%连接频谱xn=2*real(c*exp(i*2*pi*n*t));%x=x+xn;%求分量%累加分量set(h,'YData',x)%设合成线

plot(t,xn,'Color',rand(3,1))%画分支线yn=input('还要吗(y/n)?','s');%键盘输入iflower(yn)=='n',break,end%判断退出n=n+1;%波的个数加1end%结束循环An(An<1e-5)=nan;%太小的数改为非数figure%创建图形窗口stem(1:n,An)%画杆图gridon%加网格xlabel('\it\omega','FontSize',fs)%x标签ylabel('\itA_n/A','FontSize',fs)%y标签title('方波的频谱图','FontSize',fs)%标题text(1:n,An,num2str(An',4))%显示频谱方法三：用快速傅立叶变换指令fft。假设有N个离散的数列xk，(k=1,2,…，N),可按傅立叶公式展开TOC\o"1-5"\h\z1令n—1x二工Cexp[i2n(k—1)](k=1，2,...,N)(8.1.14)kNnNn=1其中Nn—1C=》xexp[—i2n(k—1)](n=1，2,N)(8.1.15)nkNk=1[算法]对于离散的数列，可用快速傅立叶变换指令fft求系数Co当N是2的整数幕的n时候，fft的效率最高。在傅立叶级数中,系数关系为IC/N(n=0,1,...,N/2)(8.1.15*)C=2n+1(8.1.15*)n[C/N(n=—1,...,—N/2+1)N+n+1这里，C1/N等于平均值c0,C2/N、C』N、…、Cn/2+1/N等于&c2、…、cN/2,而Cn/2+2/N、…、CnJN、CIN等于…、c2、j。由于c=c*,在变换之后只要取C的前半部数据N1NN/2+121nnn就行了。[程序]zqy8_l_3fourier.m如下。%方波的傅立叶分析(用快速傅立叶变换fft指令)clear%清除变量n=212;t=linspace(0,1,n);x=zeros(size(t));x(t<=0.5)=1;Cn=fft(x);cn=Cn/n;n=212;t=linspace(0,1,n);x=zeros(size(t));x(t<=0.5)=1;Cn=fft(x);cn=Cn/n;a0=cn(1);tm=2;%时间向量(周期的倍数)%函数值取全零向量%函数符合条件的部分取1%离散序列的傅立叶变换(2)%求复系数%平均值%周期个数也是最大时间

70t0=[0,0,0.5,0.5,1];%矩形波自变量x0=[0,1,1,0,0];%矩形波函数figure%创建图形窗口plot([t0,t0+1],[x0,x0],'LineWidth',2)%画矩形波gridon%加网格axis([0,tm,-0.65,1.2])%曲线范围fs=16;%字体大小xlabel('\itt\rm/\itT','FontSize',fs)%x标签ylabel('\itx\rm/\itA','FontSize',fs)%y标签title('方波的傅立叶分析','FontSize',fs)%标题holdon%保持图像plot([0,tm],[a0,a0])%画水平线h=plot(t,zeros(size(t)),'r','EraseMode,'Xor');%画线并取句柄text(0,a0,num2str(a0),'FontSize',fs)%显示平均值文本n=1;%求和下标x=a0;%常量An=[];%频谱变量置空while1QZ-%无限循环%an=2*real(cn(n+1));%取实部(4)bn=-2*imag(cn(n+1));%取虚部(4)An=[An,2*abs(cn(n+1))];%连接频谱(4)xn=an*cos(n*2*pi*t)+bn*sin(n*2*pi*t);%求分量(4)/Ox=x+xn;%累加分量set(h,'YData',x)%设合成线plot(t,xn,'Color',rand(3,1))%画分支线yn=input('还要吗(y/n)?','s');%键盘输入iflower(yn)=='n',break,end%判断退出n=n+1;%波的个数加1end%结束循环figure%创建图形窗口stem(1:n,An)%画杆图gridon%加网格xlabel('\it\omega','FontSize',fs)%x标签ylab