◆例1设方波的数学模型为：]5sin 513sin 31[sin 4)(000 t t t E t f T ，基频： T20用MATLAB 软件完成该方波的合成设计◆ MATLAB 源程序t=-10:0.1:10; %设定一个数组有201个点，方波周期为20e=5;w=pi/10; %设定方波幅值为5，w 代表w0m=-5*sign(t); %给定幅值为5的方波函数y1=(-4*e/pi)*sin(w*t); %计算1次谐波y3=(-4*e/pi)*(sin(w*t)+sin(3*w*t)/3); %计算3次谐波y5=(-4*e/pi)*(sin(w*t)+sin(3*w*t)/3+sin(5*w*t)/5); %计算5次谐波plot(t,y1,'y');hold; grid; %用黄色点线画出1次谐波及网格线,并在同一张图上画其余曲线plot(t,y3,'g'); %用绿色点线画出3次谐波plot(t,y5,'b'); %用蓝色点线画出5次谐波plot(t,m,'-k'); %用黑色实线画方波title('方波合成');xlabel('t');ylabel('f(t)'); %为图形加上标题n=50; %合成任意次方波，n 决定方波的合成次数，在此给定50yn=0; %设置初始值for i=1:nyn=yn+(-4*e/pi)*(1/(2*i-1))*sin((2*i-1)*w*t);end; %计算n 次谐波合成plot(t,yn,'r') %用红色实线画出n 次谐波合成◆ 从图中我们可以看到Gibbs 现象。

在函数的间断点附近，增加傅里叶级数的展开次数，虽然可以使其间断点附近的微小振动的周期变小，但振幅却不能变小。

此现象在控制系统表现为：当求控制系统对阶跃函数的响应时，超调量总是存在的。

◆例2（P110）MATLAB中函数FFT应用举例。

%MATLAB中函数FFT应用举例t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(size(t))subplot(2,1,1)plot(y(1:50))xlable(‘时间轴t’)ylable(‘信号值f(t)’)title(‘正弦波＋随即噪声’,’FontSize’,10)y=fft(y,512);f=1000*(0:256)/512;subplot(2,1,2)plot(f,Y(1:257))set(gca,’Xtick’,[0,50,100,150,200,250,300,350,400,450,500])set(gca,’XtickLabel’,’0|50|100|150|200|250|300|350|400|450|500|’) xlabel(‘频率轴\omega’)ylabel(‘频谱幅值F(\omega)’)title(‘信号频谱’,’FontSize’,10)◆例3例 3.8.3 有一二阶系统，阻尼比 ＝0.47，固有频率500 n Hz ，采样间隔0004.0 t s,采样点数N ＝256。

试计算理论幅频特性与由系统阶跃响应计算出的幅频特性数据值，并画出两个计算结果的幅频特性曲线。

％example 3.8.3 MA TLAB PROGRAMN=256;dt=0.0004wn=500;seta=0.47;dw=2*pi/(N*dt);a=wn^2;b=[1,2*seta*wn,a];t=[0:dt:(N-1)*dt];c=step(a,b,t);w=[0:dw:(N-1)*dw];[mag,phase]=bode(a,b,w);ycw=fft(c);Re=real(ycw);Im=imag(ycw);for i=1:NRw(i,1)=1-Im(i,1)*(i-1)*dw*dt;Iw(i,1)=Re(i,1)*(i-1)*dw*dt;endffw=Rw+Iw*sqrt(-1);Aw=abs(ffw)semilogx(w,20*log10(mag),'r-')axis([100,10000,-30,10])text(600,12,'幅频特性')hold onsemilogx(w,20*log10(Aw))axis([100,10000,-30,10])grid on◆例4例6.2.4 用MATLAB 中的函数XCORR 求出下列两个周期信号的互相关函数，式中的f=10Hz 。

),2sin()(ft t x )902sin(5.0)( ft t y％例6.2.4 中计算两个周期信号互相关函数的MA TLAB 程序N=500;Fs=500;n=0:N-1;t=n/Fs;Lag=200;x=sin(2*pi*10*t);y=0.5*sin(2*pi*10*t+90*pi/180);[c,lags]=xcorr(x,y,Lag,'unbiased');subplot(2,1,1)plot(t,x,'r')hold onplot(t,y,'b')xlabel('t');ylabel('x(t)y(t)');title('原周期信号')gridhold offsubplot(2,1,2);plot(lags/Fs,c,'k');xlabel('t');ylabel('Rxy(t)');title('互相关函数');grid◆例5例 6.2.5 若有信号为)()2sin(2)2sin()(21t t f t f t x式中，,501Hz f Hz f 1002 ，)(t 为白噪声（用MATLAB 中的函数产生）。

设采样频率2000 Fs ；试用周期法并应用MA TLAB 编程计算，当数据长度分别为2561 N 和 10242 N 两种情况下上述信号的功率谱。

%例6.2.5中周期图法计算信号功率谱的MATLAB 程序clfFs=2000;%情况1:数据长度N1=256N1=256;N1fft=256;n1=0:N1-1;t1=n1/Fs;f1=50;f2=100;xn1=sin(2*pi*f1*t1)+sin(2*pi*f2*t1)+randn(1,N1);Pxx1=10*log10(abs(fft(xn1,N1fft).^2)/N1);f1=(0:length(Pxx1)-1)*Fs/length(Pxx1);subplot(2,1,1)plot(f1,Pxx1)ylabel('功率谱(dB)');title('数据长度N1=256')grid%情况2:数据长度N2=1024N2=1024;N2fft=1024;n2=0:N2-1;t2=n2/Fs;f1=50;f2=100;xn2=sin(2*pi*f1*t2)+2*sin(2*pi*f2*t2)+randn(1,N2);Pxx2=10*log10(abs(fft(xn2,N2fft).^2)/N2);f2=(0:length(Pxx2)-1)*Fs/length(Pxx2);subplot(2,1,2)plot(f2,Pxx2)xlabel('频率(Hz)');ylabel('功率谱(dB)');title('数据长度N2=1024')◆例6（例3.8.1）分别用conv和FFT算法计算序列：x(n)为在区间[0，1]上均匀分布的N 点随机序列，表示为：x(n)＝rand(1,N),h(n)为均值为零、方差为1的N点高斯分布随机序列，表示为：h(n)＝rand(1,L),试求1≤N≤150时的卷积并比较其运算时间。

%例3.8.1 直接卷积和快速卷积的比较%conv_time=zeros(1,150);fft_time=zeros(1,150);%for N=1:150tc=0;tf=0; %初始化L=2*N-1; %加长序列长度nu=round((log10(L)/log10(2))+0.45);L=2^nu; %使点数成为2的幂次for I=1:100h=randn(1,N); %产生两个随机序列x=rand(1,N);t0=clock;yc=conv(h,x); %直接卷积计算t1=etime(clock,t0);tc=tc+t1; %直接卷积运算的时间t0=clock;y2=ifft(fft(h,L).*fft(x,L)); %快速卷积计算t2=etime(clock,t0);tf=tf+t2; %快速卷积计算的时间end%conv_time(N)=tc/100; %直接卷积计算的平均时间fft_time(N)=tf/100; %快速卷积计算的平均时间end%n=1:150;subplot(1,1,1); %图形显示上述两种卷积的计算时间plot(n(25:150),conv_time(25:150),n(25:150),fft_time(25:150))上述两种卷积的计算时间的比较如图3.8.5所示。

图3.8.5 两种卷积计算时间比较◆例7（例3.8.2 ）运用FFT 求取矩形脉冲 1,01t 0 ,1)(t t 的谱，说明采样频率低引起的混叠现象。

（1） 先编写有一定通用性的函数文件cftbyfft.m[cftbyfft.m]function [AW,f]=cftbyfft(wt,t,flag)%cftbyfft.m%本程序采用FFT 计算连续时间Fouie 变换。

输出幅频数据对(f,AW)。

%输入量(wt,t)为已经窗口化了的时间函数wt(t)，它们分别是长度为N 的向量。

%对于时限信号，应使该取值时段与窗口长度相比足够小。

以提高频率分辨率。

%对于非时限信号，窗口长度的选取应使窗口外的函数值小%到可忽略，以提高近似精度。

%输入宗量flag 控制输出CFT 的频率范围。

% flag 取非0时(缺省使用)，频率范围在[0,fs);% flag 取0时，频率范围在[-fs/2,fs/2)。

if nargin==2;flag=1;endN=length(t); %采样点数，应为2的幂次，以求快速。

T=t(length(t))-t(1); %窗口长度dt=T/N; %时间分辨率W0=fft(wt); %实施FFT 变换W=dt*W0; %算得[0,fs]上的N 点CFT 值df=1/T; %频率分辨率n=0:1:(N-1)%把以上计算结果改写到[-fs/2,fs/2]范围if flag==0n=-N/2:(N/2-1);W=fftshift(W); %产生周期序列的频谱endf=n*df; %频率分量向量AW=abs(W); %幅频谱数据向量if nargout==0plot(f,AW);grid,xlabel('频率f');ylabel('|w(f)|')end（2）运行以下指令，绘制时域波形和幅频谱M=5; %做2的幂次用。