数字信号处理实习报告.doc

中国地质大学（武汉）

数字信号处理上机实习

学生姓名：
班 级：071132
学 号：2013100
指导老师：王晓莉
题目一 离散卷积计算
一、 实验题目
设线性时不变（LTI）系统的冲激响应为h(n)，输入序列为x(n)
1、h(n)=(0.8)n，0≤n≤4； x(n)=u(n)-u(n-4)
2、h(n)=(0.8)nu(n), x(n)=u(n)-u(n-4)
3、h(n)=(0.8)nu(n), x(n)=u(n)
求以上三种情况下系统的输出y(n)，显示输入和输出波形。
二、 实验目的
1.理解和掌握离散卷积计算；
2.学习如何用Mtalab实现离散卷积计算。
三、算法设计
离散卷积定义为：

1、 ，，，

(a) 当 时，；
(b) 当时,(0.8)n ；
(c) 当时，(0.8)n；
(d) 当时，；

2、，，
(a) 当 时，；
(b) 当时,(0.8)n；
(c) 当时，(0.8)n；
(d) 当时，(0.8)n；
(e) 当时，；
3、，，
(a) 当 时，；
(b) 当时,(0.8)n；
(c) 当时，(0.8)n；
(d) 当时，；
四、 程序分析
所用到的函数：
（1）y=conv（x.,h）：卷积运算函数，计算;
（2）n1=0:4：n1取0~4；
（3）subplot(m,n,p):subplot()函数是将多个图画到一个平面上的工具。其中，m表示是图排成m行，n表示图排成n列，也就是整个figure中有n个图是排成一行的，一共m行，如果m=2就是表示2行图。p表示图所在的位置，p=1表示从左到右从上到下的第一个位置。
（4）title(‘content ’):title()函数的功能是为当前坐标系添加标题“content ”。
五、程序设计
n=0:4;
h=0.8.^n;
x=[1 1 1 1]
subplot(331);stem(x); title('x(n)');
y=conv(x,h);
subplot(332);stem(h); title('h(n)');
subplot(333);stem(y); title('y(n)');

n=0:40;
h=0.8.^n;
x=[1 1 1 1]
subplot(334);stem(x); title('x(n)');
y=conv(x,h);
subplot(335);stem(h); title('h(n)');
subplot(336);stem(y); title('y(n)');

n=0:40;
h=0.8.^n;
x=[zeros(1,0),ones(1,40)];
subplot(337);stem(x); title('x(n)');
y=conv(x,h);
subplot(338);stem(h); title('h(n)');
subplot(339);stem(y); title('y(n)');
六、运行结果

图中从左至右三列依次对应、及卷积结果

题目二 离散傅立叶变换
一、 实验题目
设有离散序列
分析下列三种情况下的幅频特性。
（1） 采集数据长度N=16，分析16点的频谱，并画出幅频特性。采集数据长度N=16，并补零到64点，分析其频谱，并画出幅频特性。
（2） 采集数据长度N=64，分析64点的频谱，并画出幅频特性。
观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。

二、 实验目的
1、理解掌握DFT及FFT算法；
2、利用FFT算法计算信号的频谱。

三、算法设计
当抽样数N=2M时，以下为蝶形算法图。

1、 当N=2M时，则要进行M次分解，即进行M级蝶形单元的计算；
2、按自然顺序输入，输出是码位倒置；
3、每一级包含N/2个基本蝶形运算；
4、第L级有2L-1个蝶群，蝶群间隔为N/2L-1；
四、程序分析
所用到的函数：
（1） fft（）：求x的一维傅里叶变换，计算X（k）；
使用方法Xk=fft(xn)，例如：
N=8;
n=0:N-1;
xn=[4 3 2 6 7 8 9 0];
Xk=fft(xn)
（2）stem（）:显示函数图像。
（3）xlabel():横坐标轴的名称；
（4）ylabel():纵坐标轴的名称；

五、程序设计
%x(n)=cos(0.48πn)+cos(0.52πn)
n=0:15; %产生序列x(n) 取16点
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(231);stem(n,x);title('采集数据长度N=16'); %显示x(n)
x1k=fft(x,16);%进行16点傅里叶变换
subplot(234);stem(0:15,abs(x1k),'.'); %显示X（k）
xlabel('n');ylabel('X1(k)'); title('16点傅立叶变换频谱');

