王实验一 离散信号的matlab 实现
一、实验目的
1、 熟悉matlab 软件，学会matlab 语言的编写
2、 使用matlab 软件产生一些常见的离散信号
3、 掌握用matlab 软件作信号的相关分析
二、实验环境
计算机操作系统、matlab 软件
三、实验内容
1、用matlab 程序产生下列离散信号或连续信号，并画出其波形。
a 单位抽样序列的产生)(n δ
参考程序：N=100;
x=zeros(1,N); 产生一个1行N 列值全为0的矩阵，如看成数组x （1）-x （100）
都为0
x(1)=1;
n=0:N-1;
stem(n,x);
产生序列)20(-n δ
参考程序：N=100;
x=zeros(1,N);
k=20;
x(k+1)=1;
xn=0:N-1;
stem(xn,x);
b 单位阶跃序列的产生)(n u
参考程序：N=32;
x=ones(1,N);产生一个1行N 列值全为1的矩阵
n=0:N-1;
stem(n,x);
产生序列)20(-n u
参考程序：N=32;
k=20;
x1=zeros(1,k);
x2=ones(1,N-k);
x=[x1,x2];
xn=0:N-1;
stem(xn,x);
c 模拟信号)8cos(5)4sin(2)(t t t x ππ+=，以t=0.01n （n=0：N-1）进行采样后的离散信号。 参考程序： N=128；n=[0:N-1];
t=0.01*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
figure(1);
subplot(211);
stem(t,x);
subplot(212);
stem(n,x);
d产生一个sinc(t)=sint/t抽样函数
参考程序：n=200;
step=4*pi/n;
t=-2*pi:step:2*pi;
y=sinc(t);
plot(t,y,t,zeros(size(t)));%同时画出y(t）和横轴
grid on;
plot(t,y,t,zeros(size(t)),zeros(size(y)),y);%同时画出y(t）和横轴、纵轴e方波信号square(t) square(t,duty) 产生周期是2pi，幅度为正负1的方波，duty占空比,高电平跟整个周期的比值
参考程序：t=0:0.01:2*pi;
y=square(t,3);
plot(t,y);
试产生一个周期为1，高低电平分别为半个周期的方波信号
2、相关分析去除噪声
x(n)=sin(2*pi*n)+u(n) 噪声为高斯分布白噪声，使用相关分析去除噪声，噪声1功率为1，噪声2功率为0.1
%rxy=xcorr(x,y);
%rx=xcorr(x,Mlag,'flag') Mlag表示rx的单边长度，总长度为2Mlag+1, flag---'biased' rx(m)/N --unbiased rx(m)/(N-abs(m))
参考程序：N=500;
p1=1;
p2=0.1;
f=1/8;
Mlag=60;
u=randn(1,N);
u2=u*sqrt(p2);
n=[0:N-1];
s=sin(2*pi*f*n);
x1=u(1:N)+s;
rx1=xcorr(x1,Mlag,'biased');
subplot(211);
plot(-Mlag:Mlag,rx1);
x2=u2(1:N)+s;
rx2=xcorr(x2,Mlag,'biased');
subplot(212);
plot(-Mlag:Mlag,rx2);
3.书后习题1.15
实验二、离散信号的傅里叶变换
一、实验目的
1、进一步熟悉matlab软件的使用，熟悉matlab的编程语言
2、用matlab语言编写程序进行离散信号的傅里叶分析
二、实验原理
设离散序列)
(n
x，长度为N，其DTFT定义为：
∑∞
-∞=-
=
n jwn
j e
n
x
e
X)
(
)
(ω
在实际计算中无法取到无限长序列，通常通过无限长序列加窗作有限长序列的DTFT。
三、实验内容
试求序列)6/
cos(
)
(π
n
n
x=，其中n=0，1，……，N-1，N=12的DTFT。并画出其幅频特性曲线和相频特性曲线。
若N=24,36,120,其幅频特性曲线和相频特性曲线如何变化，为什么？
参考程序：N=12;
n=[0:1:N-1];
xn=cos(n*pi/6);
w=0:0.01:2*pi;
M=length(w);
xejw=zeros(1,M);
for i=0:1:N-1
xejw=xejw+cos(i*pi/6)*exp(-j*w*i);
end;
xr=abs(xejw);
xphase=angle(xejw);
subplot(211);
plot(w,xr);
subplot(212);
plot(w,xphase);
实验三、FFT 频谱分析及应用
一.实验目的
1、通过实验，加深对FFT 的理解，熟悉FFT 子程序
2、熟悉应用FFT 对典型信号进行频谱分析的方法
二、实验原理与方法
在各种信号序列中，有限长序列占有重要地位，对有限长序列，可以利用离散傅里叶变换（DFT ）进行分析。DFT 不但可以很好地反映序列的频谱特性，而且易于用快速算法（FFT ）在计算机上实现。
设离散序列)(n x ，长度为N ，其DFT 定义为：
kn N j kn N N n kn N e W W n x k X π
210,)()(--===∑
k=0，1，2……,N-1
有限长序列的DFT 是其Z 变换在单位圆上的等距采样，因此可以用于序列的谱分析。FFT 是DFT 的一种快速算法，它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。常用的FFT 是以2为基数的，其长度L N 2=。它的效率高，程序简单，使用方便。当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT ，可以用末尾补零的方法，使其长度延长至2的整数次方。
在MA TLAB 中，可以用函数FFT 来实现。
格式：y=FFT(x)
y=FFT(x,n)
说明：y=FFT(x)为利用FFT 算法计算矢量的离散傅里叶变换，当x 为矩阵时，y 为矩阵x 每一列的FFT 。当x 的长度为2的整数次方时，则FFT 函数采用基2的FFT 算法，否则采用稍慢的混合基算法。
y=FFT(x,n)采用n 点FFT 。当x 的长度小于n 时，FFT 函数在x 的尾部补零，以构成n 点数据；当x 的长度大于n 时，FFT 会截断序列x 。
三、实验内容
1、模拟信号)8cos(5)4sin(2)(t t t x ππ+=，以)1:0(01.0-==N n n t 进行采样。 求（1）N=40点FFT 的幅度频谱，从图中，能否观察出信号的2个频率分量？
（2）提高采样点数，如N=128,再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT 频谱分析结果与理论上是否一致？
2、研究高密度频谱与高分辨率频谱
设有连续信号
)10*9*2cos()10*7*2cos()10*5.6*2cos()(3
33t t t t x a πππ++= 以采样频率kHz f s 32=对信号)(t x a 采样，分析下列三种情况的幅频特性。
（1） 采样数据长度N=16点，做N=16的FFT ，并画出幅频特性
（2） 采集数据长度N=16点，补零到256点，做256点的FFT ，并画出幅频特
性
（3） 采样数据长度N=256点，做N=256的FFT ，并画出幅频特性。
参考程序：
N=40;
n=[0:N-1];
t=0.01*n;
x=2*sin(4*pi)+5*cos(8*pi*t);
k=[0:1:N/2];
w=2*pi/N*k;
X=fft(x,N);
magx=abs(X(1:1:N/2+1);
subplot(211);stem(n,x);title(‘signal x(n)’); subplot(212);plot(w/pi,magx);title(‘fft N=40’); xlabel(‘频率（单位：pi）’);
ylabel(‘| X|’);
grid;
取N=128时,程序参考以上程序.
实验内容2的程序参考以上程序.
实验四、IIR 数字滤波器的设计
一实验目的
1、掌握脉冲响应不变法和双线性变换法设计IIR 数字滤波器的具体方法和原理,熟悉双线性变换法和脉冲响应不变法设计低通,带通IIR 数字滤波器的计算机编程.
2、观察双线性变换法和脉冲响应不变法设计的数字滤波器的频域特性,了解双线性变换法和脉冲响应不变法的特点和区别.
3、熟悉Butterworth 滤波器和Chebyshev 滤波器的频率特性.
二实验原理与方法
本实验利用模拟滤波器设计IIR 数字滤波器,这是IIR 数字滤波器设计最常用的方法.利用模拟滤波器设计,需要将模拟域的)(s H a 转换为数字域的)(z H ,最常用的转换方法为脉冲响应不变法和双线性变换法.
1、脉冲响应不变法
用数字滤波器的单位脉冲响应序列)(n h 模仿模拟滤波器的冲激响应)(t h a ,让)(n h 正好等于)(t h a 的采样值,即:
)()(nT h n h a =
其中T 为采样间隔,如果以)(s H a 及)(z H 分别表示)(t h a 的拉氏变换及)(n h 的z 变换,则:
∑
∞-∞==+=m a e z m T j s H T z H sT )2(1
)(π
在MA TLAB 中，可用函数impinvar 实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：
[b,a]=impinvar(c,d,fs)
[b,a]=impinvar(c,d)
其中，c 、d 分别为模拟滤波器的分子和分母多项式系数向量，fs 为采样频率（Hz ），缺省值fs=1Hz ，b 、a 分别为数字滤波器分子和分母多项式系数向量。 例：已知3212)(2+++=
s s s s H a ，T=0.1，利用脉冲响应不变法求H(z) MA TLAB 程序：
c=[2 1];d=[1 2 3];T=0.1;
[b,a]=impinvar(c,d,1/T);
执行结果：b=0.2 -0.1881
a=1 -1.7916 0.8187
数字滤波器为：2118187.07916.111881.02.0)(---+--=z z z
z H
2、 双线性变换法
S 平面与z 平面之间满足下列映射关系：
11112--+-=z z
T s
s
T
s T
z -+=22
jw re z j s =Ω+=,σ
S 平面的虚轴映射在z 平面的单位圆上，s 平面的左半平面完全映射到z 平面的单位圆内。双线性变换不存在频率混叠问题。
在MA TLAB 中，提供了一个叫做bilinear 的函数实现这种映射，调用格式为：
[b,a]=bilinear(c,d,fs)
双线性变换是一种非线性变换，即22T tg
Ω=ω，这种非线性引起的幅频特性畸变可通过
预畸得到校正。
3、 一般设计步骤
（1） 给定技术指标转换为模拟低通原型设计性能指标。
（2） 估计满足性能指标的模拟低通原型阶数和截止频率。
利用MA TLAB 中buttord 、cheb1ord 、cheb2ord 、ellipord 等函数，调用格式如：
[n,wn]=buttord(wp,ws,rp,rs,’s ’)
（3）设计模拟低通原型
利用MA TLAB 中buttap 、cheb1ap 、cheb2ap 、ellipap 等函数，调用格式如：
[z,p,k]=buttap(n)
采用上述函数所得到原型滤波器的传递函数为零点、极点、增益形式，需要和函数
[c,d]=zp2tf(z,p,k)配合使用，以转化为多项式形式。
（4）由模拟低通原型经频率变换获得模拟低通、高通、带通或带阻滤波器。
利用MA TLAB 中lp2lp 、lp2hp 、lp2bp 、lp2bs 等函数，调用格式如：
[c1,d1]=lp2lp(c,d,wn)
（5）利用脉冲响应不变法或者双线性不变法，实现模拟滤波器到数字滤波器的映射。
说明：MA TLAB 信号处理工具箱还提供了模拟滤波器设计的完全工具函数：butter 、cheby1、cheby2、ellip 、besself 。用户只需一次调用就可自动完成以上设计步骤中的3-4步，调用格式如：
[c,d]=butter(n,wn,’ftype ’,’s ’)
三、实验内容
1、已知fp=0.3kHz ，rp=1.2db ，fs=0.2kHz,rs=20db,T=1ms,利用双线性变换法设计一个chebyshevI 型数字高通滤波器，观察通带损耗和阻带衰减是否满足要求。
按照要求写出matlab 程序，并附上所设计的滤波器传递函数H(z)及相应的幅频特性曲线。
2、要求fp=0.2kHz ，rp=1db ，fs=0.3kHz,rs=25db ，T=1ms ，分别用脉冲响应不变法和双线性变换法设计一个butterworth 数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。
程序1：
rp=1.5;rs=40;T=0.001;fp=0.3;fs=0.2;wp=2*pi*fp;ws=2*pi*fs;%wp,ws为数字频率wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);%预畸为模拟频率
[n,wn]=cheb1ord(wp1,ws1,rp,rs,’s’); %求阶数n，3db截止频率wn
[z,p,k]=cheb1ap(n,rp)；%模拟低通原型
[c,d]=zp2tf(z,p,k); %转化为多项式形式
[c1,d1]=lp2hp(c,d,wn); %模拟低通到模拟高通
[b,a]=bilinear(c1,d1,1/T); %双线性变换法设计出数字滤波器[db,mag,pha,grd,w]=freqz_m(b,a)；%作数字滤波器的幅频响应和相频响应subplot(211);
plot(w/pi,mag);
title('幅频特性');
xlabel('w(/pi)');
ylabel('|H(jw)|');
subplot(212);
plot(w/pi,pha/pi);
title('相频特性');
xlabel('w(/pi)');
ylabel('pha(/pi)');
freqz_m函数如下：
function[db,mag,pha,grd,w]=freqz_m(b,a)
%滤波器的幅值响应(相对、绝对)
%Usage:[db,mag,pha,grd,w]=freqz_m(b,a)
%w 采样频率；b系统函数H(z)的分子项(对FIR，b=h)
%a 系统函数H(z)的分母项(对FIR，a=1)
[H,w]=freqz(b,a); %500点的复频响应
mag=abs(H); %绝对幅值响应
db=20*log10(mag/max(mag)); %相对幅值响应
pha=angle(H); %相位响应
grd=grpdelay(b,a,w); %群延迟响应
程序2：
fp=0.2;fs=0.3;
wp=2*pi*fp;ws=2*pi*fs;%数字频率
rp=1;rs=25;T=0.001;
wp1=tan(wp/2);ws1=tan(ws/2);%预变形为模拟频率
[n,wn]=buttord(wp1,ws1,rp,rs,'s');%设计模拟滤波器阶数
[c,d]=butter(n,wn,'s');
[bz,az]=bilinear(c,d,1/2); %双线性变换法由s到z
[h,w]=freqz(bz,az,1024,1000);
plot(w,abs(h));grid on;
程序3：
fp=200;fs=300;
rp=1;rs=25;T=0.001;
wp=2*pi*fp;ws=2*pi*fs;% 模拟频率
[n,wn]=buttord(wp,ws,rp,rs,’s’);%设计模拟滤波器的阶数
[z,p,k]=buttap(n);
[bp,ap]=zp2tf(z,p,k);
[bs,as]=lp2lp(bp,ap,wp);
[bz,az]=impinvar(bs,as,1000);
%[bz,az]=butter(n,wn/pi);%设计数字滤波器
[db,mag,pha,grd,w]=freqz_m(b,a);%作出数字滤波器的频率特性图subplot(211);
plot(w/pi,mag);
title('幅频特性');
xlabel('w(/pi)');
ylabel('|H(jw)|');
subplot(212);
plot(w/pi,pha/pi);
title('相频特性');
xlabel('w(/pi)');
ylabel('pha(/pi)');
实验五、FIR 数字滤波器的设计
一、实验目的：
1、掌握用窗函数法设计FIR 滤波器的原理和方法，熟悉相应的计算机编程。
2、了解用频率采样法设计FIR 滤波器的原理和实现。
3、熟悉线性相位FIR 滤波器的幅频特性和相频特性。
4、了解不同窗函数对滤波器性能的影响。
二、实验原理与方法
1、窗函数法：
窗函数法设计线性相位FIR 滤波器步骤：
（1）确定数字滤波器的性能要求：临界频率{wk}，滤波器单位脉冲响应长度N 。
（2）根据性能要求，合理选择h(n)的奇、偶对称性，从而确定理想频率响应)(ωj d e H 的幅频特性和相位特性。
（3）求理想单位脉冲响应)(n h d 。在实际计算中，可对)(ωj d e H 按M 等间隔采样，并对其
求IDFT 得)(n h M ，用)(n h M 代替)(n h d 。
（4）选择合适的窗函数w(n),由)()()(n w n h n h d =求所设计的FIR 滤波器单位脉冲响应。
（5）求)(ωj e H ，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N ，重复上述设计过程，以得到满意的结果。
窗函数的傅氏变换)(ωj e W 的主瓣决定了)(ωj e H 的过渡带宽。)(ωj e W 的旁瓣大小和多少决定了)(ωj e H 在通带和阻带范围内的波动幅度。
2.频率采样法
频率采样法是从频域出发，将给定的)(ωj d e
H 加以等间隔采样： )()(2k H e H d k N j d ==π
ωω
然后以此)(k H d 作为实际数字滤波器的频率特性的采样值)(k H ，即令
1......2,1,0),()()(2-====N k k H e H k H d k N j d π
ωω
由H(k)通过IDFT 可得有限长序列h(n)
代入到z 变换中可得：
∑-=-----=1011)(1)(N k k N N z W k H N z
z H
用频率采样法设计线性相位FIR 滤波器的一般步骤为：
（1）由设计要求选择滤波器的种类
（2）根据线性相位的约束条件，确定k H 和k θ，进而得到)(k H
（3）将)(k H 代入)(ωj e H 内插公式得到所设计滤波器的频率响应。关于第3步，在MA TLAB 中可由函数
h=real(ifft(H,N))和[db,mag,pha,grd,w]=freqz_m(h,l)实现。
三、实验内容：
1、用窗函数法设计一线性相位FIR 低通滤波器，设计指标为：
db R db R s p s p 50,25.05.0,3.0====π
ωπω
（1）选择一个合适的窗函数，取N=15，确定脉冲响应，并给出所设计的滤波器的频率响应图，分析是否满足设计要求。
（2）若取N=45，重复这一设计，观察幅频和相位特性变化，分析长度N 变化的影响。
（3）保持N=45不变，改变窗函数，（如由hamming 窗变为blackman 窗），观察并记录窗函数对滤波器幅频特性的影响，比较两种窗的特点。
2、用kaiser 窗设计一个数字带通滤波器，设计指标为：
低阻带：db R s s 60,2.01==πω
低通带：db R p p 1,35.01==πω
高通带：db R p p 1,65.02==πω
高阻带：db R s s 60,8.02==πω
3、 用频率采样法设计一个低通滤波器：
db R db R s p s p 50,1,
35.0,2.0====πωπω
问：（1）采样点数N=33，过渡带设置1个采样点，H(k)=0.5，最小阻带衰减为多少，是否满足设计要求？
（2）采样点数N=34，过渡带设置2个采样点，1099.0)(,5925.0)(21==k H k H ,最小阻带
衰减为多少，是否满足设计要求？
参考程序1：
N=15;M=128;
b1=fir1(N,0.15,boxcar(N+1));
b2=fir1(N,0.15,hamming(N+1));
h1=freqz(b1,1,M);h2=freqz(b2,1,M);
t=0:N;
subplot(221);stem(t,b2,'.');grid;
f=0:0.5/M:0.5-0.5/M;
subplot(222);
plot(f,abs(h1),'b-',f,abs(h2),'g-');grid;
db1=20*log10(abs(h1)/1);db2=20*log10(abs(h2)/1);
subplot(224);
subplot(224);plot(f,db1,'b-',f,db2,'g-');
参考程序2：
N=15;M=128;
b2=fir1(N,[0.35,0.65],kaiser(N+1));
h2=freqz(b2,1,M);
t=0:N;
subplot(221);stem(t,b2,'.');grid;
f=0:0.5/M:0.5-0.5/M;
subplot(222);
plot(f*2,abs(h2),'g-');grid;
db2=20*log10(abs(h2)/1);
subplot(224);
subplot(224);plot(f*2,db2,'g-');
参考程序3：
%过渡带设置1个采样点
N=33;
alpha=(N-1)/2;
k=0:N-1;
wk=(2*pi/N)*k;
Hk=[ones(1,4),0.5,zeros(1,24),0.5,ones(1,3)];
angH=-alpha*(2*pi)/N*k;
H=Hk.*exp(j*angH);
h=real(ifft(H,N));
[hejw,w]=freqz(h,N);
subplot(211);
plot(w/pi,abs(hejw));
%过渡带设置2个采样点
N=34;
alpha=(N-1)/2;
k=0:N-1;
wk=(2*pi/N)*k;
Hk=[ones(1,4),0.5925,0.1099,zeros(1,23),-0.1099,-0.5925,-ones(1,3)]; angH=-alpha*(2*pi)/N*k;
H=Hk.*exp(j*angH);
h=real(ifft(H,N));
[hejw,w]=freqz(h,N);
subplot(212);
plot(w/pi,abs(hejw));