1、1例 1设方波的数学模型为： ，基频： 5sin13sisin4)( 000 tttEtfT 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)/
2、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:n yn=yn+(-4*e/pi)*(1/(2*i-1)*sin(2*i-1)*w*t);end; %计算 n次谐波合
3、成plot(t,yn,r) %用红色实线画出 n次谐波合成 从图中我们可以看到 Gibbs 现象。在函数的间断点附近，增加傅里叶级数的展开次数，虽然可以使其间断点附近的微小振动的周期变小，但振幅却不能变小。此现象在控制系统表现为：当求控制系统对阶跃函数的响应时，超调量总是存在的。23例 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)ylab
4、le(信号值 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
5、，固有频率 Hz，采样间隔50ns,采样点数 N256。试计算理论幅频特性与由系统阶跃响应计算出的幅频特04.t性数据值，并画出两个计算结果的幅频特性曲线。example 3.8.3 MATLAB PROGRAMN=256;dt=0.0004wn=500;seta=0.47;dw=2*pi/(N*dt);a=wn2;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-I
6、m(i,1)*(i-1)*dw*dt;Iw(i,1)=Re(i,1)*(i-1)*dw*dt;end4ffw=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()fttx)902sin(5.0(ftty5例 6.2
7、.4 中计算两个周期信号互相关函数的 MATLAB 程序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(互
8、相关函数);grid例 5例 6.2.5 若有信号为)(2sin()2sin()1tftftx式中， ， 为白噪声（用 MATLAB 中的函数产生） 。设采样,501Hzfz)t频率 ；试用周期法并应用 MATLAB 编程计算，当数据长度分别为 和Fs 2561N两种情况下上述信号的功率谱。24N%例 6.2.5 中周期图法计算信号功率谱的 MATLAB 程序clfFs=2000;%情况 1:数据长度 N1=256N1=256;N1fft=256;n1=0:N1-1;6t1=n1/Fs;f1=50;f2=100;xn1=sin(2*pi*f1*t1)+sin(2*pi*f2*t1)+randn
9、(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)
10、.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), 试求 1N150 时的卷积并比较其运算时间。%例 3.8.1 直接卷积和快速卷积的比较%conv_ti
11、me=zeros(1,150);fft_time=zeros(1,150);%for N=1:150 tc=0;tf=0; %初始化 L=2*N-1; %加长序列长度nu=round(log10(L)/log10(2)+0.45);L=2nu; %使点数成为 2 的幂次for I=1:100h=randn(1,N); %产生两个随机序列x=rand(1,N);7t0=clock;yc=conv(h,x); %直接卷积计算t1=etime(clock,t0);tc=tc+t1; %直接卷积运算的时间t0=clock;y2=ifft(fft(h,L).*fft(x,L); %快速卷积计算t2=et
12、ime(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 ,0t)(t引起的混叠
13、现象。（1） 先编写有一定通用性的函数文件 cftbyfft.mcftbyfft.mfunction AW,f=cftbyfft(wt,t,flag)8%cftbyfft.m%本程序采用 FFT 计算连续时间 Fouie 变换。输出幅频数据对(f,AW)。%输入量(wt,t)为已经窗口化了的时间函数 wt(t)，它们分别是长度为 N 的向量。%对于时限信号，应使该取值时段与窗口长度相比足够小。以提高频率分辨率。%对于非时限信号，窗口长度的选取应使窗口外的函数值小%到可忽略，以提高近似精度。%输入宗量 flag 控制输出 CFT 的频率范围。% flag 取非 0 时(缺省使用)，频率范围在0,
14、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(
15、W); %幅频谱数据向量if nargout=0plot(f,AW);grid,xlabel(频率 f);ylabel(|w(f)|)end（2） 运行以下指令，绘制时域波形和幅频谱M=5; %做 2 的幂次用。本例把 M 设得较小，是为了观察混叠。tend=1; %波形取非零值的时间长度。T=10; %窗口化长度应足够大，以减小窗口化引起的泄露“旁瓣”效应。N=2M; %采样点数，取 2 的幂是为了使 FFT 运算较快。dt=T/N; %以上 T、N 的取值应使 N/Tfs 采样频率大于两倍时间波形%带宽，以克服采样引起的频谱混叠。%在本例中，据理论分析知 W(f=7.5)=Sa(7.5*p
16、i)=1/(7.5%*pi)0); %产生非零波形时段的相应序列w(Tow,1)=ones(length(Tow),1); %在窗口时段内定义的完整波形plot(t,w,b,LineWidth,2.5),title(Time Waveform);xlabel(t-)由上述程序画出的被窗口化的时间波形如图 3.8.12 所示：9图 3.8.12 被窗口化的时间波形AW,f=cftbyfft(w,t,0);ff=f+eps; %为避免下面指令出现 0/0 而采取的措施AWW=abs(sin(pi*ff)./(pi*ff);plot(f,AW,b-,ff,AWW,r:);title(Aliasing
17、 caused by undersampling)xlabel(f);ylabel(|W(f)|),legend(by FFT,Theoretical)由上述程序画出的“欠”采样时引起的混叠曲线如图 3.8.13 所示：图 3.8.13 “欠”采样时引起的混叠10例 8（例 3.8.5） 试分别求出含白噪声干扰的正弦信号与白噪声的自相关函数，并对它们的结果进行比较。%example 3.8.5 MATLAB PROGRAMN=1000;n=0:N-1;Fs=500;t=n/Fs;Lag=100; %相关信号的最大延迟量x1=sin(2*pi*10*t)+0.6*randn(1,length(t
18、); %含白噪声的正弦信号 x1c,lags=xcorr(x1,Lag,unbiased) %无偏自相关函数的计算subplot(2,2,1) %画出 x1 曲线plot(t,x1)xlabel(t);ylabel(x1(t);title(含白噪声的正弦信号 x1);gridsubplot(2,2,2) %画出 x1 自相关曲线plot(lags/Fs,c)xlabel(t);ylabel(Rxx1(t);title(x1 的自相关函数);gridx2=randn(1,length(t); %发生白噪声 x2c,lags=xcorr(x2,Lag,unbiased) %白噪声 x2 的无偏自相
19、关函数subplot(2,2,3) %画出 x2 曲线plot(t,x2)xlabel(t);ylabel(x2(t);title(白噪声 x2);gridsubplot(2,2,4) %画出 x2 自相关曲线plot(lags/Fs,c)xlabel(t);ylabel(Rxx2(t);title(x2 的自相关函数);grid由上述程序画出的信号与自相关函数的曲线如图 3.8.19 所示：11图 3.8.19 信号与自相关函数由计算结果和和图 3.8.19 的曲线可以看出：同时具有周期性和白噪声干扰的信号，其自相关函数不但在 0 时具有最大值，而且在 较大时自相关函数仍有明显的周期性，它的
20、频率和原周期信号的周期相同；而无周期的白噪声，自相关函数在 0 时也具有最大值，但当 稍大时迅速衰减至零附近，利用自相关函数的这一特性可用来识别随机信号中是否含有周期成分和它的频率。 例 9（例 4.4.3） 试用 MATLAB 语言绘制巴特沃思低通模拟滤波器的平方幅频相应曲线，滤波器的阶数分别为 2，5，10，20。解例 4.4.3 的 MATLAB 程序n=0:0.01:2for i=1:4switch icase 1N=2;case 2N=5;case 3N=10; case 4N=20;end z,p,k=buttap(N); %巴特沃思滤波器原型设计函数，n:阶数;z,p,k:12%
21、滤波器零点，极点和增益b,a=zp2tf(z,p,k); % 零极点增益转换为传递函数H,w=freqs(b,a,n); %模拟滤波器频率响应函数magH2=(abs(H).2; %幅度平方函数hold onplot(w,magH2);axis(0201);endxlabel(w/wc) ;ylabel(H(jw)2 );grid图 4.4.4 不同阶次巴特沃思滤波器的幅度平方函数 例 10（例 4.4.6 ）试用 MATLAB 程序，确定一个模拟低通滤波器的阶数 N 和截止频率 Wc。设计指标为：通带边界频率 200，阻带边界频率 300，通带波纹为 1dB，在 Z 处，幅度衰减大于 18
22、dB。解 设计切比雪夫低通滤波器的 MATLAB 程序wp=200*pi;ws=300*pi;Rp=1; %通带波纹Rs=16; %阻带衰减%计算滤波器阶数 ebs=sqrt(10(Rp/10)-1); %波纹系数 A=10(Rs/20); %A 为参变量 Wc=wp; %截止频率 Wr=ws/wp; g=sqrt(A*A-1)/ebs; %g 为参变量N1=log10(g+sqrt(g*g-1)/log10(Wr+sqrt(Wr*Wr-1);%滤波器阶数计算N=ceil(N1); %N 取整数运行上述程序后，可得滤波器的截止频率 Wc 和阶数 N 为Wc=628.3185N413 例 11（
23、例 4.4.4） 绘制出阶数分别为 2，4，6，8 的切比雪夫模拟低通滤波器的平方幅频响应曲线。解 写出其 MATLAB 程序如下：绘制切比雪夫低通滤波器幅度平方函数的 MATLAB 程序n=0：0.01：2；for i=1:4switch icase 1N=2；case 2N=4；case 3N=6；case 4N=8；endRp=1; 滤波器通带波纹系数z,p,k=cheblap(N,Rp); 设计切比雪夫模拟低通滤波器原型函数，z,p,k 分别为滤波器的零点、极点和增益b,a=zp2tf(z,p,k); 零点、极点和增益转换为传递函数H,w=freqs(b,a,n); 模拟滤波器的频率响
24、应magH2=(abs(H).2; 幅度平方函数Outputposplot=10num2str(i); 定义字符串变量subplot(posplot)plot(w,magH2,k);ylabel(|H(jw)2|);title(N=num2str(N);gridend其幅度平方函数曲线如图 4.4.6 所示。由图 4.4.6 可以看出：切比雪夫滤波器在通带内具有等波纹起伏特性，而在阻带内则单调下降。阶次越高，特性越接近矩形。 例 12（例 4.4.7） 试设计一巴特沃思模拟带通滤波器，设计要求为：通带频率2kHz3kHz,两边的过渡带宽为 0.5kHz，若通带波纹 1dB，阻带衰减大于 100
25、dB.解 设计巴特沃思带通滤波器 MATLAB 程序滤波器设计指标wp=2000 3000*2*pi ；ws=1500 3500*2*pi ；Rp=1;14Rs=100;%计算阶数和截止频率N,Wn=buttord (wp,ws,Rp,Rs,s);Fc=Wn/(2*pi);计算滤波器传递函数多项式系数b,a=butter(N,Wn,s);%画出滤波器幅频特性w=linspace(1,4000,1000)*2*pi; %生成线性等间隔向量H=freqs(b,a,w);magH=abs(H);phaH=unwrap(angle(H);plot(w/(2*pi),20*log10(magH),k);
26、xlabel(频率（Hz）)；ylabel(幅度（dB）);title(巴特沃思模拟滤波器 )grid on由上述程序画出的巴特沃思模拟滤波器的幅频特性如图 4.4.11 所示。图 4.4.11 巴特沃思带通滤波器的幅频特性 例 13（例 4.4.8 ） 试设计一个切比雪夫模拟带阻滤波器。要求的指标为：阻带上下边界频率为 2KHz 与 3KHz，两侧过渡带为 0.5KHz，通带波纹 1dB,阻带衰减大于 60 dB。解设计切比雪夫带阻滤波器 MATLAB 程序滤波器设计指标wp=2000 3000*2*pi ；ws=1500 3500*2*pi ；Rp=1;Rs=60;%计算阶数和截止频率N,
27、Wn=cheblord (wp,ws,Rp,Rs,s);计算滤波器传递函数多项式系数b,a=chebyl(N,Rp,Wn,stop,s);%画出滤波器幅频特性w=linspace(1,4000,1000)*2*pi; %生成线性等间隔向量H=freqs(b,a,w); %相位展开magH=abs(H);phaH=unwrap(angle(H);plot(w/(2*pi),20*log10(magH),k);xlabel(频率（Hz）)；ylabel(幅度（dB）);15title(切比雪夫模拟滤波器 )grid on由上述程序画出的切比雪夫模拟滤波器的幅频特性如图 4.4.12 所示。图 4.
28、4.12 切比雪夫带通滤波器的幅频特性 例 14（例 5.3.4） 一低通滤波器，其通带和阻带的技术指标分别是dB 25.0R ,2.0pp,3.zz试用频率抽样法设计一 FIR 具有线性相位的滤波器，取 N=20 。解 由 19,0k ,)()()(202 jdkNjdeHekH当 k=2 时，正好是在通带边界频率 处有一个频率抽样点，即 p20.下一个抽样点为 k=3 ，是阻带上边界频率 ，阻带域通带间无过渡带，即z3.0z则在通带 内有三个抽样点，阻带 上共有七个抽样点。从而有p0 p 1,0 ,1 ,kH由 N=20 ，相位常数 其相位可表示为5.92/)(/)(N1k0 , )20(
29、95kk根据式（5.3.37）得到 H(k) ，利用离散傅里叶变换，求得 h(n) ，并有频响内插公式最后可得 FIR 滤波器的频响 。其 MATLAB 程序如下：)(jeH%例 5.3.4 中设计 FIR 滤波器 MATLAB 程序N=20； ; 1*N)/pi2( wl;1N: 0;2/)1(alph； 1:)2/1flor(2k ;)/1flor( :0k ; ,5. ,wdl ,Hd5zes rs 16angH = -alpha * ( 2 * pi ) / N * kl ,alpha * (2 * pi ) / N * ( N-k2 ) ;H = Hrs. * exp (j * an
30、gH) ;= real (ifft (H , N ) ; db , mag , pha ,grd , w = freqz m ( h , 1 ) ;subplot (1 , 1 , 1 )subplot ( 2 , 2 , 1 ) ; plot (wl ( 1 : 11 ) , o , wdl , Hdr ) ;axis ( 0 , 1 , - 0.1 , 1.1 ) ; title ( 频率样本： N20 )xlabel (频率（单位： pi） ) ;ylabel (Hr (k)set (gca ,XtickMode,manual,XTick,0 , 0.2 , 0.3 , 1 );set
31、(gca ,YtickMode,manual,Ytick,0 , 1); girdsubplot(2 , 2 , 2) ; stem(l , h) ; axis(-1 , N ,-0.1 , 0.3)title (单位抽样响应 ) ; ylabel(h(n) ; text (N+1 , -0.1 , n)subplot (2 , 2 , 3 ) ; plot (ww/pi , Hr , wl (1:11)/pi , Hrs (1:11)/pi , Hrs(1:11) , o);axis ( 0 , 1 , -0.2 , 1.2 ) ; title (振幅响应)xlabel (频率（单位 : p
32、i ） ) ; ylable (Hr(w)set (gca , XTickMode , Xtick , 0 , 0.2 , 0.3 , 1)set (gca ,YtickMode , manual , Ytick , 0 , 1) ; girdsubplot (2, 2, 4 ) ; plot (w / pi , db) ; axis (0 , 1 200 , 10) ; girdtitle (幅度响应) ; xlabel ( 频率（单位：pi） ) ;ylable (分贝数 ) ;set (gca , XtickMode,Manual, XTick, 0 ; 0.2 ; 0.3 ; 1) ;
33、set (gca , YtickMode,Manual,YTick, -16 ; 0) ;set (gca , YtickLableMode,Manual,YTickLabels, 16 ; 0) ;上面程序中的 fteqz_m 和 Hr_Type2 为扩展函数，不在 MATLAB 所带的工具箱内。需要另行编制，函数 fteqz_m 为function db , mag , pha , grd , w = freqz_m (b , a);%db , mag , pha , grd , w = freqz_m (b , a);% db : 0- (dB)% mag : 0-% pha : 0-%
34、 grd : 0-% w : 0-% b : % a : %H , w = freqz (b , a , 1000 , whole) ;H = ( H (1 : 1 : 501 ) ) ;ag = abs ( H ) ;db = 20 * log 10 ( ( mag + eps ) / max ( mag ) ) ；pha = angle ( H ) ;grd = grpdelay (b , a , w ) ;17扩展函数 Hr_Type2 的 MATLAB 程序为% Hr , w , b , L = Hr_Type2 (h ) % Hr : % w : 0-% b :L : Hr% h :%M = length ( h ) ;L = M / 2 ;= 2 * h (L :-1 : 1) ;n = 1 : 1 : L ; n = n 0.5 ;w = 0 : 1 : 500* pi / 500 ;H = cos ( w * n )* b ;相应的单位抽样响应 h (n)，频响特性 （分别用相对幅度和 dB 表示）如图)(jeH5.3.13 所示。参考文献：周敏浩编著.信号处理技术基础.北京航空航天大学出版社，2001.