x=[x zeros(1,48)]; %取16点，补零到64点
x2k=fft(x,64);%进行64点傅里叶变换
subplot(232);stem(x);title('采集N=16,补零到64');
subplot(235);stem(0:63,abs(x2k),'.');
xlabel('n');ylabel('X2(k)');
title('64点傅立叶变换频谱');

n=0:63; %产生序列x(n)，取64点
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(233);stem(n,x);title('采集数据长度N=64');
x3k=fft(x);%做64点傅里叶变换
subplot(236);stem(0:63,abs(x3k),'.');
xlabel('n');ylabel('X3(k)');
title('64点傅立叶变换频谱');

六、运行结果

七．结果分析
　 N点DFT的频谱分辨率是2 π/N。一节指出可以通过补零观察到更多的频点，但是这并不意味着补零能够提高真正的频谱分辨率。这是因为x[n] 实际上是x(t) 采样的主值序列，而将x[n]补零得到的x'[n] 周期延拓之后与原来的序列并不相同，也不是x(t) 的采样。因此是不同离散信号的频谱。对于补零至M点的x'的DFT，只能说它的分辨率2 π/M仅具有计算上的意义，并不是真正的、物理意义上的频谱。频谱分辨率的提高只能通过提高采样频率实现。第三幅图形取样点数较多，分辨率最高，而第一幅图形则显得较为稀疏，主要是因为取样点数太少。

题目三 IIR滤波器的设计
一、 实验题目
1、设计一个切比雪夫数字低通滤波器，设计指标如下：
通带截止频率：0.2π，幅度衰减不大于3分贝
阻带截止频率：0.3π，幅度衰减大于20分贝
2、分析不同滤波器的特点和结果。
3、编程设计实现IIR滤波器。
二、实验目的
1.理解和掌握不同IIR滤波器的性质、特点。
2.掌握用Matlab设计切比雪夫滤波器和椭圆滤波器的方法。
三、算法设计
1、滤波器类型
切比雪夫滤波器：切比雪夫滤波器，又名“车比雪夫滤波器”，是在通带或阻带上频率响应幅度等波纹波动的滤波器。
2、变换方法
（a）冲激响应不变法
冲激响应不变法的基本原理是从滤波器的冲激响应出发，对模拟滤波器冲激响应h(t)进行取样，所得到的离散序列h(nT)作为数字滤波器的单位取样响应。
H(z)是由H(s)通过下式的对应关系得到。

（b）双线性变换是在所得到满足性能指标要求的模拟滤波器的基础上,通过变换
,从而得到相应的数字滤波器。

四、 程序分析
所用到的函数：
（1）[N,wn]=buttord(Wp,Ws,ap,as,'s')：[阶数，截止频率]=buttord（通带截止频率，阻带截止频率，通带衰减，阻带衰减，模拟）；
（2）[num,den]=bilinear(b,a,fs)：双线性变换可选频率预畸变，在IIR数字滤波器设计中有把模拟的原型滤波器转变成相应的数字滤波器，即通过已知S变换中模拟滤波器的系数求数字滤波器的系数。其中有两种方法来进行这种变换，一种是激冲不变法，另一种是双线性Z变换法。bilinear 函数便是双线性Z变换法。函数中的b和a是原模拟滤波器S变换中的分子和分母系数，num和den是数字滤波器的系数。
（3）grid on：增加主要网格线为当前轴；
（4）[h,w]=freqz(num,den)：求频率响应。函数的输出：a.滤波器的频率响应H（N点） b.频率向量W(N点，且单位为弧度)，频率向量W是均匀分布在滤波器的上半区，即：0：pi，这些点上的频率响应都将通过此函数计算出来。；滤波器的系数：分子为num，分母为den。  

五、程序设计
clc;
Rs=20; Rp=3;
Wp1=0.2*pi;
Ws1=0.3*pi;
[N,Wn]=cheb1ord(Wp1,Ws1,Rp,Rs,'s');
[num1,den1]=cheby1(N,Rp,Wn,'s');
[num,den]=bilinear(num1,den1,0.5);
subplot(2,2,1);
zplane(num,den);
title('零极点图')
w=0:pi/256:pi;
h=freqz(num,den,w);
g=abs(h);
g1=angle(h);

subplot(223);
plot(w/pi,g);grid
axis([0 1 0 1]);
xlabel('频率'); ylabel('幅度'); title('幅频响应');
subplot(224);
plot(w/pi,g1);grid on
axis([0 1 -5 5]);
xlabel('频率'); ylabel('相位'); title('相频响应');

%测试
t=1:300;
I=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t);%设计正弦波
figure;subplot(311);
plot(I);title('原信号x(n)=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t)');
A=filter(num,den,I);%正弦波通过滤波器
subplot(312);
plot(A);title('滤波后y(n)');
六、 运行结果
（1） 零极点图、幅频响应、相频响应

（2）输入为I=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t)时的滤波结果

七、 结果分析
由幅频特性可知，椭圆滤波器在通带内的频率特性是平坦的,并且随着频率的增加而衰减。通带截止频率：0.2π，幅度衰减不大于3分贝，阻带截止频率：0.3π，幅度衰减大于20分贝。正弦信号在经过IIR滤波器滤波后，由傅里叶变换后的频谱图可看出高频信号被滤除，低频信号被保留了下来。

题目四 FIR滤波器的设计
一、实验内容
选取合适窗函数设计一个线性相位FIR低通滤波器，使它满足如下性能指标：
通带截止频率：ωp=0.66π,通带截止频率处的衰减不大于3分贝；
阻带截止频率：ωs=0.5π，阻带衰减不小于40分贝
二、实验目的
1、掌握用窗函数法设计FIR滤波器的原理和方法。
2、熟悉线性相位滤波器特性。
3、了解各种窗函数对滤波器特性的影响。
三、原理及算法概要
1、算法
通过其通带截止频率ωp与阻带截止频率ωs算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到hd（n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于40分贝，我选择最接近的汉宁窗，最后调用函数h=hd.*win 及freqz(h,1,512)得到实际汉宁窗的响应和实际滤波器的幅度响应。

2、原理概要
利用窗函数法设计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滤波器，最简单的窗函数为矩形窗：
1 n<(N-1)/2
W（n）=
0 n>(N-1)/2
加窗后的频谱 加窗后使实际频响偏离理想频响，影响主要有两个方面：
（1） 通带和阻带之间存在过渡带，过渡带宽度取决于窗函数频响的主瓣宽度。
（2） 通带和阻带区间有纹波，这是由窗函数的旁瓣引起的，旁瓣越多，纹波越多。
增加窗函数的宽度N，其主瓣宽度减小，但不改变旁瓣的相对值。为了改善滤波器的性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣衰减尽可能大，数量尽可能大，从而改善纹波状况，使实际频响H（ejω）更好地逼近理想频响Hd（ejω）。
除了矩形窗外，一般还可以采用以下几种窗函数
汉宁窗：
海明窗
布来克曼窗
四、 程序设计
clc;
wp=0.66*pi;ws=0.5*pi;wd=wp-ws;
%hanning窗
N_hann=ceil(8*pi/wd); wc=(wp+ws)/2;
if mod(N_hann,2)==0
N_hann=N_hann+1;
end
hann_win=hanning(N_hann);
b_hann=fir1(N_hann-1,wc/pi,'high',hann_win);
freq_axis=[0:pi/512:pi-pi/512];
freqn=[0:511]/512;
H_hann=freqz(b_hann,1,512);
figure(1);subplot(311);
plot(freqn,20*log10(abs(H_hann)));
xlabel('归一化频率w/pi');ylabel('幅度谱（dB）');title('采用hanning窗设计fir滤波器幅度谱');
subplot(312);
plot(freqn,angle(H_hann));
xlabel('归一化频率w/pi');ylabel('相位谱');title('采用hanning窗设计fir滤波器相位谱');
subplot(313);
stem(b_hann);
xlabel('n');ylabel('h(n)');title('单位脉冲响应序列');
%hamming窗
N_hamm=ceil(8*pi/wd);wc=(wp+ws)/2;
if mod(N_hamm,2)==0
N_hamm=N_hamm+1;
end
hamm_win=hamming(N_hamm);
b_hamm=fir1(N_hamm-1,wc/pi,'high',hamm_win);
H_hamm=freqz(b_hamm,1,512);
figure(2); subplot(311);
plot(freqn,20*log10(abs(H_hamm)));
xlabel('归一化频率w/pi');ylabel('幅度谱（dB）');title('采用hamming窗设计fir滤波器幅度谱');
subplot(312);
plot(freqn,angle(H_hamm));
xlabel('归一化频率w/pi');ylabel('相位谱');title('采用hamming窗设计fir滤波器相位谱');
subplot(313);
stem(b_hann);
xlabel('n');ylabel('h(n)');title('单位脉冲响应序列');
五． 程序运行结果
（1）采用hanning窗设计FIR滤波器

（2）采用hamming窗设计FIR滤波器

