数字信号处理上机实验报告
姓名：马菲 班级：030812 学号：03081161
一．实验内容
a) 利用傅立叶级数展开的方法，自由生成所需的x(t)；
b) 通过选择不同的采样间隔T（分别选T>或<1/2fc），从x(t)获得相应的x(n)（作出x(n)图形）；
c) 对获得的不同x(n)分别作傅立叶变换，分析其频率响应特性（给出幅频与相频特性曲线）；
d) 利用巴特沃思、切比雪夫或椭圆滤波器设计数字滤波器（滤波特性自定），要求通过改变滤波器参数或特性（低通、高通、带通或带阻）设计至少两种数字滤波器，分析所设计滤波器（画出频率特性曲线），并对上述给出的不同x(n)分别进行滤波（画出滤波结果），然后加以讨论；
e) 利用窗函数设计法或频率采样法设计数字滤波器（滤波特性自定），要求通过改变滤波器参数或特性（低通、高通、带通或带阻等）设计至少两种数字滤波器，分析所设计滤波器（画出频率特性曲线），并对上述给出的不同x(n)分别进行滤波（画出滤波结果），然后加以讨论。
二．实验方法
通过MATLAB编程实现生成x(t)，以及对其进行上述变换等。
1. 生成x(t)
利用傅里叶级数展开的方法
傅里叶级数：f(t)=A0/2+∑cos(nwn+rn)
2 通过对N点采样，得到频率响应曲线
3 利用巴特沃斯设计低通数字滤波器和高通数字滤波器：MATLAB信号处理工具箱函数buttap,buttord和butter是巴特沃斯滤波器设计函数。
4 通过汉宁窗和哈明窗设计数字滤波器，调用工具箱函数firl实现窗函数法的设计。
三．部分代码及结果分析
1． 生成x（t），x（n）及频率响应特性
（1）通过傅里叶级数展开方法自定义了一个周期信号y=x（t）如下：
%周期信号生成
y=3*sin(2*pi*t*10)+0.5*sin(2*pi*t*15)+sin(2*pi*t*40)+sin(2*pi*t*150)
（2）选取采样间隔Ts=1400得到相应的序列x（n）：
%采样序列
Ts=1/1400;
n=0:279;
t=n*Ts;
Y=3*sin(2*pi*t*10)+0.5*sin(2*pi*t*15)+sin(2*pi*t*40)+sin(2*pi*t*150);
subplot(2,2,2),stem(n,Y);grid on;
title('N点采样序列');
xlabel('N');
ylabel('Y');
根据采样间隔取的不同，得到的序列点的密集程度就不一样
（3）通过对x（n）作傅里叶变换，得到其幅频特性曲线与相频特性曲线
这里以横坐标为采样频率：f=(0:length(Y)-1)'/(length(Y)*Ts);
通过傅里叶变换的函数，得到
幅频特性曲线：F=abs(fft(Y));相频特性曲线：A=unwrap(angle(fft(Y)));
2. 通过巴特沃斯设计低通数字滤波器和高通数字滤波器
（1）分析：[N,wc]=buttord(wp,ws,rp,rs); 计算阶数N和3dB截止频率。调用参数wp和ws分别为数字滤波器的通带边界频率和阻带边界频率的归一化值。当ws≤wp时，为高通滤波器，当ws≥wp时，为低通滤波器。
低通滤波器：
wp=8*pi/140/pi;ws=pi/7/pi;rp=2;rs=100;%数字滤波器指标参数
[N,wc]=buttord(wp,ws,rp,rs); %//计算阶数N和3dB截止频率。
[B,A]=butter(N,wc); %//求数字滤波器系统函数H（z）,B,A分别为分子和分母多项式的系数向量。
[H,W]=freqz(B,A,1000); %//计算数字滤波器的频率响应。1000为采样点数。得到的H。W都为1000个。
m=abs(H);p=unwrap(angle(H));
高通滤波器：
wp=pi/7/pi;ws=8*pi/140/pi;rp=2;rs=100;%//数字滤波器指标参数
[N,wc]=buttord(wp,ws,rp,rs); %//计算阶数N和3dB截止频率。
[B,A]=butter(N,wc,'high'); %//求数字滤波器系统函数H（z）,B,A分别为分子和分母多项式的系数向量。
[H,W]=freqz(B,A,1000); %//计算数字滤波器的频率响应。1000为采样点数。得到的H。W都为1000个。
m=abs(H);p=unwrap(angle(H));
（2）分析：
对信号滤波，选取采样间隔，对信号进行N点采样，并对信号Y进行低通滤波
低通：
f=(0:length(Y)-1)'/(length(Y)*Ts);
ft_Y=filter(B,A,Y); %对信号Y进行低通滤波
subplot(2,1,2),plot(n,ft_Y);%滤波后的信号
title('滤波后序列');
xlabel('N');
ylabel('ft-Y');
高通：
subplot(2,1,2),plot(n,ft_Y);%滤波后的信号
title('滤波后序列');
xlabel('N');
ylabel('ft-Y');
（2）描绘出其幅频特性曲线和相频特性曲线
低通：
subplot(2,1,1),
stem(f,abs(fft(ft_Y)));
title ('信号Y滤波后的幅度特性曲线');
xlabel('f(Hz)');
ylabel('Y(e^{j\omega})');
A=unwrap(angle(fft(ft_Y)));
subplot(2,1,2),stem(f,A);grid on;
title('信号Y滤波后的相频特性曲线');
xlabel('f(Hz)');
ylabel('\phi(e^{j\omega})');
高通：
subplot(2,1,1),
stem(f,abs(fft(ft_Y)));title ('对信号Y进行低通滤波');
title ('信号Y滤波后的幅度特性曲线');
xlabel('f(Hz)');
ylabel('Y(e^{j\omega})');
A=unwrap(angle(fft(ft_Y)));
subplot(2,1,2),stem(f,A);grid on;
title('信号Y滤波后的相频特性曲线');
xlabel('f(Hz)');
ylabel('\phi(e^{j\omega})');
3. 利用汉宁窗和海明窗设计数字滤波器
（1） 汉宁窗
分析：计算wp和ws；然后代入公式Bt=wp-ws来计算过渡带宽度；已知汉宁窗的精确值是6.2，计算所需h（n）的长度N0；计算理想高通滤波器通带截止频率，归一化；调用firl计算高通滤波器的h（n）。（幅频相频曲线同前分析）
wp=pi/7;ws=8*pi/140;
Bt=wp-ws;%过渡宽度
N0=ceil(6.2*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc,'high',hanning(N));
（2） 哈明窗
分析：哈明窗的过渡带宽度精确值是6.6。
wp=pi/7;ws=8*pi/140;
Bt=wp-ws;
N0=ceil(6.6*pi/Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc,'high',hamming(N));
四．总结
通过这次实验，我又一次熟悉了MATLAB，并通过MATLAB编程实现了信号的生成，以及数字滤波器的设计。也通过跟同学之间的讨论学习让我发现了以前不扎实的地方。在今后的学习中，我一定会更加努力，进一步提升自己的能力。
第二篇：数字信号实验报告
数字信号处理导论实验报告
姓名：张伦裕
学号：200905090205
实验一 信号、系统及系统响应
当系统的输入输出差分方程为：
Y(n)-0.8y(n-1)-0.5y(n-2)=0.7x(n)+0.3x（n-1）
通过 MATLAB 编程实现并考虑如下问题：
（1）当系统的输入为单位冲激函数时，分别利用filter 函数和impz 函数得到的系统单位冲激响应曲线。
（2）当系统的输入为单位阶跃函数时，分别利用filter 函数和impz 函数得到的系统单位阶跃响应曲线。
（3）对于不同的输入，系统的输出有什么变化，试讨论之。
一、实验原理
系统的零状态响应就是在系统初始状态为零的条件下微分方程的解。
y=filter（b，a，f） 其中b，a分别为系统差分方程左右端系数向量，f为数输入向量。
impz函数可直接求出离散系统的单位冲激响应，为：impz（b，a）。
conv函数可用于计算离散序列的线性卷积。为：y=conv（x，h）。
可用impz函数和conv函数计算对任意输入的离散系统的响应，即由impz求出单位冲激响应，再用单位冲激应与输入序列卷积。
二、程序
n=0:20;
A=[1,-0.8,-0.5];
B=[0.7,0.3];
y1=filter(B,A,u1);
[y2]=impz(B,A,21);
Y1=filter(B,A,u2);
[h]=impz(B,A,21);
Y22=conv(h,u2);
Y2=Y22(1:21);
subplot(221);stem(n,y1);
title('filter function');
grid on;
subplot(222);stem(n,y2);
title('impz function');
grid on;
subplot(223);stem(n,Y1);
title('filter function');
grid on;
subplot(224);stem(n,Y2);
title('impz function');
grid on;
三、结果及分析：
实验二 使用FFT作频谱分析
编程实现下列问题并讨论所得到的结果。
（1） 使用FFT 对MATLAB 中的数据noissin.dat 进行谱分析。
（2） 使用FFT 对语音数据noisyspeech.wav 进行谱分析。
一、实验原理
（1）离散傅里叶变换(DFT)公式为：
X(k)=∑x(n)W^nk;
x(n)=∑X(k)W^-nk;
其中w=e^(2∏nk/N),N为离散序列的长度。
（2）快速傅里叶变换(FFT)是利用w因子的取值特点，减少DFT的复数乘法的次数。其中一种是时间抽取基2算法，它将时间按奇偶逐级分开，直到两点的DFT。MATLAB提供了fft函数可用于计算FFT，器调用形式为;X=fft（x）或X=fft（x，N），N为2的整数次幂，若x的长度小于N，则补零，若超过则舍去N后的数据。
（3）函数形式[y,fs,bits]=wavread('Blip',[N1 N2]);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1 N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。
函数形式s=noissin（n1，n2）用于读取MALAB的噪声信号。
二、程序
xx=wavread('noisyspeech.wav');
fs=100;N=128;
x=xx(1:N);
n=1:N;
X=abs(fft(x,N));
subplot(221);plot(n,x);
xlabel('n');ylabel('x(n)');
grid on;
subplot(222);plot(n,X);
xlabel('k');ylabel('|X(k)|');
grid on;
load noissin;
s=noissin(1:20);
S=fft(s);
subplot(223);plot(abs(s));
xlabel('n');ylabel('|s(n)|');
grid on;
subplot(224);plot(abs(S));
xlabel('k');ylabel('|S(k)|');
grid on;
三、结果与分析
实验三 使用双线性Z变换设计IIR滤波器
编程实现下列问题并讨论所得到的结果。
（1） 使用FFT 对MATLAB 中的数据noissin.dat 进行谱分析。
（2） 使用FFT 对语音数据noisyspeech.wav 进行谱分析。
一、实验原理
（1）设计滤波器就是要设计一个系统是其能让一定频率的波段通过或滤去，对IIR滤波器，器转移函数是：H(Z)=(∑bz^(-r))/(1+∑az^(-k))。
（2）设计的一般原则：若使滤波器拒绝某个频率，应在单位园上相应的频率处设置一个零点，反之则设置一个极点。
(3)低通数字IIR滤波器设计步骤：
给出数字低通滤波器的技术要求→映射为模拟低通的技术要求→归一化为模拟低通滤波器的技术要求→设计出G(P)→G(S)→映射到数字滤波器的转移函数G(Z)。
(3)双线性Z变换，即S平面到Z平面的映射关系：S=2（Z-1）/Ts（Z+1）。
二、程序
fp=100;fs=300;Fs=1000;
rp=3;rs=20;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
Fs=Fs/Fs;
wap=tan(wp/2);was=tan(ws/2);
[n,wn]=buttord(wap,was,rp,rs,'s');
[z,p,k]=buttap(n);
[bp,ap]=zp2tf(z,p,k);
[bs,as]=lp2lp(bp,ap,wap);
[bz,az]=bilinear(bs,as,Fs/2);
[h,w]=freqz(bz,az,256,Fs*1000);
plot(w,abs(h));
xlabel('w');ylabel('|h(w)|');
grid on;
bp
ap
bs
as
bz
az
三，结果及分析
bp =
0 0 1
ap =
1.0000 1.4142 1.0000
bs =
0.1056
as =
1.0000 0.4595 0.1056
bz =
0.0675 0.1349 0.0675
az =
1.0000 -1.1430 0.4128
实验四使用窗函数发设计FIR滤波器
题目：根据下列指标，设计一个FIR数字低通滤波器：ω=2.0∏，ω=4.0∏
Ap=0.25dB，As=50dB。
（1） 分别考查矩形窗、三角窗、Hanning窗、Hamming窗并分析这些不同的窗函数对滤波器的幅度响应有什么影响。
（2） 选择一个合适的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。
一、实验原理
（1）窗函数法：将理想低通数字滤波器Hd（z）的单位抽样响应hd（n）截短并位移，可得到对理想滤波器的近似的转移函数H（z）。此过程即为在hd(n)上施加矩形窗口。
（2）FIR DF设计的窗含函数法仅能给出通带截止频率wp，其他几个参数是靠h（n）的长度及所使用的窗含数的性能来决定的。常用的窗函数有：矩形窗、三角窗、汉宁窗、汉明窗。
MATLAB窗函数法设计FIR数字低通滤波器用函数firl其调用形式为：b=firl（N，Wn）。式中N为滤波器的阶次，Wn是通带截止频率，b是设计好的滤波器系数h（n）。
（3）freqz函数，用来求系统的频率响应，其调用形式为：[H,w]=freqz(b,a,N,’whole’,Fs)。
二、程序
N=128;
Wp=0.2*pi;Ws=0.4*pi;
Wc=(Ws+Wp)/2;
b=fir1(N,Wc,hamming(N+1))
[h,w]=freqz(b,1,N,1);
plot(w,abs(h));
xlabel('w');ylabel('|h(w)|');
grid on;
三、结果及分析
自编FFT程序：
（用MATLAB实现）
function Xk= ffts(xn)
N=length(xn);
M=log2(N);
n1=fliplr(dec2bin(0:N-1));
n2=bin2dec(n1);
Xk=ones(1,N);
for h=1:N
Xk(h)=xn(n2(h)+1);
end
for i=1:M
s=2^i;
g=N/s;
b=s/2;
for k=0:g-1;
for m=0:b-1;
O=-2*pi*m*g/N;
w=cos(O)+1i*sin(O);
n=k*s+m;
y=w*Xk(n+b+1);
G=Xk(n+1);
H=Xk(n+1);
Xk(n+1)=G+y;
Xk(n+b+1)=H-y;
end
end
end
例：
xn=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15];
y1=fft(xn);
y2=ffts(xn);
subplot(121);stem(abs(y1));
subplot(122);stem(abs(y2));
结果