发布时间 : 星期四 文章数字信号处理实验，滤波器设计更新完毕开始阅读20dcd3523c1ec5da50e270de
数字信号处理实验指导书（修订版1）
f=fs/N*(0:N/2-1); % subplot(311);plot(x); % subplot(312);plot(f,abs(Y(1:N/2))); % subplot(313);plot(f,PSD(1:N/2)); % 1)画出图形窗口显示的图形，并注名每个图形的含义。 2)回答下列问题：
i)观察幅频谱图，可以发现，信号x(n)含有的两个频率分量分别是 Hz和 Hz。
ii)在该程序中的“f=fs/N*(0:N/2-1);”下添加“k=0:N/2-1;”， “plot(f,abs(Y(1:N/2)));”改为“plot(k,abs(Y(1:N/2)));” 重新运行该程序并观察幅频谱图，图中两峰值对应的下标分别是 和 。它们的含义为 。
再将该程序中的N改为512, 重新运行该程序并观察幅频谱图，这时图中
两峰值对应的下标分别是 和 。结果是否和上面的相同？ 为什么？ 。
iii)本例的频谱分辨率F是 Hz,改变f2=60Hz,问：在幅频谱中，能否分辨f1和f2 对应的频率分量？ 。为什么？ 。 再改变f2=52Hz,问：在幅频谱中，能否分辨f1和f2对应的频率分量？ 。 为什么？ 。
再改变f2=600Hz，在幅频谱中，f2对应的频率分量出现在 Hz； 问：在fs=1000Hz的情况下，能否正确检测f2对应的频率分量？ 。
为什么？ 。
为了正确检测f2对应的频率分量，则fs至少取多少Hz? Hz。在该程序中改变fs，验证你的结论。
iv)比较幅频谱和功率谱，可以发现功率谱具有 的特性。 3． FFT实现任意两个序列的快速卷积。 %程序 fftjuanji.m
clear clc
x1=input('x1=');x2=input('x2='); % N1=length(x1);N2=length(x2); %序列x1(n),x2(n)的长度 E=ceil(log2(N1+N2-1)); ?il---向+∞方向取整 N=2^E; % x1=[x1,zeros(1,N-N1)]; % x2=[x2,zeros(1,N-N2)];
X1=fft(x1,N); % X2=fft(x2,N);
Y=X1.*X2; % y=ifft(Y,N) % 结果分析： 1)回到MATLAB窗口，键入：
x1=[1 1 1], x2=[1 2],回车。
结果:y= 2)问：可用Matlab中的什么函数验算上述卷积结果？ 9
数字信号处理实验指导书（修订版1）
4. 利用谱分析观察太阳黑子周期性。
以100年中记录到的太阳黑子出现次数为信号x(n),对x(n)作功率谱，从中观察太阳黑子周期性。
%程序 taiyangheizi.m clear
clc
x=[101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 ... 132 131 118 90 67 60 47 41 21 16 6 4 7 14 34 45 43 48 ... 42 28 10 8 2 0 1 5 12 14 35 46 41 30 24 16 7 4 2 8 ... 17 36 50 62 67 71 48 28 8 13 57 122 138 103 86 63 37 24 ... 11 15 40 62 98 124 96 66 64 54 39 21 7 4 23 55 94 96 ... 77 59 44 47 30 16 7 37 74];0年中太阳黑子出现的次数
subplot(211);plot(x) %画x(n)
N=128; fs=1; %fs=1Hz,N=128点
s=x-mean(x); %对x作零均值化处理（去除直流分量） Y= ; %对s做N点fft
PSD= ; %做功率谱PSD f= ; %将频率定标为实际频率f subplot(212); ; %画功率谱（N/2点） 1)填写空格中的画图语句并绘出结果图形。
2) 从s的功率谱观察到，其幅度最高处对应的横坐标f= Hz, 则太阳黑子每隔 年出现一次最高峰。 3)在对s做FFT时，为何可取fs=1Hz,N=128点？
10
数字信号处理实验指导书（修订版1）
实验三 IIR数字滤波器设计
一、实验目的
(1)掌握用双线性变换法设计IIR数字低通和高通滤波器。 (2)设计低通滤波器对实际心电图信号进行滤波。
(3)设计低通滤波器对含有啸叫噪声的音乐信号进行消噪。
(4)设计IIR数字低通和高通滤波器对某个DTMF(双音多频)信号进行频带分离。 二、实验环境
1.Windows98以上操作系统 2.安装MATLAB6.0以上版本 三、实验原理
1．选频型数字滤波器的种类有 、 、 和 滤波器。 2. 从实现方法上，数字滤波器通常分为 和 滤波器。
3．IIR滤波器的设计目的是根据技术指标，找到 ；IIR滤波的MATLAB语句为y= ;
4．双线性变换法设计低通IIR滤波器的步骤是
(1) (2) (3) 四、实验内容
1. 人体心电图信号在测量过程中往往受到工业高频干扰， 必须经过低通滤波处理后才能作为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(n)， 其中存在高频干扰。 试以x(n)作为输入序列， 滤除其中的干扰成分。
x(n)= {-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0} 。 低通滤波器设计指标：ωp=0.2πrad,ωs=0.3πrad,ap=1dB,as=15dB。
H(z)?0.0007378(1?z)(1?1.268z?1?16?2?1?0.705z?2)(1?1.0106z?1?0.3583z)(1?0.904z?0.215z?2)已设计出H(z) (p300)
3??K?1Hk(z)x(n)H1(z)y1(n)H2(z)H(z)y2(n)H3(z)y3(n)＝y(n)Hk(z)?A(1?2z1?Bkz?1?1?z?2)?Ckz?2,k?1,2,3其中
A=0.09036 ；B1=1.2686,C1=-0.7051 ；B2=1.0106,C2=-0.3583；B3=0.9044,C3=-0.2155 *IIR滤波的Matlab语句：y=filter(b,a,x)
b,a----Hk(z)分子/分母系数；x---输入信号x(n)；y---滤波结果y(n)。 **求频响特性的Matlab语句：[H,w]=freqz(b,a,N)
b,a----Hk(z)分子/分母系数；N---频率点数；H—滤波器的频响特性H(w); w---数字频率轴w。
(1)用IIR低通滤波器(原理设计)对实际的心电信号进行滤波 (IIRECG.m)
11
数字信号处理实验指导书（修订版1）
%程序：(IIRECG.m)
clear clc
x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6, ... -4, -4, -6, -6,-2,6,12,8,0,-16, ... -38,-60,-84,-90,-66,-32,-4, -2, -4, 8, ... 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, ...
0, 0, -2, -4, 0, 0, 0, -2, -2, 0, ...
0, -2, -2, -2, -2, 0];%心电信号x(n)[含高频噪声] A=0.09036; %IIR低通滤波器系数 B1=1.2686;C1=-0.7051; B2=1.0106;C2=-0.3583; B3=0.9044;C3=-0.2155; b=[A 2*A A]; a1=[1 -B1 -C1]; N=128;
[H1,w]=freqz(b,a1,N);%频响特性H1(w)
y1=filter(b,a1,x); %y1(n)是x(n)经H1(z)滤波的结果 a2=[1 -B2 -C2];
[H2,w]=freqz(b,a2,N);%频响特性H2(w)
y2=filter(b,a2,y1); %y2(n)是y1(n)经H2(z)滤波的结果 a3=[1 -B3 -C3];
[H3,w]=freqz(b,a3,N);%频响特性H3(w)
H=H1.*H2.*H3; %总的频响特性H(w)=H1(w)H2(w)H3(w) mag=abs(H); %滤波器的幅频特性
db=20*log10((mag+eps)/max(mag));%幅频特性(dB)
y3=filter(b,a3,y2); %y3(n)是y2(n)经H3(z)滤波的结果 X=fft(x,N); %对x(n)做N点fft,得其频谱X
wx=2*pi*(0:N/2-1)/N; %将坐标轴从频率点k转换为数字频率wx(wx=2*pi*k/N) %%%绘图%%%
subplot(221);plot(x);grid on;title('x(n)');
subplot(222);plot(wx/pi,abs(X(1:N/2)));grid on;title('|X(wx)|'); subplot(223);plot(w/pi,db);grid on;title('|H(w)|(db)'); subplot(224);plot(y3);grid on;title('y3(n)'); 要求：
绘出并分析结果图形。
(2)用现有的Matlab函数设计上述IIR低通滤波器,对实际的心电信号进行滤波，同样绘出题1的结果图形。 Matlab函数:
*根据技术指标ωp ,ωs ,αp ,αs计算巴特沃思滤波器的阶数N和3dB截止频率ωc语句：[N,wc]=buttord(wp/pi,ws/pi,ap,as) 注：数字频率以π为单位。
**根据N,ωc确定数字滤波器H(z)分子/分母多项式的系数b,a 低通 [b,a]=butter(N,wc)
12