MATLAB信号处理例题.docx
《MATLAB信号处理例题.docx》由会员分享，可在线阅读，更多相关《MATLAB信号处理例题.docx（25页珍藏版）》请在冰点文库上搜索。
MATLAB信号处理例题
◆例1设方波的数学模型为：
，基频：
用MATLAB软件完成该方波的合成设计
◆MATLAB源程序
t=-10:
0.1:
10;%设定一个数组有201个点，方波周期为20
e=5;w=pi/10;%设定方波幅值为5，w代表w0
m=-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决定方波的合成次数，在此给定50
yn=0;%设置初始值
fori=1:
n
yn=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，固有频率
Hz，采样间隔
s,采样点数N＝256。
试计算理论幅频特性与由系统阶跃响应计算出的幅频特性数据值，并画出两个计算结果的幅频特性曲线。
％example3.8.3MATLABPROGRAM
N=256;
dt=0.0004
wn=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）;
fori=1:
N
Rw（i,1）=1-Im（i,1）*（i-1）*dw*dt;
Iw（i,1）=Re（i,1）*（i-1）*dw*dt;
end
ffw=Rw+Iw*sqrt（-1）;
Aw=abs（ffw）
semilogx（w,20*log10（mag）,'r-'）
axis（[100,10000,-30,10]）
text（600,12,'幅频特性'）
holdon
semilogx（w,20*log10（Aw））
axis（[100,10000,-30,10]）
gridon
◆例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+90*pi/180）;
[c,lags]=xcorr（x,y,Lag,'unbiased'）;
subplot（2,1,1）
plot（t,x,'r'）
holdon
plot（t,y,'b'）
xlabel（'t'）;ylabel（'x（t）y（t）'）;
title（'原周期信号'）
grid
holdoff
subplot（2,1,2）;
plot（lags/Fs,c,'k'）;
xlabel（'t'）;ylabel（'Rxy（t）'）;
title（'互相关函数'）;
grid
◆例5
例6.2.5若有信号为
式中，
，
为白噪声（用MATLAB中的函数产生）。
设采样频率
；试用周期法并应用MATLAB编程计算，当数据长度分别为
和
两种情况下上述信号的功率谱。
%例6.2.5中周期图法计算信号功率谱的MATLAB程序
clf
Fs=2000;
%情况1:
数据长度N1=256
N1=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=1024
N2=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）;
%
forN=1:
150
tc=0;tf=0;%初始化
L=2*N-1;%加长序列长度
nu=round（（log10（L）/log10
（2））+0.45）;
L=2^nu;%使点数成为2的幂次
forI=1:
100
h=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）先编写有一定通用性的函数文件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）。
ifnargin==2;flag=1;end
N=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]范围
ifflag==0
n=-N/2:
（N/2-1）;
W=fftshift（W）;%产生周期序列的频谱
end
f=n*df;%频率分量向量
AW=abs（W）;%幅频谱数据向量
ifnargout==0
plot（f,AW）;grid,xlabel（'频率f'）;ylabel（'|w（f）|'）
end
（2）运行以下指令，绘制时域波形和幅频谱
M=5;%做2的幂次用。
本例把M设得较小，是为了观察混叠。
tend=1;%波形取非零值的时间长度。
T=10;%窗口化长度应足够大，以减小窗口化引起的泄露“旁瓣”效应。
N=2^M;%采样点数，取2的幂是为了使FFT运算较快。
dt=T/N;%以上T、N的取值应使N/T＝fs采样频率大于两倍时间波形
%带宽，以克服采样引起的频谱混叠。
%在本例中，据理论分析知W（f=7.5）=Sa（7.5*pi）=1/（7.5
%*pi）<5%.因此，可近似认为本例时间信号带宽为7.5Hz.
n=0:
N-1;%采样序列
t=n*dt;%采样点时间序列
w=zeros（size（t,2）,1）;
Tow=find（（tend-t）>0）;%产生非零波形时段的相应序列
w（Tow,1）=ones（length（Tow）,1）;%在窗口时段内定义的完整波形
plot（t,w,'b','LineWidth',2.5）,title（'TimeWaveform'）;xlabel（'t--->'）
由上述程序画出的被窗口化的时间波形如图3.8.12所示：
图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（'Aliasingcausedbyundersampling'）
xlabel（'f-->'）;ylabel（'|W（f）|'）,legend（'byFFT','Theoretical'）
由上述程序画出的“欠”采样时引起的混叠曲线如图3.8.13所示：
图3.8.13“欠”采样时引起的混叠
◆例8（例3.8.5）试分别求出含白噪声干扰的正弦信号与白噪声的自相关函数，并对它们的结果进行比较。
%example3.8.5MATLABPROGRAM
N=1000;
n=0:
N-1;
Fs=500;
t=n/Fs;
Lag=100;%相关信号的最大延迟量
x1=sin（2*pi*10*t）+0.6*randn（1,length（t））;%含白噪声的正弦信号x1
[c,lags]=xcorr（x1,Lag,'unbiased'）%无偏自相关函数的计算
subplot（2,2,1）%画出x1曲线
plot（t,x1）
xlabel（'t'）;ylabel（'x1（t）'）;
title（'含白噪声的正弦信号x1'）;
grid
subplot（2,2,2）%画出x1自相关曲线
plot（lags/Fs,c）
xlabel（'t'）;ylabel（'Rxx1（t）'）;
title（'x1的自相关函数'）;
grid
x2=randn（1,length（t））;%发生白噪声x2
[c,lags]=xcorr（x2,Lag,'unbiased'）%白噪声x2的无偏自相关函数
subplot（2,2,3）%画出x2曲线
plot（t,x2）
xlabel（'t'）;ylabel（'x2（t）'）;
title（'白噪声x2'）;
grid
subplot（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的曲线可以看出：
同时具有周期性和白噪声干扰的信号，其自相关函数不但在τ＝0时具有最大值，而且在τ较大时自相关函数仍有明显的周期性，它的频率和原周期信号的周期相同；而无周期的白噪声，自相关函数在τ＝0时也具有最大值，但当τ稍大时迅速衰减至零附近，利用自相关函数的这一特性可用来识别随机信号中是否含有周期成分和它的频率。
◆例9（例4.4.3）试用MATLAB语言绘制巴特沃思低通模拟滤波器的平方幅频相应曲线，滤波器的阶数分别为2，5，10，20。
解
％例4.4.3的MATLAB程序
n=0:
0.01:
2
fori=1:
4
switchi
case1
N=2;
case2
N=5;
case3
N=10;
case4
N=20;
end
[z,p,k]=buttap（N）;%巴特沃思滤波器原型设计函数，n:
阶数;z,p,k:
%滤波器零点，极点和增益
[b,a]=zp2tf（z,p,k）;%零极点增益转换为传递函数
[H,w]=freqs（b,a,n）;%模拟滤波器频率响应函数
magH2=（abs（H））.^2;%幅度平方函数
holdon
plot（w,magH2）;
axis（[0201]）;
end
xlabel（‘w/wc’）;
ylabel（‘∣H（jw）∣^2’’）;
grid
图4.4.4不同阶次巴特沃思滤波器的幅度平方函数
◆例10（例4.4.6）试用MATLAB程序，确定一个模拟低通滤波器的阶数N和截止频率Wc。
设计指标为：
通带边界频率200π，阻带边界频率300π，通带波纹为1dB，在ΩZ处，幅度衰减大于18dB。
解
％设计切比雪夫低通滤波器的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.3185
N＝
4
◆例11（例4.4.4）绘制出阶数分别为2，4，6，8的切比雪夫模拟低通滤波器的平方幅频响应曲线。
解写出其MATLAB程序如下：
％绘制切比雪夫低通滤波器幅度平方函数的MATLAB程序
n=0：
0.01：
2；
fori=1:
4
switchi
case1
N=2；
case2
N=4；
case3
N=6；
case4
N=8；
end
Rp=1;％滤波器通带波纹系数
[z,p,k]=cheblap（N,Rp）;％设计切比雪夫模拟低通滤波器原型函数，
％z,p,k分别为滤波器的零点、极点和增益
[b,a]=zp2tf（z,p,k）;％零点、极点和增益转换为传递函数
[H,w]=freqs（b,a,n）;％模拟滤波器的频率响应
magH2=（abs（H））.^2;％幅度平方函数
％Output
posplot=[‘10’num2str（i）];％定义字符串变量
subplot（posplot）
plot（w,magH2,‘k’）;
ylabel（‘|H（jw）^2|’）;
title（[‘N=’num2str（N）]）;
grid
end
其幅度平方函数曲线如图4.4.6所示。
由图4.4.6可以看出：
切比雪夫滤波器在通带内具有等波纹起伏特性，而在阻带内则单调下降。
阶次越高，特性越接近矩形。
◆例12（例4.4.7）试设计一巴特沃思模拟带通滤波器，设计要求为：
通带频率2kHz~3kHz,两边的过渡带宽为0.5kHz，若通带波纹1dB，阻带衰减大于100dB.
解
％设计巴特沃思带通滤波器MATLAB程序
％滤波器设计指标
wp=[20003000]*2*pi；
ws=[15003500]*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）,‵k‵）;
xlabel（‵频率（Hz）‵）；
ylabel（‵幅度（dB）‵）;
title（‵巴特沃思模拟滤波器‵）
gridon
由上述程序画出的巴特沃思模拟滤波器的幅频特性如图4.4.11所示。
图4.4.11巴特沃思带通滤波器的幅频特性
◆例13（例4.4.8）试设计一个切比雪夫模拟带阻滤波器。
要求的指标为：
阻带上下边界频率为2KHz与3KHz，两侧过渡带为0.5KHz，通带波纹1dB,阻带衰减大于60dB。
解
％设计切比雪夫带阻滤波器MATLAB程序
％滤波器设计指标
wp=[20003000]*2*pi；
ws=[15003500]*2*pi；
Rp=1;
Rs=60;
%计算阶数和截止频率
[N,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）‵）;
title（‵切比雪夫模拟滤波器‵）
gridon
由上述程序画出的切比雪夫模拟滤波器的幅频特性如图4.4.12所示。
图4.4.12切比雪夫带通滤波器的幅频特性
◆例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=[-alpha*（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（‘频率样本：
N＝20’）
xlabel（‘频率（单位：
pi）’）;ylabel（‘Hr（k）’）
set（gca,‘XtickMode’,‘manual’,‘XTick’,[0,0.2,0.3,1]）;
set（gca,‘YtickMode’,‘manual’,‘Ytick’,[0,1]）;gird
subplot（2,2,2）;stem（l,h）;axis（[-1,N,-0.1,0.3]）
title（‘单位抽样响应’
升级会员 