MATLAB信号处理例题

1、共享知识分享快乐卑微如蝼蚁、坚强似大象例1设方波的数学模型为:,基频:4E11f (t)sin 0t Sin3 0tSin5 0t35用MATLA软件完成该方波的合成设计 MATLAB源 程序t=-10:0.1:10; % 设定一个数组有201个点，方波周期为20e=5;w=pi/10; % 设定方波幅值为5, W代表w0m=-5*sign(t); % 给定幅值为5的方波函数y1=(-4*epi)*si n(w*t); % 计算 1 次谐波y3=(-4*epi)*(sin(w*t)+sin(3*w*t)3); %计算 3 次谐波y5=(-4*epi)*(sin(w*t)+sin(3*w*t)3

2、+sin(5*w*t)5); %计算 5 次谐波plot(t,y1,'y');hold; grid; %用黄色点线画出1次谐波及网格线,并在同一张图上画其余曲线plot(t,y3,'g'); % plot(t,y5,'b'); % plot(t,m,'-k'); %用绿色点线画出3次谐波 用蓝色点线画出5次谐波 用黑色实线画方波title(' 方波合成');XIabeI('t')ylabel('f(t)'); %为图形加上标题n=50; %合成任意次方波，n决定方波的合成次数，在此给

3、定50yn=0; %设置初始值for i=1: nyn=yn+(-4*epi)*(1(2*i-1)*si n(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=Si n( 2*pi*5

4、0*t)+si n(2*pi*120*t);y=x+2*ra ndn( size(t)SUbPlot(2,1,1)plot(y(1:50)xlable(时间轴 t')ylable(信号值 f(t)')title(芷弦波 + 随即噪声','Fo ntSize ',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,'XtickL

5、abel ',X50100150200250300350400450500') xlabel(频率轴 omega')ylabel(频谱幅值 F(omega)title(信号频谱',XontSize)0)例3例3.8.3 有一二阶系统，阻尼比 =0.47,固有频率 n 500 HZ ,采样间隔 t0.0004s采样点数N = 256。试计算理论幅频特性与由系统阶跃响应计算出的幅频特性数据值，并画出两个计算结果的幅频特性曲线。% example 3.8.3 MATLAB PRoGRAMN=256;dt=0.0004Wn=500;Seta=0.47;dw=2*pi(N

6、*dt);a=w n2;b=1,2*seta*w n,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(6

7、00,12,'幅频特性')hold Onsemilogx(w,20*log10(Aw)axis(100,10000,-30,10)grid On幅频特性例4例6.2.4用MATLAB中的函数XCoRR求出下列两个周期信号的互相关函数，式中的 f=10Hz。x(t) Sin(2 ft), y(t) 0.5sin(2 ft 90 )%例6.2.4中计算两个周期信号互相关函数的MATLAB程序N=500;Fs=500;n=0:N-1;t=n/Fs;Lag=200;X=S in (2*pi*10*t);y=0.5*si n(2*pi*10*t+90*pi180);c,lags=xcor

8、r(x,y丄ag,' Un biased');SUbPlot(2,1,1)plot(t,x,'r')hold Onplot(t,y,'b')XIabel('t')ylabel('x(t)y(t)');title('原周期信号')gridhold offsubplot(2,1,2);Plot(IagS/Fs,c,'k');XIabel('t')ylabel('Rxy(t)'); title('互相关函数');grid例5例6.2.5若有信

9、号为x(t) sin(2 fit) 2sin(2 f2t)(t)。设采样频256和式中，f150Hz, f2100Hz ,(t)为白噪声(用 MATLAB中的函数产生)率FS 2000;试用周期法并应用MATLAB编程计算，当数据长度分别为N1N21024两种情况下上述信号的功率谱。%例6.2.5中周期图法计算信号功率谱的MATLAB程序CIfFs=2000;%情况1：数据长度N1=256N1=256;NIfft=256;n1=0:N1-1;t1= n1Fs;f1=50;f2=100;Xn 1=si n(2*pi*f1*t1)+si n( 2*pi*f2*t1)+ra ndn (1,N1);P

10、XX仁 10*log10(abs(fft(x n1,N1fft).2)N1);f1=(0:Ie ngth(Pxx1)-1)*Fsle ngth(Pxx1);subplot(2,1,1)plot(f1,Pxx1)ylabel('功率谱(dB)');title('数据长度N仁256')grid%情况2:数据长度N2=1024N2=1024;N2fft=1024;n2=0:N2-1;t2=n2Fs;f1=50;f2=100;Xn 2=si n(2*pi*f1*t2)+2*si n( 2*pi*f2*t2)+ra ndn (1,N2);Pxx2=10*log10(abs

11、(fft(x n2,N2fft).2)N2); f2=(0:Ie ngth(Pxx2)-1)*Fsle ngth(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

12、N 150时的卷积并比较其运算时间。9例3.8.1直接卷积和快速卷积的比较%Con v_time=ZeroS(1,150);fft_time=ZeroS(1,150);%for N=1:150tc=0;tf=0;%初始化L=2*N-1;%加长序列长度n u=rou nd(log10(L)log10(2)+0.45);L=2nu;%使点数成为2的幕次for 1=1:100h=ra ndn (1,N);%产生两个随机序列x=ran d(1,N);t0=clock;yc=c On v(h,x);%直接卷积计算t仁 etime(clock,t0);tc=tc+t1;%直接卷积运算的时间t0=clock

13、;y2=ifft(fft(h,L).*fft(x,L);%快速卷积计算t2=etime(clock,t0);tf=tf+t2;%快速卷积计算的时间end%con v_time(N)=tc/100;%直接卷积计算的平均时间fft_time(N)=tf/100;%快速卷积计算的平均时间end%n=1:150;SUbPlot(1,1,1);%图形显示上述两种卷积的计算时间plot( n(25:150),co nv_time(25:150), n(25:150),fft_time(25:150)上述两种卷积的计算时间的比较如图3.8.5所示。图385 两种卷积计算时间比较例7 (例3.8.2 )运用F

14、FT求取矩形脉冲(t)1,0 t 10, t 1的谱，说明采样频率低引起的混叠现象。(1)先编写有一定通用性的函数文件Cftbyfft.mcftbyfft.mfun ctio n AW,f=cftbyfft(wt,t,flag)%cftbyfft.m%本程序采用FFT计算连续时间FoUie变换。输出幅频数据对(f,AW)。%俞入量(wt,t)为已经窗口化了的时间函数wt(t),它们分别是长度为N的向量。%寸于时限信号，应使该取值时段与窗口长度相比足够小。以提高频率分辨率。%寸于非时限信号，窗口长度的选取应使窗口外的函数值小%到可忽略，以提高近似精度。%俞入宗量flag控制输出CFT的频率范围。

15、% flag 取非0时(缺省使用)，频率范围在0,fs);% flag 取0时，频率范围在-fs2,fs2)采样点数，应为2的幕次，以求快速。 窗口长度时间分辨率实施FFT变换算得0,fs 上的N点CFT值频率分辨率if nargin=2;flag=1;e ndN=le ngth(t);%T=t(le ngth(t)-t(1);%dt=TN;%W0=fft(wt);%df=1T;%W=dt*W0;%n=0:1:(N-1)%e以上计算结果改写到if flag=On=-N/2:(N/2-1);W=fftshift(W);%end f=n *df;%AW=abs(W);%if n argout=0p

16、lot(f,AW);grid,xlabel('-fs2,fs2范围产生周期序列的频谱频率分量向量幅频谱数据向量频率 f)ylabel('w(f)')end(2)运行以下指令,绘制时域波形和幅频谱M=5;%做2的幕次用。本例把 M设得较小，是为了观察混叠。ten d=1;%波形取非零值的时间长度。T=10;%窗口化长度应足够大，以减小窗口化引起的泄露“旁瓣”效应。N=2M;%采样点数，取2的幕是为了使FFT运算较快。dt=TN;%以上T、N的取值应使N/T = fs采样频率大于两倍时间波形%带宽，以克服采样引起的频谱混叠。%在本例中，据理论分析知W(f=7.5)=Sa (

17、7.5*pi)=1(7.5%*pi)<5%.因此，可近似认为本例时间信号带宽为7.5H z.n=0:N-1;%采样序列t=n *dt;%采样点时间序列W=ZerOS(SiZe(t,2),1);Tow=fi nd(te nd-t)>O);%产生非零波形时段的相应序列w(Tow,1)=o nes(le ngth(Tow),1);%在窗口时段内定义的完整波形plot(t,w,'b','Li neWidth',2.5),title('Time WaVeform');XIabeI('t->')由上述程序画出的被窗口化的时间

18、波形如图3.8.12所示：图3.8.12被窗口化的时间波形AW,f=cftbyfft(w,t,O);ff=f+eps;%为避免下面指令出现0/0而采取的措施AWW=abs(s in (pi*ff).(pi*ff);plot(f,AW,'b - ',ff,AWW,'r:');title('Aliasing CaUSed by Undersampling')XIabel('f->')ylabel('W(f)'),lege nd('by FFTTTheoreticaI')由上述程序画出的“欠”采样时

19、引起的混叠曲线如图3.8.13所示：图3813 “欠”采样时引起的混叠例8 （例3.8.5 ）试分别求出含白噪声干扰的正弦信号与白噪声的自相关函数，并对它 们的结果进行比较。%example 3.8.5 MATLAB PRoGRAM N=1000;n=0:N-1;Fs=500;t=n/Fs;Lag=100;%x1=si n(2*pi*10*t)+0.6*ra ndn (1,le ngth(t);%c,lags=xcorr(x1,Lag,' Un biased')%subplot(2,2,1)%plot(t,x1)XIabel('t')ylabel('x1

20、(t)');title('含白噪声的正弦信号x1');gridsubplot(2,2,2)%plot(lagsFs,c)XIabel('t')ylabel('Rxx1(t)');相关信号的最大延迟量 含白噪声的正弦信号x1无偏自相关函数的计算 画出x1曲线画出x1自相关曲线title('x1 的自相关函数');gridx2=ra ndn (1,le ngth(t);%c,lags=xcorr(x2,Lag,' Un biased') %SUbPlot(2,2,3)%plot(t,x2)XIabel('

21、;t')ylabel('x2(t)');title(' 白噪声 x2');gridsubplot(2,2,4)%plot(lagsFs,c)XIabel('t')ylabel('Rxx2(t)');title('x2 的自相关函数');发生白噪声x2白噪声x2的无偏自相关函数画出x2曲线画出x2自相关曲线grid由上述程序画出的信号与自相关函数的曲线如图3.8.19所示:含白噪声的正弦信号Jd图3.8.19信号与自相关函数由计算结果和和图 3.8.19的曲线可以看出：同时具有周期性和白噪声干扰的信号，其自相

22、关函数不但在T= 0时具有最大值，而且在T较大时自相关函数仍有明显的周期性，它 的频率和原周期信号的周期相同；而无周期的白噪声，自相关函数在T= 0时也具有最大值，但当T稍大时迅速衰减至零附近，利用自相关函数的这一特性可用来识别随机信号中是否含 有周期成分和它的频率。 例9（例4.4.3）试用MATLAB语言绘制巴特沃思低通模拟滤波器的平方幅频相应曲线，滤波器的阶数分别为2, 5, 10, 20。解%例443的MATLAB 程序 n=0:0.01:2 for i=1:4SWitCh iCaSe 1N=2;CaSe 2N=5;CaSe 3N=10;CaSe 4N=20;%巴特沃思滤波器原型设计函

23、数，n:阶数;z,p,k:%滤波器零点，极点和增益%零极点增益转换为传递函数%模拟滤波器频率响应函数%幅度平方函数end z,p,k=buttap(N);b,a=zp2tf( z, p,k); H,w=freqs(b,a, n); magH2=(abs(H).2; hold On plot(w,magH2);axis(0201);endxlabel( W/wc 'ylabel( H(jw) 2'' grid图444 不同阶次巴特沃思滤波器的幅度平方函数 例10（例446 ）试用MATLAB程序，确定一个模拟低通滤波器的阶数N和截止频率WC。设计指标为：通带边界频率200

24、 ,阻带边界频率 300 ,通带波纹为1dB,在 Z处，幅度衰减大于 18 dB。解%设计切比雪夫低通滤波器的MATLAB程序wp=200*pi;ws=300*pi;Rp=1;%通带波纹Rs=16;%计算滤波器阶数ebs=sqrt(10(Rp10)-1);A=10(Rs20);Wc=wp;Wr=WSwp;g=sqrt(A*A-1)ebs;%阻带衰减%波纹系数%A为参变量%截止频率%g为参变量N仁 log10(g+sqrt(g*g-1)log10(Wr+sqrt(Wr*Wr-1);%滤波器阶数计算 N=CeiI(N1);%N 取整数运行上述程序后，可得滤波器的截止频率WC和阶数N为WC =628

25、.3185N =4 例11 (例444)绘制出阶数分别为 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);b,a=zp2tf(z,p,k); H,w=freqs(b,a, n); magH2=(abs(H).2;% OUtPUt%滤波器通带波纹系数%设计切比雪夫模拟低通滤波器原型函数,%乙p,

26、k分别为滤波器的零点、极点和增益%零点、极点和增益转换为传递函数%模拟滤波器的频率响应%幅度平方函数%定义字符串变量posplot= 1'0 'u m2str(i); subplot(posplot) plot(w,magH2, k'); ylabel( H(jw)2 ' title( N= 'u m2str(N);gridend其幅度平方函数曲线如图446所示。由图446可以看出：切比雪夫滤波器在通带内具有等波纹起伏特性，而在阻带内则单 调下降。阶次越高，特性越接近矩形。 例12(例447)试设计一巴特沃思模拟带通滤波器，设计要求为：通带频率2kHz3

27、kHz,两边的过渡带宽为0.5kHz ,若通带波纹IdB ,阻带衰减大于100dB.解%设计巴特沃思带通滤波器MA TLAB程序%滤波器设计指标wp=2000 3000*2*pi;ws=1500 3500*2*pi;Rp=1;RS=100;%计算阶数和截止频率N,Wn=buttord (wp,ws,Rp,Rs, ' s');FC=Wn/(2*pi);%计算滤波器传递函数多项式系数b,a=butter(N,Wn, ' S ' );%画出滤波器幅频特性w=li nspace(1,4000,1000)*2*pi;% 生成线性等间隔向量H=freqs(b,a,w);ma

28、gH=abs(H);PhaH=Un wrap(a ngle(H);Plot(W/(2*pi),20*log10(magH),' k ' );xlabel( '频率(Hz)');ylabel( '幅度(dB)' );title( '巴特沃思模拟滤波器')grid On由上述程序画出的巴特沃思模拟滤波器的幅频特性如图4411所示。图4411巴特沃思带通滤波器的幅频特性 例13 (例448 )试设计一个切比雪夫模拟带阻滤波器。要求的指标为：阻带上下边界频率为2KHz与3KHz ,两侧过渡带为 0.5KHZ ,通带波纹1dB,阻带衰减大于

29、 60 dB。 解%设计切比雪夫带阻滤波器MATLAB程序%滤波器设计指标wp=2000 3000*2*pi;ws=1500 3500*2*pi;Rp=1;Rs=60;%计算阶数和截止频率N,Wn=cheblord (wp,ws,Rp,Rs, ' s');%计算滤波器传递函数多项式系数b,a=chebyl(N,Rp,Wn, ' stop ' , ' S ' );%画出滤波器幅频特性W=Ii nspace(1,4000,1000)*2*pi;% 生成线性等间隔向量H=freqs(b,a,w);% 相位展开magH=abs(H);PhaH=Un wr

30、ap(a ngle(H);Plot(W/(2*pi),20*log10(magH),' k ' );xlabel( '频率(Hz)');ylabel( '幅度(dB)' );title( '切比雪夫模拟滤波器' )grid On由上述程序画出的切比雪夫模拟滤波器的幅频特性如图4412所示。图4.4.12切比雪夫带通滤波器的幅频特性 例14 （例5.3.4）一低通滤波器，其通带和阻带的技术指标分别是0.2 IRP 0.25 dB0.3 IRZ 50dB试用频率抽样法设计一FIRN=20 。具有线性相位的滤波器，取jkjk由 H(k)

31、 Hd(e N ) Hd(e 20 ) , k 0, 1 ,19当k=2时,正好是在通带边界频率处有一个频率抽样点，即P 0.2下一个抽样点为k=3 ,是阻带上边界频率Z ，阻带域通带间无过渡带，即Z 03则在通带0P内有三个抽样点，阻带上共有七个抽样点。从而有Hk 1,1,1,0,由N=20 ,相位常数(N 1)/2(201）/29.5,其相位可表示为9k0.95 k ,0.95 (20 k),0 k 910 k 19h(n)，并有频响内插公式最根据式(5.3.37)得到H(k),利用离散傅里叶变换，求得程序如下:后可得FlR滤波器的频响H(ej )。其MATLAB%例5.3.4中设计FIR

32、滤波器MATLAB程序N=20； alpha (N 1)/2;10: N 1 ; Wl (2*piN)*1HrS 1,1,1,ZerOSIl ,15),1,1;Hdr 1,1,0,0;Wdl 0,0.25,0.25,1;kl 0:floor(N 1)/2); k2 floor(N 1)/2) 1: N 1 ；angH = -alpha * ( 2 * Pi ) / N * kl ,alpha * (2 * Pi ) / N * ( N-k2 );H = Hrs. * exp (j * an gH);=real (ifft (H , N );db , mag , Pha ,grd , W = f

33、reqz -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 ('频率样本： N = 20')xlabel (频率(单位：Pi) ');ylabel ( Hr (k) JSet (gca , 'XtiCkMOde ' ,' manual',' XTiCk ',0,0.2,0.3,1 );Set (gca , YtiCk

34、MOde ' ,manual', Ytick',0,1); girdsubplot(2,2,2) ; Stem(I , 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) ,6');axis ( 0 , 1 , -0.2,1.2 ) ; title (振幅响应')xlabel (频率(单位:Pi ) '); ylable ( Hr(W)')Set (gca , XTiCkMOde ', XtiCk ', 0,0.2,0.3 ,