MATLAB信号处理仿真实验.docx
《MATLAB信号处理仿真实验.docx》由会员分享,可在线阅读,更多相关《MATLAB信号处理仿真实验.docx(27页珍藏版)》请在冰点文库上搜索。
![MATLAB信号处理仿真实验.docx](https://file1.bingdoc.com/fileroot1/2023-5/7/b49b3d88-c68c-48a5-a9dd-f209a5e7da68/b49b3d88-c68c-48a5-a9dd-f209a5e7da681.gif)
MATLAB信号处理仿真实验
3信号处理仿真实验(教师)
3.1MATLAB信号处理基础实验
一、实验目的
1.掌握MATLAB常用信号处理波形;
2.学习信号序列的各种操作;
3.学习离散傅立叶变换的MATALB相关操作。
二、实验内容
1.产生下列各种波形,并记录结果。
(1)单位抽样序列
x=[1zeros(1,n-1)]
>>n=7;
>>x=[1zeros(1,n-1)]
x=
1000000
(2)单位阶跃矩阵
x=ones(1,N)
>>N=6;
>>x=ones(1,N)
x=
111111
(3)实指数序列
n=0:
N-1;
x=a.^n;
>>a=6;
>>N=6;
>>n=0:
N-1;
>>x=a.^n
x=
163621612967776
(4)复指数序列
n=0:
N-1;
x=ex((lu+j*w0)*n);
>>N=6;
lu=3;
w0=40;
n=0:
N-1;
x=exp((lu+j*w0)*n)
x=
1.0e+006*
0.0000-0.0000+0.0000i-0.0000-0.0004i0.0066+0.0047i-0.1588+0.0357i1.5926-2.8548i
(5)随机序列
rand(1,N);
randn(1,N);
>>N=6;
>>rand(1,N)
ans=
0.45650.01850.82140.44470.61540.7919
>>randn(1,N)
ans=
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);
xlabel('时间');
ylabel('幅值');
(8)锯齿波
Fs=10000;
t=0:
1/Fs:
1.5;
x=sawtooth(2*pi*50*t);
plot(t,x);
axis([00.2-11]);
(9)基本非周期波形
t=0:
1/1000:
2;
x=chirp(t,0,1,150);
specgram(x,256,1000,256,250);
(10)sinc信号
t=linspace(-5,5);
x=sinc(t);
plot(t,x);
(11)pulstran信号
t=0:
1/50E3:
10E-3;
d=[0:
1/1E3:
10E-3;0.8.^(0:
10)]';
x=pulstran(t,d,'gauspuls',10E3,0.5);
plot(t,x);
(12)diric信号
t=[-4*pi:
0.1/pi:
4*pi];
x=diric(t,7);
y=diric(t,8);
subplot(1,2,1);
title('n为奇数');
plot(t,x);
subplot(1,2,2);
plot(t,y);
title('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:
50));
title('Xh');
subplot(2,3,4);
plot(t(1:
50),x(1:
50));
title('x');
subplot(2,3,2);
plot(y1,'x');
title('抽样和');
subplot(2,3,3);
plot(y2,'x');
title('抽样积');
subplot(2,3,5);
plot(Ex,'x');
title('能量');
subplot(2,3,6);
plot(Px,'x');
title('功率');
结果:
3.求信号
和
之和(其中
,
),并计算和的FFT变换。
有限长序列的离散傅立叶变换对定义:
正变换:
反变换:
程序如下:
t=(0:
1/99:
1);
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);
subplot(1,2,1);
plot(f,m)
set(gca,'XTick',[15406085]);grid;
subplot(1,2,2);
plot(f,p*180/pi)
set(gca,'XTick',[15406085]);grid;
结果:
4.模拟信号
,求N点DFT的幅值谱和相位谱。
程序如下:
N=64;
n=0:
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);
plot(t,x);
title('sourcesignal');
subplot(3,1,2);
plot(q,abs(y));
title('magnitude');
subplot(3,1,3);
plot(q,angle(y));
title('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低通滤波器原型
[A,B,C,D]=zp2ss(z,p,k);
%表达形式从零极点增益形式转换成状态方程形式
[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
xlabel('频率/Hz');
ylabel('幅值');
结果:
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]);
n=0:
149;t=n/fs;
x=5*sin(2*pi*fr*t)+2*cos(2*pi*fp*t);
subplot(412);plot(t,x);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];
end
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:
1:
i+1,:
);
Arow=real(poly(Arow));
A(fix(i+1)/2,:
)=Arow;
end
function[db,mag,pha,w]=freqz_m(b,a);
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:
1:
501))';
w=(w(1:
1:
501))';
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)));grid;
xlabel('频率(归一化)'),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型
n=0:
1:
N-1;
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:
1:
wp/delta_w+1)));%实际通带波动
Ar=-round(max(db(wr/delta_w+1:
1:
501)));%最小阻带衰减
subplot(221);stem(n,hd);title('理想单位脉冲响应hd(n)')
subplot(222);stem(n,w_ham);title('海明窗w(n)')
subplot(223);stem(n,h);title('实际单位脉冲响应h(n)')
subplot(224);plot(w/pi,db);title('幅度响应(dB)')
axis([0,1,-100,10]);
functionhd=ideal_lp(wc,N);
alpha=(N-1)/2;
n=0:
1:
N-1;
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;%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));
grid
xlabel('归一化频率(Nyquist频率=1)');
ylabel('幅度');
%对xn进行滤波
xo=filter(b,1,xn);
figure;
nn=500:
750;
tt=nn/Fs;
subplot(311);
plot(tt,x(nn));
ylabel('原始信号'),grid;
subplot(312);
plot(tt,xn(nn));
ylabel('污染信号'),grid;
subplot(313);
plot(tt,xo(nn));
ylabel('滤波信号'),grid;
xlabel('时间(秒)');
三、实验报告要求
完成上述各题给定要求并做好实验记录。