一、数字信号处理(确定性信号)
1、对于一个LTI 系统,设其输入序列为矩形冲激信号x(n)=u(n)-u(n-10),而冲激相应为)(9.0)(n u n h n ,用MATLAB 求解输出信号。可以直接调用卷积函数来实现。
解:
clear all
x=[1,1,1,1,1,1,1,1,1,1]; n=[0:9]; y=0.9.^n; z=conv(x,y); N=length(z); stem(0:N-1,z);
图像如下:
2、编程求两个序列之间的相关系数。设序列x (k )={3,11,7,0,-1,4,2},n=[-3,-2,-1,0,1,2,3],将x 进行移位再加上一个白噪声信号,即y(k)=x(k-2)+w(k),其中k 属于n ,需要计算x 序列与y 序列之间的相关系数,可以使用卷积来实现。
解:
clear all
>> x=[3,11,7,0,-1,4,2]; >> nx=[-3:3];
>> [y,ny]=sigshift(x,nx,2); >> w=randn(1,length(y)); >> nw=ny;
>> [y,ny]=sigadd(y,ny,w,nw); >> [x,nx]=sigfold(x,nx);
>> [rxy,nrxy]=conv_m(y,ny,x,nx); >> subplot(1,1,1); >> stem(nrxy,rxy)
>> axis([-5,10,-50,250]); >> xlabel('延迟量1'); >> ylabel('rxy');
>> title('噪声序列的互相关')
图像如下:
3、利用filter 函数计算冲激相应和单位阶跃响应。设离散系统由下列差分方程表示:)()2(9.0)1()(n x n y n y n y =-+--。
解:
冲激响应: clear all
a1=[1,-1,0.9]; b1=[1]; n=0:100; x1=[1 zeros(1,100)]; %补零 y1filter=filter(b1,a1,x1); stem(n,y1filter); title('冲激响应'); xlabel('x'); ylabel('y');
阶跃响应:Array clear all
a1=[1,-1,0.9]; b1=[1];
n=0:100;
x2=ones(1,101); %全一矩阵y1filter=filter(b1,a1,x2); stem(n,y1filter);
title('阶跃响应');
xlabel('x');
ylabel('y');
图像如下:
4、编程求解有限时宽复指数序列n j e n x )9.0()(3/π=,n=[0,10]的频谱及能量谱。
解:
n=0:10;x=(0.9*exp(j*pi/3)).^n; k=-200:200;w=(pi/100)*k; X=x*(exp(-j*pi/100)).^(n'*k); magX=abs(X);angX=angle(X); subplot(2,1,1);
plot(w/pi,magX);grid
xlabel('以w/pi 为单位的频率'); ylabel('|X|') title('幅度部分') subplot(2,1,2);
plot(w/pi,angX/pi);grid
xlabel('以w/pi 为单位的频率'); ylabel('弧度/pi') title('相角部分') figure;
Y=sum(abs(x).^2); stem(Y);
title('能量谱');
给定信号x (n)
对x 进行离散时间傅里叶变换得到X
设定X 的幅值和相角
函数
绘图
图像如下:
5、已知信号()2sin(4)5cos(8)x t t t ππ=+,随机正态白噪声为()w t 的均值为0,标准误差为1。试比较采样点数分别为45、64点时()()x t w t +的频谱(幅度谱)图像的差异。
解:
clear all; fs=100; N=45; n=0:N-1; t=n/fs;
x=2*sin(4*pi*t)+5*cos(8*pi*t); y=x+randn(size(x)); X=fft(y,N); magX=abs(X); f=n*fs/N; subplot(2,1,1); plot(f,magX); xlabel('频率Hz'); ylabel('45点幅度谱');
clear all; fs=100; N=64; n=0:N-1; t=n/fs;
x=2*sin(2*pi*2*t)+5*cos(2*pi*4*t);
y=x+randn(size(x)); X=fft(y,N);
magX=abs(X);
f=n*fs/N;
subplot(2,1,2);
plot(f,magX); xlabel('频率Hz'); ylabel('64点幅度谱');
6、设序列)52.0cos()48.0cos()(n n n x ππ+=,分别考虑n=[0,10]和n=[0,100]两种情况频谱特征。对比当n=[0,10]时,取10点及取100点(补90个0点)时情况,并说明高密度和高分辨率频谱之间的差异。
解:
n=[0:100];
x=cos(0.48*pi*n)+cos(0.52*pi*n); n1=[0:10]; x1=x(1:1:10); Y1=fft(x1,10); magY1=abs(Y1); subplot(3,1,1); plot(magY1);
title('10点的幅度谱'); n2=[0:100];
x2=[x(1:1:10),zeros(1,90)]; Y2=fft(x2,100); magY2=abs(Y2); subplot(3,1,2); plot(magY2);
title('补90个0的幅度谱(高密度)');
n3=[0:100];
x3=x(1:1:100);
Y3=fft(x3,100);
magY3=abs(Y3);
subplot(3,1,3);
plot(magY3);
title('100点的幅度谱(高分辨率)'); 图像如下:
7、设计具有下列指标的低通FIR 滤波器:通带截频0.2p ωπ=,阻带截频
0.3s ωπ=,通带波动0.25p R =dB ,最小的阻带衰减50s A =dB 。
解:
wp=0.2*pi;ws=0.3*pi; tr_width=ws-wp
%过渡带宽 M=ceil(6.6*pi/tr_width)+1 n=[0:1:M-1];
wc=(ws+wp)/2,%理想LPF 截止频率 hd=ideal_lp(wc,M);
%理想低通滤波器的单位冲激响应
w_ham=(hamming(M))'; h=hd.*w_ham; %截取得到实际的单位脉冲响应
[db,mag,pha,grd,w]=freqz_m(h,[1]); %计算实际滤波器的幅度响应 delta_w=2*pi/1000;
Rp=-(min(db(1:1:wp/delta_w+1))) %实际通带波动
As=-round(max(db(ws/delta_w+1:1:501))) %最小阻带衰减
subplot(1,1,1)
subplot(2,2,1);stem(n,hd);
axis([0 M-1 -0.1 0.3]);xlabel('n'); ylabel('hd(n)');
subplot(2,2,2);stem(n,w_ham);
title('哈明窗')
axis([0 M-1 0 1.1]);
xlabel('n');
ylabel('w(n)');
subplot(2,2,3);stem(n,h);
title('实际脉冲响应')
axis([0 M-1 -0.1 0.3]);xlabel('n'); ylabel('h(n)');
subplot(2,2,4);plot(w/pi,db);
title('幅度响应(单位:dB)');grid
axis([0 1 -100 10]);
xlabel('频率(单位:pi)');
ylabel('分贝数');
8、设带通滤波器的指标为:低阻带:
1
0.2,60
s s
A dB
ωπ
==;低通带:
1
0.35,1
p p
R dB
ωπ
==;高通带:
2
0.65,1 p p
R dB
ωπ
==;高阻带:
2
0.8,60
s s
A dB
ωπ
==。解:
clear all;
wpl=0.35*pi;
wph=0.65*pi;
wsl=0.2*pi;
wsh=0.8*pi;
tr_width=min((wpl-wsl),(wsh-wph)); %过渡带宽度
N=ceil(11*pi/tr_width)+1; %滤波器长度
n=0:1:N-1;
wcl=(wpl+wsl)/2;
%理想带通滤波器的下截止频率
wch=(wsh+wph)/2;
%理想带通滤波器的上截止频率
hd=ideal_bp(wcl,wch,N); %理想带通滤波器的单位冲激响应
w_bman=(blackman(N))'; %布莱克曼窗
h=hd.*w_bman;
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2(h,[1]);
%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
rp=-(min(db(wpl/delta_w+1:1:wph/delta_w+1))); %实际通带纹波
as=-round(max(db(wsh/delta_w+1:1:501))); %实际阻带纹波
subplot(2,2,1);
stem(n,hd);
title('理想单位脉冲响应');
subplot(2,2,2);
stem(n,w_bman);
title('布莱克曼窗');
subplot(2,2,3);
stem(n,h);
title('实际单位脉冲响应');
subplot(2,2,4);
plot(w/pi,db);
title('幅度响应');
图像如下:
9、用频率设计法设计具有下列指标的低通FIR 滤波器:通带截频0.2p ωπ=,阻带截频0.3s ωπ=,通带波动0.25p R =dB ,最小的阻带衰减50s A =dB 。
解:
clear all;
M=20;
alpha=(M-1)/2; l=0:M-1;
wl=(2*pi/M)*l;
Hrs=[1,1,1,zeros(1,15),1,1]; Hdr=[1,1,0,0];
wdl=[0,0.25,0.25,1];
k1=0:floor((M-1)/2);
k2=floor((M-1)/2)+1:M-1; angH=[-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)]; H=Hrs.*exp(j*angH); h=real(ifft(H,M));
[db,mag,pha,w]=freqz_m2(h,l);
[Hr,w,a,L]=Hr_Type2(h);
subplot(2,2,1);
plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);
title('频率样本:M=20');
xlabel('频率(单位:pi )');
ylabel('Hr(k)'); subplot(2,2,2); stem(l,h);
title('脉冲响应');
xlabel('n'); ylabel('h(n)'); subplot(2,2,3);
plot(w/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o'); title('振幅响应');
xlabel('频率(单位:pi )'); ylabel('Hr(w)'); subplot(2,2,4); plot(w/pi,db); title('幅度响应');
xlabel('频率(单位:pi )'); ylabel('分贝数');
图像如下:
10、
用最优设计法设计具有下列指标的低通FIR 滤波器:通带截频0.2p ωπ=,
阻带截频0.3s ωπ=,通带波动0.25p R =dB ,最小的阻带衰减50s A =dB 。
解:
clear all;
M=40;
alpha=(M-1)/2;
l=0:M-1;
wl=(2*pi/M)*l;
Hrs=[ones(1,5),0.5,zeros(1,29),0.5,ones(1,4)]; Hdr=[1,1,0,0];
wdl=[0,0.25,0.25,1];
k1=0:floor((M-1)/2);
k2=floor((M-1)/2)+1:M-1;
angH=[-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)]; H=Hrs.*exp(j*angH);
h=real(ifft(H,M));
[db,mag,pha,w]=freqz_m2(h,l);
[Hr,w,a,L]=Hr_Type2(h);
subplot(2,2,1);
plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);
title('频率样本:M=40,T=0.5');
xlabel('频率(单位:pi)');
ylabel('Hr(k)');
subplot(2,2,2);
stem(l,h);
title('脉冲响应');
xlabel('n');
ylabel('h(n)');
subplot(2,2,3);
plot(w/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');
title('振幅响应');
xlabel('频率(单位:pi)');
ylabel('Hr(w)');
图像如下:
11、
基于Butterworth 模拟滤波器原形使用冲激不变转换设计低通数字滤波器,其中参数指标为:通带截频0.2p πΩ=,通带波动值1p R =dB ,阻带截频Ωs=0.3π,阻带衰减值15s A =dB 。
解:
clear all
wp=0.2*pi;
%数字通带频率(弧度)
ws=0.3*pi;
%数字阻带频率
Rp=1;
%通带波纹(db)
As=15;
%阻带衰减(db)
T=1;
%性能指标
OmgP=wp/T;
%原型通带频率
OmgS=ws/T;
%原型阻带频率
[N,OmgC]=buttord(OmgP,OmgS,Rp,As,'s');
%选取模拟滤波器的阶数
[cs,ds]=butter(N,OmgC,'s');
%设计出所需的模拟低通滤波器
[b,a]=impinvar(cs,ds,T);
%应用脉冲响应不变法进行转换
%求得相对、绝对频响及相位、群迟延响应[db,mag,pha,grd,w]=freqz_m(b,a);
%下面绘出各条曲线
subplot(2,2,1);
plot(w/pi,mag);
title('幅频特性');
xlabel('w(/pi)');
ylabel('|H(jw)|');
subplot(2,2,2);
plot(w/pi,db);
title('幅频特性(dB)');
xlabel('w(/pi)');
ylabel('dB');
subplot(2,2,3);
plot(w/pi,pha/pi);
title('相频特性');
xlabel('w(/pi)');
ylabel('pha(/pi)');
subplot(2,2,4);
plot(w/pi,grd);
title('群延迟');
xlabel('w(/pi)');
ylabel('样本');
图像如下:
12、
基于Butterworth 模拟滤波器原形应用双线性变换设计低通数字滤波器,其中参数指标为:通带截频0.2p πΩ=,通带波动值1p R =dB ,阻带截频
0.2s πΩ=,阻带衰减值15s A =dB 。
解:
clear all
wp=0.2*pi; %数字通带频率(弧度) ws=0.3*pi; %数字阻带频率
Rp=1; %通带波纹(db)
As=15; %阻带衰减(db)
T=1;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
%预修正原型阻带频率
[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s');
%选取模拟滤波器的阶数
[cs,ds]=butter(N,OmegaC,'s');
%设计出所需的模拟低通滤波器
[b,a]=bilinear(cs,ds,T);
%应用双线性变换进行转换
[db,mag,pha,grd,w]=freqz_m(b,a);
subplot(2,2,1);
plot(w/pi,mag);
title('幅频特性');
xlabel('w(/pi)');
ylabel('|H(jw)|');
subplot(2,2,2);
plot(w/pi,db);
title('幅频特性(dB)');
xlabel('w(/pi)');
ylabel('dB');
subplot(2,2,3);
plot(w/pi,pha/pi);
title('相频特性');
xlabel('w(/pi)');
ylabel('pha(/pi)');
subplot(2,2,4);
plot(w/pi,grd);
title('群延迟');
xlabel('w(/pi)');
ylabel('样本');
图像如下: