数字信号处理实验报告二与三

实验二 用FFT 进行谱分析
一.实验目的：
1 进一步加深DFT 算法原理和基本性质的理解（因为FFT 只是DFT 的一种快速算法，所以FFT 的运算结果必然满足DFT 的基本性质）。

熟悉FFT 程序结构及编程方法。

2 熟悉应用FFT 对确定信号进行谱分析方法，熟悉FFT 算法原理和FFT 子程序的应用。

3 学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应有FFT 。

二.实验内容：
（1）用matlab 编程产生并画出信号x1(n)、x2(n)、x3(n)、x4(n)、x5(n)。

（2）用matlab 编制FFT 函数对上述信号进行频谱分析，并画出上述信号谱图。

三.实验结果（1）
1.%This programm is to generate signal x1(n)=R4(n).
k=-6:6;
x=[zeros(1,6),ones(1,4),zeros(1,3)];
stem(k,x); （信号图如图1） title('图1');
2.n=-5:1:10;
x=(n+1).*(n>=0 & n<=3)+(8-n).*(n>=4 & n<=7)+0; stem(n,x)； title('图2');
3.n=-5:10;
x=(4-n).*(n>=0 & n<=3)+(n-3).*(n>=4 & n<=7); stem(n,x); title('图3');
-6
-4-20246
00.10.20.30.40.50.60.70.80.91⎪⎩⎪⎨⎧≤≤-≤≤+==n n n n n n x n R n x 其它,074,830,1)()()(241⎪⎩⎪⎨⎧≤≤-≤≤-=n n n n n n x 其它,074,330,4)(3n n x 4cos )(4
π=n n x 8
sin )(5π=图1
-5
0510
00.511.522.53
3.5
4
-5
0510
00.5
1
1.5
2
2.5
3
3.5
4
图3
4.n=-10:10; x=cos(pi/4*n); stem(n,x); title('图4');
5.n=-10:10;
x=sin(pi/8*n); stem(n,x); title('图5');
实验结果（2）: FFT 算法
function y=myditfft(x) % y=myditfft(x)
% 本程序对输入序列 x 实现DIT-FFT 基2算法，点数取大于等于x 长度的2的幂次 % x 为给定时间序列
% y 为x 的离散傅立叶变换
m=nextpow2(x);N=2^m; % 求x 的长度对应的2的最低幂次m if length(x)<N;
% 若x 的长度不是2的幂，补零到2的整数幂 x=[x,zeros(1,N-length(x))]; end
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m 数列的倒序 y=x(nxd); % 将x 倒序排列作为y 的初始值 for mm=1:m; % 将DFT 作m 次基2分解,从左到右，对每次分解作DFT 运算 Nmr=2^mm;u=1; % 旋转因子u 初始化为WN^0=1
WN=exp(-i*2*pi/Nmr); % 本次分解的基本DFT 因子WN=exp(-i*2*pi/Nmr) for j=1:Nmr/2; % 本次跨越间隔内的各次蝶形运算
for k=j:Nmr:N; % 本次蝶形运算的跨越间隔为Nmr=2^mm kp=k+Nmr/2; % 确定蝶形运算的对应单元下标 t=y(kp)*u; % 蝶形运算的乘积项 y(kp)=y(k)-t; % 蝶形运算 y(k)=y(k)+t; % 蝶形运算 end
u=u*WN; % 修改旋转因子,多乘一个基本DFT 因子WN end
-10
-8
-6
-4
-2
2
4
6
8
10
-1-0.8-0.6-0.4-0.200.20.40.60.81图4
-10
-8
-6
-4
-2
2
4
6
8
10
-1-0.8-0.6-0.4-0.200.20.40.60.81图5
end 1.
k=-6:6;
x=[zeros(1,6),ones(1,4),zeros(1,3)]; y=myditfft(x); k=-6:9; stem(k,y); xlabel('m'); ylabel('X[M]');
title('FFT 图');
2.
n=-5:1:10;
x=(n+1).*(n>=0 & n<=3)+(8-n).*(n>=4 & n<=7)+0;
y=myditfft(x); stem(n,y); xlabel('n'); ylabel('X[M]'); title('FFT 图'); 3.
n=-5:10;
x=(4-n).*(n>=0 & n<=3)+(n-3).*(n>=4 & n<=7); y=myditfft(x); stem(n,y); xlabel('n'); ylabel('X[M]'); title('FFT3'); 4.
n=-10:10;
x=cos(pi/4*n); y=myditfft(x); n=-10:21; stem(n,y); xlabel('n'); ylabel('X[M]'); title('FFT4'); 5.
n=-10:10;
x=sin(pi/8*n); y=myditfft(x); n=-10:21; stem(n,y); xlabel('n'); ylabel('X[M]'); title('FFT5');
-6
-4-20
246810
-4-3-2-101234m X [M ]
FFT 图
-5
5
10
-20-15-10-505
10
15
20
n
X [M ]
FFT 图
-5
510
-10-5
5
10
15
20
n
X M FFT3
-10
-505
10152025
-4-3-2-1012345n X [M ]
FFT4
-10
-505
10152025
-6-4
-2
24
6
8
n
X [M ]
FFT5
四.简要回答以下问题：
①在N=8时，x2(n)和x3(n)的幅频特性会相同吗？为什么？N=16呢？
答：不相同。

当N=8时，序列x1(n)和x2(n)中相同的元素值对应的 n 值是不同的，所乘的旋转因子的值也不同，因而得到的最终结果也是不同的。

同理，N=16时，所得的幅频特性也是不同的。

②FFT在什么条件下也可以用来分析周期信号序列的频谱？如果正弦信号系统sin(2πf0k),f0=0.1Hz，用16点FFT来做DFT运算，得到的频谱是信号本身的真实谱吗？为什么？
答：由于FFT算法对序列长度的要求是 N=2^M，M为正整数。

所以，当周期信号序列一个周期的长度满足 N=2^M（M为正整数）的条件时，FFT可以用来分析周期信号的频谱。

不是真实的频谱。

因为序列的周期N=10不是 2 的整数次幂，所以不是真实的。

实验三：用双线性变换法设计IIR 数字滤波器
一.实验目的：
1. 掌握用双线性变换法设计IIR DF 的原理及具体设计方法，熟悉用双线性变换法设计IIR DF 的计算机编程。

2. 观察用双线性变换法设计的DF 的频响特性，了解双线性变换法的特点。

3. 熟悉用双线性变换法设计BW 和CB 型DF 的全过程。

二，实验原理：
为了克服冲激响应不变法产生的频率混叠现象，这是从S 平面到Z 平面的标准变换z ＝esT 的多值对应关系导致的，为了克服这一缺点，产生了双线性变换法。

下式就是模拟滤波器和经采样后的数字滤波器之间的变换关系：
三.实验内容：
1.读懂所给参考程序，熟悉程序的整体结构和功能。

2.设计一个CBI 型低通DF ，通带截频fp=3000Hz ，衰耗满足Apmax=3dB ，阻带截频fS=3400Hz,衰耗ASmin=31dB,取样频率fs=8000Hz 。

画出其模拟滤波器及数字滤波器频响特性曲线。

程序：FS=8000;
Fl=3000;Fh=3400; %通带、阻带截止频率 Rp=3;Rs=31;
wp=Fl*2*pi; %临界频率采用角频率表示 ws=Fh*2*pi; %临界频率采用角频率表示 wp1=wp/FS; %求数字频率 ws1=ws/FS; %求数字频率
OmegaP=2*FS*tan(wp1/2);%频率预畸
0σ-11]Im[z j Ωj ]Re[z 0
s 平面z 平面
()()
11
211a z s T z H z H s --+=
-=
OmegaS=2*FS*tan(ws1/2);%频率预畸 %选择滤波器的最小阶数
[n,Wn]=cheb1ord(OmegaP,OmegaS,Rp,Rs,'s'); %此处是代入经预畸变后获得的归一化模拟频率参数
[bt,at]=cheby1(n,Rp,Wn,'s'); % 设计一个n 阶的契比雪夫I 型模拟滤波器 subplot(2,1,1);
[h,w]=freqz(bt,at);
plot(w*FS/(2*pi),abs(h)); grid;
xlabel('频率/Hz');ylabel('幅值'); title('模拟滤波器频响特性');
[bz,az]=bilinear(bt,at,FS); %双线性变换为数字滤波器 subplot(2,1,2);
[H,W] = freqz(bz,az); %求解数字滤波器的频率响应 plot(W*FS/(2*pi),abs(H));grid;
xlabel('频率/Hz');ylabel('幅值'); title('数字滤波器频响特性');
3
、设计一个BW 型低通DF ，满足：通带截频fp=100Hz ，衰耗满足Apmax=3dB ，阻带截频
fS=400Hz,衰耗ASmin=15dB,取样频率
fs=2000Hz 。

画出其模拟滤波器及数字滤波器频响特性曲线。

FS=2000;
Fl=100;Fh=400; %通带、阻带截止频率 Rp=3;Rs=15;
wp=Fl*2*pi; %临界频率采用角频率表示 ws=Fh*2*pi; %临界频率采用角频率表示
05001000150020002500300035004000
0.7079
0.70790.708频率/Hz
幅值
模拟滤波器频响特性
0500100015002000250030003500400000.51频率/Hz
幅值
数字滤波器频响特性
wp1=wp/FS; %求数字频率 ws1=ws/FS; %求数字频率
OmegaP=2*FS*tan(wp1/2);%频率预畸
OmegaS=2*FS*tan(ws1/2);%频率预畸 %选择滤波器的最小阶数
[n,Wn]=buttord(OmegaP,OmegaS,Rp,Rs,'s'); %此处是代入经预畸变后获得的归一化模拟频率参数
[bt,at]=butter(n,Wn,'s'); % 设计一个n 阶的契比雪夫I 型模拟滤波器 subplot(2,1,1);
[h,w]=freqz(bt,at);
plot(w*FS/(2*pi),abs(h)); grid;
xlabel('频率/Hz');ylabel('幅值'); title('模拟滤波器频响特性');
[bz,az]=bilinear(bt,at,FS); %双线性变换为数字滤波器 [H,W] = freqz(bz,az); %求解数字滤波器的频率响应 subplot(2,1,2);
plot(W*FS/(2*pi),abs(H));grid;
xlabel('频率/Hz');ylabel('幅值'); title('数字滤波器频响特性');
4
、设计一个BW 型高通DF ，满足：通带截频fp=400Hz ，衰耗满足Apmax=3dB ，阻带截频
fS=350Hz,衰耗ASmin=15dB,取样频率fs=1000Hz。

画出其模拟滤波器及数字滤波器频响特性曲线。

fp=400;fs=350;ft=1000;
010020030040050060070080090010000.9981
1.0021.004频率/Hz
幅值模拟滤波器频响特性
0100200300400500600700800900100000.51
频率/Hz
幅值
数字滤波器频响特性
rp=3;rs=15;
wp=fp/(ft/2);ws=fs/(ft/2); %利用Nyquist 频率进行归一化
[n,wc]=buttord(wp,ws,rp,rs,'s'); %求模拟滤波器的最小阶数和截止频率 [b,a]=butter(n,wc, 'high','s'); %设计高通模拟滤波器系数b ，a %Omega1=linspace(0,1000,5); %[h,w]=freqs(b,a,Omega1); [h,w]=freqs(b,a); subplot(2,1,1);
plot(w*ft/(2*pi),abs(h));grid;
title('模拟滤波器');
xlabel('频率/Hz');ylabel('幅值');
[N,Wc]=buttord(wp,ws,rp,rs); %求数字滤波器的最小阶数和截止频率 [B,A]=butter(n,wc, 'high'); %设计高通数字滤波器系数b ，a [H,W]=freqz(B,A); %绘出频率响应曲线 subplot(2,1,2);
plot(W*ft/(2*pi),abs(H));grid;
title('数字滤波器');
xlabel('频率/Hz');ylabel('幅值');
设计一个
CB 型带通DF ，满足：通带边界频率为100Hz ～500Hz ，通带衰耗小于3dB ，过渡带
宽20Hz ，阻带衰耗大于15dB ，取样频率fs=2000Hz。

画出其模拟滤波器及数字滤波器频响特性曲线。

fp=[100,500];fs=[80,520];ft=2000; rp=3;rs=15;
wp=fp/(ft/2);ws=fs/(ft/2);
02004006008001000120014001600
00.5
1
1.5模拟滤波器
频率/Hz
幅值
05010015020025030035040045050000.511.5
数字滤波器频率/Hz
幅值
[n,wn]=cheb2ord(wp,ws,rp,rs,'s'); [b,a]=cheby2(n,rs,wn,'s'); [h,w] = freqs(b,a); subplot(2,1,1);
plot(w*ft/(2*pi),abs(h));grid; title('模拟滤波器');
xlabel('频率/Hz');ylabel('幅值'); [N,Wn]=cheb2ord(wp,ws,rp,rs); [B,A]=cheby2(N,rs,Wn); [H,W] = freqz(B,A); subplot(2,1,2);
plot(W*ft/(2*pi),abs(H));grid;
title('数字滤波器');
xlabel('频率/Hz');ylabel('幅值');
●求出上述AF 低通原型的H(s)以及H(z)，检查所设计的DF 是否满足要求。

答：由实验结果分析可知，所设计的DF 满足要求。

●与脉冲响应不变设计法相比较，简述双线性变换设计法的优缺点。

答：双线性变换法能克服脉冲响应不变法的频谱混叠，适合于设计幅度响应为分段常数的数字滤波器，不适合设计幅度响应为非常数的数字滤波器。

0500100015002000
250030003500
0.511.5模拟滤波器
频率/Hz 幅值01002003004005006007008009001000
00.511.5数字滤波器
频率/Hz 幅值。