实验一 熟悉matlab
一、实验内容
练习matlab的使用，熟悉离散卷积运算，产生复指数序列
二、实验目的
1、掌握离散卷积计算机实现。
2、进一步对离散信号卷积算法的理解。
三、原理及算法概要
算法：把冲激响应h(n)与输入序列x(n)分别输入到程序中，然后调用离散卷积函数y=conv(x.,h)即可得到所要求的结果。
原理：离散卷积定义为
当序列为有限长时，则
四、程序与运行结果
% 实现复指数序列程序
n=0:30;
x=exp(.05+i*pi/6).^n;
xr=real(x);
xi=imag(x);
xm=abs(x);
xa=angle(x);
figure;
subplot(221);stem(n,xr);title('实部');
subplot(222);stem(n,xi);title('虚部');
subplot(223);stem(n,xm);title('模');
subplot(224);stem(n,xa);title('相角');
程序运行结果如图1所示。从图中可以看出，复指数序列的实部和虚部都是幅度按指数增长的正弦序列。
图1 复指数序列波形
x1=[1 1 1 1 ];nx1=0:3;% 实现卷积程序
h1=[1 0.8 0.64 0.8^3 0.8^4];nh1=0:4;
y1=conv(x1,h1);
subplot(2,3,1);stem(nx1,x1);title('序列x1');
xlabel('n');ylabel('x1(n)');
subplot(2,3,2);stem(nh1,h1);title('序列h1');
xlabel('n');ylabel('h1(n)');
subplot(2,3,3);stem(y1);title('序列y1');
xlabel('n');ylabel('y1(n)');
x2=[1 1 1 1];nx2=0:3;
nh2=0:1:20;
h2=(0.8).^nh2;
y2=conv(x2,h2);
subplot(2,3,4);stem(x2);title('序列x2');
xlabel('n');ylabel('x2(n)');
subplot(2,3,5);stem(h2);title('序列h2');
xlabel('n');ylabel('h2(n)');
subplot(2,3,6);stem(y2);title('序列y2');
xlabel('n');ylabel('y2(n)')
图2 卷积程序波形
五、结果分析
有限长序列的离散卷积计算结果与理论值一致，而存在无限长序列做卷积时，由于在程序处理时是用比较长有限长序列代替的，所以与理论值基本相同。
六、实习心得
主要进行的是有matlab实现离散卷积的计算，分为三个类型：冲激响应h(n)与输入序列x(n)都为有限长，一个序列为有限长一个序列为无限长和两个序列均为无限长。这部分实习内容不难，是原来做过的，主要是为了拾拣起原来所学的知识。第一种类型是最简单，也是最基本的，直接调用函数y=conv(x.,h)即可。在第二种和第三种类型的计算中遇到了一些困难，在输入序列时，由于存在无限长的序列，不知道该怎么输入，查过了相关的题目才弄明白。通过本专题的实习，让我对数字信号处理中的一些基础知识有了一些回顾，对原来所学过的知识也熟悉了一些。
实验二 FFT
一、实验内容
设有离散序列 x(n)= 0.51sin(2πn/32)+sin(4πn/32)-0.1sin(6πn/32);
采集数据长度N=32，分析32点的频谱，并画出幅频特性。
二、实验目的
1、了解DFT及FFT的性质和特点
2、利用FFT算法计算信号的频谱。
三、关键算法
算法：读入离散序列x(n)= 0.51sin(2πn/32)+sin(4πn/32)-0.1sin(6πn/32)，采集长度为N=32的数据，调用matlab中的函数fft（x）对其作离散傅里叶变换得到32点的频谱。
原理概要：当抽样数N=2M时，以下为算法蝶形图。
一般规律如下：
1、 当N=2M时，则要进行M次分解，即进行M级蝶形单元的计算
2、按自然顺序输入，输出是码位倒置。
3、每一级包含N/2个基本蝶形运算
4、第L级有2L-1个蝶群，蝶群间隔为N/2L-1
如果是Matlab实现的话，可用以下两种方法计算信号频谱
1、调用库函数为：fft（），直接计算X（k）
2、进行矩阵运算
四、程序及运行结果
n = 0:31;
x = 0.51*sin(2*pi*n/32)+sin(4*pi*n/32)-0.1*sin(6*pi*n/32);
y=fft(x);
y1=abs(y);
y2=angle(y);
subplot(2,1,1);stem(y1);
subplot(2,1,2);
stem(y2);
title('Hzhao');
typedef struct{float real;float imag;
}COMPLEX;
extern void fft(COMPLEX*x,int m);
extern void ifft(COMPLEX*x,int m);
void fft(COMPLEXS*x,int m)
{ static COMPLEX*w;
static int mstore=0;
static int n=1;
COMPLEX u,temp,tm;
COMPLEX*xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_reaL;
if(m! =mstore)
{if(mstore!=0)free(w);
mstore=m;
if(m==0)return;
n=1<<m;
le=n/2;
w=(COMPLEX*)colloc(le-1,sizeof(COMPLEX));
if(! w)
{printf("\nUnable to allocate complex w array\n");
exit(1); }
arg=4.0*atan(1.0)/le;
wrecur_real=w_real=cos(arg);
wrecur_imag=w_imag=-sin(arg);
xj=w;
for(j=1;j<le;j++)
{ xj->real=(float)wrecur_real;
xj->imag=(float)wrecur_imag;
xj++;
wtemp_real=wrecur_real*w_real_wrecur_imag*w_imag;
wrecur_imag=wrecur_real*w_imag+wrecur_imag*w_real;
wrecur_real=wtemp_real; }}
le=n;
windex=1;
for(l=0;1<m;1++)
{le=le/2;
for(i=0;i<n;i=i+2*le)
{xi=x+i;
xip=xi+le;
temp.real=(xi->real+xip->real);
temp,imag=(xi->imag+xip->imag);
xip->real=(xi->real-xip->real);
xip-imag=(xi->imag-xip->imag);
xi=temp; }
wptr=w+windex-1;
for(j=;j<le;j++)
{u=*wptr
for(i=j;i<n;i=i+2*le)
{ xi=x+i;
xip=xi+le;
temp.real=(xi->real+xip->real);
temp.imag=(xi->imag+xip->imag);
tm.real=xi->real-xip->real;
tm.imag=xi->imag-xip->imag;
xip->real=(tm.real*u.real-tm.imag*u.imag);
xip->imag=(tm.real*u,imag+tm.imag*u.real);
*xi=temp; }
wptr=wptr+windex; }
windex=2*windex; }
for(i=0;i<n;++i)
{ j=0
for(k=0;k<m;++k)
j=(j<<1)|(1&(i>>k));
if(i<j)
{ xi=x+i;
xj=x+j;
temp=*xj;
*xj=*xi;
*xi=temp; }} }
void ifft(COMPLEX*x,int m)
{static COMPLEX*w)
{static int mstore=0;
static int n=1;
COMPLEX u,temp,tm;
COMPLEX*xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
float scale;
double arg,w_real,w_imag,wrecur_real.wrecur_imag,wtemp_real;
五、结果分析
N 点DFT的频谱分辨率是2 π/N。一节指出可以通过补零观察到更多的频点，但是这并不意味着补零能够提高真正的频谱分辨率。这是因为x[n] 实际上是x(t) 采样的主值序列，而将x[n]补零得到的x'[n] 周期延拓之后与原来的序列并不相同，也不是x(t) 的采样。因此是不同离散信号的频谱。对于补零至M点的x'的DFT，只能说它的分辨率2 π/M仅具有计算上的意义，并不是真正的、物理意义上的频谱。频谱分辨率的提高只能通过提高采样频率实现。
六、专题实习心得
离散傅里叶变换是一种快速算法，由于有限长序列在其频域也可离散化为有限长序列，因此离散傅里叶变换在数字信号处理中是非常有用的。DFT是重要的变换，在分析有限长序列的有用工具、信号处理的理论上有重要意义、运算方法上起核心作用，谱分析、卷积、相关都可以通DFT在计算机上实现。DFT解决了两个问题：一是离散与量化，二是快速运算。通过编程实践体会到了时域、频域信号的对应关系，也对采样频率的含义有了深刻的认识，同时也加深了对采样信号频谱周期性的理解。
专题三 数字滤波器的设计
一、实验内容
1、设计一个Butterworth和Chebyshev数字低通滤波器，并学习模拟滤波器转换数字滤波器的两种方法，练习课本169页、181页、187页编程题
2、分析不同滤波器的特点和结果。
3、编程设计实现FIR滤波器。
二、实验目的
1.掌握不同IIR滤波器的性质、特点，和用窗函数法设计FIR滤波器的原理和方法。
2.通过实验学习如何设计各种常用的IIR滤波器和FIR滤波器，以便在实际工作中能根据具体情况使用IIR和FIR滤波器。
三、原理及算法概要
IIR滤波器算法：输入通带截止频率Wp，阻带截止频率Ws，通带衰减Rp，阻带衰减Rs，通过这些数值调用[N Wn]=buttord(Wp,Ws,Rp,Rs) 函数计算巴特沃斯数字滤波器的阶数N和截止频率wn，再根据阶数N通过函数[b,a]=butter(N,Wn)，即可得到所要的巴特沃斯滤波器。
IIR滤波器原理概要：
1、滤波器类型
1.1 Butterworth滤波器
Butterworth滤波器的特点是在通带内的频率特性是平坦的,并且随着频率的增加而衰减。Butterworth滤波器又是最简单的滤波器。
N阶低通Butterworth滤波器的幅度平方函数为:
1.2 ChebyshevⅠ型滤波器
ChebyshevⅠ型滤波器的在通带内的响应是等纹波的，而在阻带内是单调下降的，或在通带内是单调下降的，在阻带内是等纹波的特性。
其幅度平方函数为
其中VN是Chebyshev多项式。
1.3椭圆滤波器
椭圆滤波器在通带和阻带内都是等纹波振荡。
椭圆滤波器的特性函数为：
其中UN是N阶雅可比函数。
2、变换方法
2.1冲激响应不变法
冲激响应不变法的基本原理是从滤波器的冲激响应出发，对模拟滤波器冲激响应h(t)进行取样，所得到的离散序列h(nT)作为数字滤波器的单位取样响应。
H(z)是由H(s)通过下式的对应关系得到。
2.2双线性变换是在所得到满足性能指标要求的模拟滤波器的基础上,通过变换
,从而得到相应的数字滤波器。
FIR滤波器算法:通过其通带截止频率ωp与阻带截止频率ωs算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到hd（n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于40分贝，我选择最接近的汉宁窗，最后调用函数h=hd.*win 及freqz(h,1,512)得到实际汉宁窗的响应和实际滤波器的幅度响应。
FIR滤波器原理概要：
1、利用窗函数法设计FIR 滤波器
FIR滤波器的最大特点是其相位特性可以设计为严格的线性，而其幅值可以任意设置，这样输出波形就不会相位失真。
理想低通滤波器的单位取样响应hd（n）是无限长的，所以要用一个有限长的因果序列h（n）进行逼近，最有效的方法是截断hd（n），即用有限长的窗函数w（n）来截取 hd（n），表示为h（n）=hd（n）w（n）
为获得线性相位的FIR滤波器，h（n）必须满足中心对称条件，序列h（n）应有一定的延迟α，且α=（N-1）/2
频率响应逼近hd（ejw）的FIR滤波器，最简单的窗函数为矩形窗
W（n）= 1 n<(N-1)/2
0 n>(N-1)/2
加窗后的频谱
加窗后使实际频响偏离理想频响，影响主要有两个方面：
(1)通带和阻带之间存在过渡带，过渡带宽度取决于窗函数频响的主瓣宽度。
(2)通带和阻带区间有纹波，这是由窗函数的旁瓣引起的，旁瓣越多，纹波越多。增加窗函数的宽度N，其主瓣宽度减小，但不改变旁瓣的相对值。
为了改善滤波器的性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣衰减尽可能大，数量尽可能大，从而改善纹波状况，使实际频响H（ejω）更好地逼近理想频响Hd（ejω）
2.除了矩形窗外，一般还可以采用以下几种窗函数
Battlet窗 ：
汉宁窗：
海明窗
布来克曼窗
其中：WR是矩形窗的频谱
四、程序及运行结果
1.w1=2*400*tan(2*pi*90/(2*400));
w2=2*400*tan(2*pi*110/(2*400));
wr=2*400*tan(2*pi*120/(2*400));
[N,wn]=buttord([w1 w2],[100 wr],3,10,'s');
[B,A]=butter(N,wn,'s');
[num,den]=bilinear(B,A,400);
[h,w]=freqz(num,den);
f=w/pi*200;
plot(f,20*log10(abs(h)));
axis([40,160,-30,10]);
grid;
xlabel('频率/kHz')
ylabel('幅度/dB')
2.[B,A]=cheby1(3,0.7,2*pi*1000,'s');
[num1,den1]=impinvar(B,A,4000);
[h1,w]=freqz(num1,den1);
[B,A]=cheby1(3,0.7,2/0.00025,'s');
[num2,den2]=bilinear(B,A,4000);
[h2,w]=freqz(num2,den2);
f=w/pi*2000;
plot(f,abs(h1),'-',f,abs(h2),'-');
grid;
xlabel('频率/Hz')
ylabel('幅值/dB')
3.w1=2*400*tan(2*pi*90/(2*400));
w2=2*400*tan(2*pi*110/(2*400));
wr=2*400*tan(2*pi*120/(2*400));
[N,wn]=cheb1ord([w1 w2],[100 wr],3,10,'s');
Rp=0.7;
[B,A]=cheby1(N,Rp,wn,'s');
[num,den]=bilinear(B,A,400);
[h,w]=freqz(num,den);
f=w/pi*200;
plot(f,20*log10(abs(h)));
axis([40,160,-30,10]);
grid;
xlabel('频率/kHz')
ylabel('幅度/dB')
4. [B,A]=butter(3,2*pi*1000,'s');%3阶滤波器、3分贝截止频率2*pi*1000
[num1,den1]=impinvar(B,A,4000);%冲激响应不变法
[h1,w]=freqz(num1,den1);%freqz是绘制一个数字滤波器的频率响应的函数，返回值h1是通过 滤波器后的响应值向量，w是输入的频率向量
[B,A]=butter(3,2/0.00025,'s');
[num2,den2]=bilinear(B,A,4000);%双线性变换
[h2,w]=freqz(num2,den2);
f=w/pi*2000;
plot(f,abs(h1),'-',f,abs(h2),'-');
grid;
xlabel('频率/Hz')
ylabel('幅值/dB')
5.deltw=0.4*pi;
Wc=0.2*pi;
As=40;
N=ceil(8*pi/deltw)+1;
fprintf('N=%.0f\n',N);
win=hanning(N);
h=fir1(N-1,Wc/pi,win);
omega=linspace(0,pi,512);
mag=freqz(h,[1],omega);
magdb=20*log10(abs(mag));
plot(omega/pi,magdb);axis([0 1 -100 0]);
grid;xlabel('归一化频率');ylabel('幅度')
五、结果分析
Butterworth滤波器在通带内的频率特性是平坦的,并且随着频率的增加而衰减。正弦信号在经过IIR滤波器滤波后，高频信号被滤除，低频信号被保留了下来。
根据滤波器特性：通带截止频率处的衰减不大于3分贝；阻带衰减不小于40分贝。选用汉宁窗。FIR滤波器的最大特点是其相位特性可以设计为严格的线性，而其幅值可以任意设置，这样输出波形就不会相位失真。
六、专题实习心得
通过本专题的学习，熟悉和巩固了Butterworth数字低通滤波器的设计方法和原理，实现滤波器设计的有关经典算法，更重要的是熟练掌握使用MATLAB语言设计各种要求的数字滤波器。这一部分的内容相对来说是有些难度了，做起来花费的精力也多了一些，不过对数字信号处理内容的掌握上又加深了一层。
FIR有限冲击响应数字滤波器由于设计灵活，在滤波效果和响应时间之间能较好的折衷，因此在数字信号处理领域应用广泛。用窗函数法设计FIR数字滤波器在截断点依旧取模拟样本冲激响应的全值更利于理解和仿真。但是，窗函数法设计FIR数字滤波器的频率特性
与模拟样本频率特性不完全一样，且具有线性相移的特点。窗函数法设计FIR数字滤波器是傅立叶变换的典型运用。通过学习有益于加深对数字滤波器的理解，并提高使用傅立叶变换分析解决实际问题的能力。