MATLAB信号处理例题.docx

1、MATLAB信号处理例题例1设方波的数学模型为： ，基频： 用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;

2、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次谐波合成plot(t,yn,r) %用红色实线画出n次谐波合成 从图中我们可以看

3、到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

4、=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，固有频率Hz，采样间隔s,采样点数N256。试计算理论幅频特性与由系统阶跃响应计算出的幅频特性数据值，并画出两

5、个计算结果的幅频特性曲线。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:N Rw(i,1)=1-Im(i,1)*(i-1)*dw*dt; Iw(i,1)=Re(i,1)*(i-1)*dw*dt;endffw=Rw+Iw*sqr

6、t(-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。 例6.2.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+9

7、0*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 若有信号为 式中，为白噪声（用MATLAB中的函数产生）。设采样频率；试用周期法并应用MATLAB编程计算，当数据长度分别为和 两种情况下上述信号的功率谱。%例6.

8、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=

9、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上均匀

10、分布的N点随机序列，表示为：x(n)rand(1,N),h(n)为均值为零、方差为1的N点高斯分布随机序列，表示为：h(n)rand(1,L),试求1N150时的卷积并比较其运算时间。%例3.8.1 直接卷积和快速卷积的比较%conv_time=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:100 h=randn(1,N); %产生两个随机序列 x=rand(1

11、,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_

12、time(25:150),n(25:150),fft_time(25:150)上述两种卷积的计算时间的比较如图3.8.5所示。图3.8.5 两种卷积计算时间比较例7（例3.8.2 ）运用FFT求取矩形脉冲 的谱，说明采样频率低引起的混叠现象。（1）先编写有一定通用性的函数文件cftbyfft.m cftbyfft.mfunction AW,f=cftbyfft(wt,t,flag)%cftbyfft.m%本程序采用FFT计算连续时间Fouie变换。输出幅频数据对(f,AW)。%输入量(wt,t)为已经窗口化了的时间函数wt(t)，它们分别是长度为N的向量。%对于时限信号，应使该取值时段与窗口长

13、度相比足够小。以提高频率分辨率。%对于非时限信号，窗口长度的选取应使窗口外的函数值小%到可忽略，以提高近似精度。%输入宗量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)%把以上计算

14、结果改写到-fs/2,fs/2范围if flag=0 n=-N/2:(N/2-1); W=fftshift(W); %产生周期序列的频谱endf=n*df; %频率分量向量AW=abs(W); %幅频谱数据向量if nargout=0 plot(f,AW);grid,xlabel(频率f);ylabel(|w(f)|)end（2）运行以下指令，绘制时域波形和幅频谱M=5; %做2的幂次用。本例把M设得较小，是为了观察混叠。tend=1; %波形取非零值的时间长度。T=10; %窗口化长度应足够大，以减小窗口化引起的泄露“旁瓣”效应。N=2M; %采样点数，取2的幂是为了使FFT运算较快。dt=

15、T/N; %以上T、N的取值应使N/Tfs采样频率大于两倍时间波形 %带宽，以克服采样引起的频谱混叠。 %在本例中，据理论分析知W(f=7.5)=Sa(7.5*pi)=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所示：图3.8.12 被窗口化的时间波形AW,f=cftbyfft(w,t,0);ff=f+eps; %为避免下面指令出

16、现0/0而采取的措施AWW=abs(sin(pi*ff)./(pi*ff);plot(f,AW,b-,ff,AWW,r:);title(Aliasing caused by undersampling)xlabel(f-);ylabel(|W(f)|),legend(by FFT,Theoretical)由上述程序画出的“欠”采样时引起的混叠曲线如图3.8.13所示：图3.8.13 “欠”采样时引起的混叠例8（例3.8.5） 试分别求出含白噪声干扰的正弦信号与白噪声的自相关函数，并对它们的结果进行比较。%example 3.8.5 MATLAB PROGRAMN=1000;n=0:N-1;Fs

17、=500;t=n/Fs;Lag=100; %相关信号的最大延迟量x1=sin(2*pi*10*t)+0.6*randn(1,length(t); %含白噪声的正弦信号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,

18、length(t); %发生白噪声x2c,lags=xcorr(x2,Lag,unbiased) %白噪声x2的无偏自相关函数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所示：图3.8.19 信号与自相关函数由计算结果和和图3.8.19的曲线可以看出：同时具有周期性和白噪声干

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

20、buttap(N); %巴特沃思滤波器原型设计函数，n:阶数;z,p,k: %滤波器零点，极点和增益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，阻带边界频率30

21、0，通带波纹为1dB，在Z处，幅度衰减大于18 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= 6

22、28.3185N 4例11（例4.4.4） 绘制出阶数分别为2，4，6，8的切比雪夫模拟低通滤波器的平方幅频响应曲线。解 写出其MATLAB程序如下：绘制切比雪夫低通滤波器幅度平方函数的MATLAB程序n=0：0.01：2；for i=1:4 switch i case 1 N=2； case 2 N=4； case 3 N=6； case 4 N=8； endRp=1; 滤波器通带波纹系数z,p,k=cheblap(N,Rp); 设计切比雪夫模拟低通滤波器原型函数， z,p,k分别为滤波器的零点、极点和增益b,a=zp2tf(z,p,k); 零点、极点和增益转换为传递函数H,w=freqs(

23、b,a,n); 模拟滤波器的频率响应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，阻带

24、衰减大于100dB.解 设计巴特沃思带通滤波器MATLAB程序滤波器设计指标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=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),

25、k);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,Wn=cheblord

26、 (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）);title(切比雪夫模拟滤波器)grid on由上述程序画出的切比雪夫模拟滤波器的幅频特性如图4.4.12所示。 图4.4.12 切比雪夫带通滤波器的幅

27、频特性例14（例5.3.4） 一低通滤波器，其通带和阻带的技术指标分别是 试用频率抽样法设计一FIR 具有线性相位的滤波器，取 N=20 。解 由 当k=2 时，正好是在通带边界频率 处有一个频率抽样点，即 下一个抽样点为k=3 ，是阻带上边界频率 ，阻带域通带间无过渡带，即 则在通带 内有三个抽样点，阻带上共有七个抽样点。从而有 由N=20 ，相位常数 其相位可表示为根据式（5.3.37）得到H(k) ，利用离散傅里叶变换，求得h(n) ，并有频响内插公式最后可得FIR滤波器的频响 。其MATLAB 程序如下：%例5.3.4中设计FIR滤波器MATLAB程序N=20；angH = -alph

28、a * ( 2 * pi ) / N * kl ,alpha * (2 * pi ) / N * ( N-k2 ) ;H = Hrs. * exp (j * angH) ;= 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 (gca ,YtickMode,manual,Ytick,0 , 1); girdsubplot(2 , 2 , 2) ; stem(l , h) ; axis(-1 , N ,-0.1 , 0.3)title (单位抽样响应 ) ; ylabe