六． 结果分析
上图为选用汉宁窗和汉明窗的FIR 滤波器的相频特性、幅频特性及单位取样响应hd（n）。
由运行得出的幅频特性可知，满足设计指标通带截止频率：ωp=0.66π；阻带截止频率：ωs=0.5π。得出的窗主瓣宽度较宽，旁瓣宽度较窄，也合乎要求。

题目五
一、实验内容
录制一段电话的拨号音，并对录制的信号进行采样；画出采样后信号的时域波形和频谱图；必要的话对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，根据信号的频谱解析出拨号的数字。最后，用MATLAB设计一信号处理系统界面。

二、实验目的
1、掌握用滤波器处理语音信号的原理和方法。
2、掌握用matlab设计界面。
三、原理及算法概要
调用函数function pushbutton1_Callback(hObject, eventdata, handles)实现一个信号处理系统界面。选择button1时，用双线性变换法设计滤波器来对信号进行处理。读取语音信号，对语音信号进行f=8000的频率进行采样，调用函数y1=fft(x1,2048)对所采集的点做2048点FFT变换。先设计butterworth模拟滤波器，再用双线性变换法实现模拟滤波器到数字滤波器的转换。最后调用函数f1=filter(bz,az,x2)对加了噪声的语音信号进行滤波，得到滤波后的频谱图。
四． 程序

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
clc;
fs=8000; %语音信号采样频率为8000
x1=wavread('D:\matlab2011\work\shua.wav');
t=(0:length(x1)-1)/8000;
y1=fft(x1,2048); %对信号做2048点FFT变换
f=fs*(0:1023)/2048;
plot(handles.axes1,t,x1); %做原始信号的时域波形
grid on;axis tight;
title(handles.axes1,'原始语音信号');xlabel(handles.axes1,'time(s)');ylabel(handles.axes1,'幅度');
plot(handles.axes2,f,abs(y1(1:1024))); %做原始信号的FFT频谱
grid on;axis tight;
title(handles.axes2,'原始语音信号的FFT频谱');
xlabel(handles.axes2,'Hz');ylabel(handles.axes2,'幅度');%双线性变换法设计的巴特沃斯滤波器
A1=0.05;A2=0.10;
d=[A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t)]';
x2=x1+d;
wp=0.8*pi;ws=0.85*pi;
Rp=1;Rs=15;
Fs=8000;Ts=1/Fs;
wp1=2/Ts*tan(wp/2); %将模拟指标转换为数字指标
ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %选择滤波器最小阶数
[Z,P,K]=buttap(N); %创建butterworth模拟滤波器
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %用双线性法实现模拟到数字的转换
[H,W]=freqz(bz,az); %绘制频率响应曲线
plot(handles.axes3,W*Fs/(2*pi),abs(H));
grid on;axis tight;
xlabel(handles.axes3,'频率（Hz）');ylabel(handles.axes3,'频率响应');
title(handles.axes3,'Butterworth');
f1=filter(bz,az,x2);
plot(handles.axes4,t,x2); %画出滤波前的时域图
grid on;axis tight;title(handles.axes4,'加噪声信号滤波前的时域波形');
plot(handles.axes5,t,f1); %画出滤波后的时域图
grid on;axis tight;title(handles.axes5,'滤波后的时域波形');
y3=fft(f1,2048);y2=fft(x2,2048);
plot(handles.axes6,f,abs(y2(1:1024))); %画出滤波前的频谱图
grid on;axis tight;title(handles.axes6,'加噪声信号滤波前的频谱');
xlabel(handles.axes6,'Hz');ylabel(handles.axes6,'幅度');
plot(handles.axes7,f,abs(y3(1:1024))); %画出滤波后的频谱图
grid on;axis tight;
title(handles.axes7,'滤波后的频谱');
xlabel(handles.axes7,'Hz');ylabel(handles.axes7,'幅度');

四． 实验结果
1. 实验界面

2. 原始信号和加噪声信号的时域、频域波形

3. 巴特沃斯滤波器

4. 滤波后信号的时域、频域波形

五．结果分析
所以从语音的频谱可以看出，通过低通滤波器之后，滤波后的时域波形对滤波前的时域波形看起来相对光滑一些，滤波后的频谱图，滤除了添加的高频噪声信号，保留了低频信号，语音效果也改善了不少。通过最后的频谱图可以看出，与之前语音信号的频谱图几乎一样。而通过比较双线性变换法的滤波器可以看出，它们的滤波效果相似。