MATLAB信号处理仿真实验Word文档格式.docx
《MATLAB信号处理仿真实验Word文档格式.docx》由会员分享，可在线阅读，更多相关《MATLAB信号处理仿真实验Word文档格式.docx（27页珍藏版）》请在冰点文库上搜索。
1.1892-0.03760.32730.1746-0.18670.7258
（6）方波
t=0:
0.1*pi:
6*pi;
y=square（t）;
axis（[07*pi-1.51.5]）;
plot（t,y）;
xlabel（'
时间'
）;
ylabel（'
幅值'
（7）正弦波
t=0:
0.01*pi:
2*pi;
x=sin（2*pi*t）;
plot（t,x）;
（8）锯齿波
Fs=10000;
1/Fs:
1.5;
x=sawtooth（2*pi*50*t）;
axis（[00.2-11]）;
（9）基本非周期波形
1/1000:
2;
x=chirp（t,0,1,150）;
specgram（x,256,1000,256,250）;
（10）sinc信号
t=linspace（-5,5）;
x=sinc（t）;
（11）pulstran信号
1/50E3:
10E-3;
d=[0:
1/1E3:
0.8.^（0:
10）]'
;
x=pulstran（t,d,'
gauspuls'
10E3,0.5）;
（12）diric信号
t=[-4*pi:
0.1/pi:
4*pi];
x=diric（t,7）;
y=diric（t,8）;
subplot（1,2,1）;
title（'
n为奇数'
subplot（1,2,2）;
n为偶数'
2．有两信号分别为
，
，其中
，编程实现此二信号的叠加，并计算它的抽样和、抽样积、信号能量和信号功率。
信号叠加：
抽样和：
抽样积：
信号能量：
信号功率：
程序如下：
t=（0:
0.001:
1）;
xh=sin（2*pi*50*t）+2*sin（2*pi*120*t）;
x=xh+0.5*randn（size（t））;
y1=sum（x（1:
50））;
y2=rod（x（1:
50）;
Ex=sum（abs（x）.^2）;
Px=sum（abs（x）.^2/n）;
subplot（2,3,1）;
plot（t（1:
50）,xh（1:
Xh'
subplot（2,3,4）;
50）,x（1:
x'
subplot（2,3,2）;
plot（y1,'
抽样和'
subplot（2,3,3）;
plot（y2,'
抽样积'
subplot（2,3,5）;
plot（Ex,'
能量'
subplot（2,3,6）;
plot（Px,'
功率'
结果：
3．求信号
和
之和（其中
），并计算和的FFT变换。
有限长序列的离散傅立叶变换对定义：
正变换：
反变换：
1/99:
x=sin（2*pi*15*t）+sin（2*pi*40*t）;
y=fft（x）;
m=abs（y）;
p=unwrap（angle（y））;
f=（0:
length（y）-1）*99/length（y）;
plot（f,m）
set（gca,'
XTick'
[15406085]）;
grid;
plot（f,p*180/pi）
4．模拟信号
，求N点DFT的幅值谱和相位谱。
N=64;
63;
t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=fft（x,N）;
subplot（3,1,1）;
sourcesignal'
subplot（3,1,2）;
plot（q,abs（y））;
magnitude'
subplot（3,1,3）;
plot（q,angle（y））;
phase'
三、实验报告要求
完成上述各题给定要求并做好实验记录。
3.2IIR数字滤波器的设计
学会用MATLAB设计IIR数字滤波器。
1．设计一个三阶的模拟低通滤波器，通带内最大衰减为3dB，在阻带内的最大衰减为40dB，截止频率为
弧度，再把它转换成截止频率为
弧度的高通滤波器，绘出它们的频率响应图。
[z,p,k]=ellipap（3,3,40）;
[A,B,C,D]=zp2ss（z,p,k）;
[At,Bt,Ct,Dt]=lp2lp（A,B,C,D,8*pi）;
[num,den]=ss2tf（At,Bt,Ct,Dt）;
figure;
freqs（num,den）;
[At1,Bt1,Ct1,Dt1]=lp2hp（A,B,C,D,50*pi）;
[num1,den1]=ss2tf（At1,Bt1,Ct1,Dt1）;
freqs（num1,den1）;
（分别为三阶椭圆低通滤波器和高通滤波器的频率响应）
2．设计一个数字信号处理系统，它的抽样频率为Fs＝100Hz，在该系统设计一个Butter-worth型高通滤波器，使其在通带内允许的最小衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。
wp=30*2*pi;
ws=40*2*pi;
rp=0.5;
rs=40;
Fs=100;
%把数字滤波器的频率特征转换成模拟滤波器的频率特征
[N,Wn]=buttord（wp,ws,rp,rs,'
s'
%选择滤波器的阶数
[z,p,k]=buttap（N）;
%创建Butterworth低通滤波器原型
%表达形式从零极点增益形式转换成状态方程形式
[At,Bt,Ct,Dt]=lp2hp（A,B,C,D,Wn）;
%实现低通到高通滤波器类型的转换
[num1,den1]=ss2tf（At,Bt,Ct,Dt）;
[num2,den2]=bilinear（num1,den1,100）;
%从模拟高通到数字高通的转换，采用双线性变换
[H,W]=freqz（num2,den2）;
%绘出频率响应
plot（W*Fs/（2*pi）,abs（H））;
grid
频率/Hz'
3．模拟信号
，由系统
处理，采样频率为1000Hz。
（1）设计一个最小阶数的IIR数字滤波器，以小于1dB的衰减通过150Hz的分量，以至少40dB抑制100Hz的分量。
滤波器应有单调的通带和等波动的阻带，求出有理函数形式的系统函数，画出幅度响应（dB）。
（2）产生上述信号
的150个样本，通过上述所设计的滤波器得到输出序列，内插此序列得到
。
画出输入和输出信号的图并解释结果。
fp=150;
fr=100;
fs=1000;
wp=2*pi*fp/fs;
wr=2*pi*fr/fs;
Ap=1;
Ar=40;
[N,wn]=cheb2ord（wp/pi,wr/pi,Ap,Ar）;
[b,a]=cheby2（N,Ar,wn,'
high'
[C,B,A]=dir2cas（b,a）;
[db,mag,pha,w]=freqz_m（b,a）;
subplot（411）;
plot（w/pi,db）;
axis（[0,1,-50,7]）;
149;
t=n/fs;
x=5*sin（2*pi*fr*t）+2*cos（2*pi*fp*t）;
subplot（412）;
axis（[0,0.15,-7,7]）;
y=filter（b,a,x）;
subplot（413）;
stem（y）;
ya=y*sinc（fs*（ones（length（n）,1）*t-（n/fs）'
*ones（1,length（t））））;
subplot（414）;
plot（t,ya）;
axis（[0,0.15,-2,2]）;
%直接型转换成级联型
function[b0,B,A]=dir2cas（b,a）;
b0=b
（1）;
b=b/b0;
a0=a
（1）;
a=a/a0;
b0=b0/a0;
M=length（b）;
N=length（a）;
ifN>
M
b=[bzeros（1,N-M）];
elseifM>
N
a=[azeros（1,M-N）];
N=M;
else
NM=0;
end
k=floor（N/2）;
B=zeros（k,3）;
A=zeros（k,3）;
ifk*2==N;
b=[b0];
a=[a0];
broots=cplxpair（roots（b））;
aroots=cplxpair（roots（a））;
fori=1:
2:
2*k
Brow=broots（i:
1:
i+1,:
Brow=real（poly（Brow））;
B（fix（i+1）/2,:
）=Brow;
Arow=aroots（i:
Arow=real（poly（Arow））;
A（fix（i+1）/2,:
）=Arow;
function[db,mag,pha,w]=freqz_m（b,a）;
[H,w]=freqz（b,a,1000,'
whole'
H=（H（1:
501））'
w=（w（1:
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
结果为：
C=
0.2942
B=
1.0000-1.62841.0000
1.0000-1.79191.0000
1.0000-1.97081.0000
A=
1.0000-1.00320.2706
1.0000-1.02130.4255
1.0000-1.13730.7517
由图可见，输出信号抑制了输入信号的100Hz分量，基本上是150Hz的余弦信号，达到了信号处理的目的。
3.3FIR数字滤波器的设计
一、实验目的
学会用MATLAB设计FIR数字滤波器。
1．设理想带阻滤波器频率响应为
利用Kaiser窗函数，设计长度为55的带阻滤波器，使阻带衰减为60dB。
参数
可由式
确定，其中
n=55-1;
w=[0.40.6];
beta=0.1102*（60-8.7）;
kai_w=kaiser（n+1,beta）;
b=fir1（n,w,'
stop'
kai_w）;
[h,w]=freqz（b,1,512,2）;
plot（w,20*log10（abs（h）））;
频率（归一化）'
）,ylabel（'
幅度（dB）'
2．根据下列技术指标，设计一个数字FIR低通滤波器。
选择一个恰当的窗函数，确定单位脉冲响应，绘出所设计的滤波器的幅度响应。
理想低通数字滤波器的频率响应为
式中，
为截止频率，
为采样延迟。
理想低通数字滤波器的单位脉冲响应为
为无限长非因果序列，关于
偶对称。
为了从
得到一个FIR数字滤波器，必须同时在两边截取
，要得到一个因果的线性相位FIR滤波器，它的
的长度为N，必须有
这种截取可看作是
矩形窗
为关于
偶对称的有限长因果序列。
N为奇数时是1型，N为偶数时是2型。
根据窗函数最小阻带衰减的特性，只有海明窗和布莱克曼窗可提供大于50dB的衰减。
故选择海明窗，它提供较小的过渡带，因此具有较小的阶数。
wp=0.2*pi;
wr=0.4*pi;
tr_width=wr-wp;
%过渡带宽度
N=ceil（6.6*pi/tr_width）+1%滤波器的长度，N＝奇数为1型；
N＝偶数为2型
wc=（wr+wp）/2;
%理想低通滤波器的截止频率
hd=ideal_lp（wc,N）;
%理想低通滤波器的单位脉冲响应
w_ham=（hamming（N））'
%海明窗
h=hd.*w_ham;
%截取得到实际单位脉冲响应
[db,mag,pha,w]=freqz_m（h,[1]）;
%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
Ap=-（min（db（1:
wp/delta_w+1）））;
%实际通带波动
Ar=-round（max（db（wr/delta_w+1:
501）））;
%最小阻带衰减
subplot（221）;
stem（n,hd）;
理想单位脉冲响应hd（n）'
）
subplot（222）;
stem（n,w_ham）;
海明窗w（n）'
subplot（223）;
stem（n,h）;
实际单位脉冲响应h（n）'
subplot（224）;
幅度响应（dB）'
axis（[0,1,-100,10]）;
functionhd=ideal_lp（wc,N）;
alpha=（N-1）/2;
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
N=34
Ap=0.0477
Ar=52
滤波器的长度为34，实际通带波动0.0477dB，最小阻带衰减52dB，满足设计要求。
3．设抽样频率为Fs＝1000Hz，已知原信号为
，由于某种原因，信号被白噪声污染，实际获得的信号为xn＝x＋randn（size（t）），要求设计一个FIR滤波器恢复出原始信号。
%生成相应的信号
Fs=1000;
t=1:
1/fs:
%2秒长度的序列
x=sin（2*pi*80*t）+2*sin（2*pi*140*t）;
xn=x+randn（size（t））;
%使用最小二乘的滤波器设计方法设计一个多带滤波器
%取滤波器的阶数为100
n=100;
f=[00.130.150.170.190.250.270.290.311];
m=[0011001100];
b=firls（n,f,m）;
[H,W]=freqz（b,1,512,2）;
plot（W,abs（H））;
归一化频率（Nyquist频率＝１）'
幅度'
%对xn进行滤波
xo=filter（b,1,xn）;
nn=500:
750;
tt=nn/Fs;
subplot（311）;
plot（tt,x（nn））;
原始信号'
）,grid;
subplot（312）;
plot（tt,xn（nn））;
污染信号'
subplot（313）;
plot（tt,xo（nn））;
滤波信号'
时间（秒）'
三、实验报告要求
升级会员 