数字信号处理课程设计实验报告

第 2 页
华 北 电 力 大 学
实 验 报 告
反应信号的任何频率信息，只是在原有频谱图的基础上变得更加密集而已，这也 就是高密度频谱图。 当取样点 N=256>64 时，可以从右上角的 256 点 FFT 变换频谱图中清楚的 看出吸纳后的频率成分（6.5KHz,7KHz,9KHz）,此时的频谱图就是高分辨率频谱 图。 当取样点足够时，只是改变 FFT 变换的点数，在原有高分辨率频谱图上增 加一些叠加的频率成分，并且在信号频率的点上又叠加了一些新的幅值，所以， 信号的的频率成分又在原有的高分辨率频谱的基础上变成了高密度的频谱图。
实验程序： N=16;L=256;fs=32000;f1=6500;f2=7000;f3=9000;N2=512; T=1/fs;ws=2*pi*fs; n=0:N-1; x=cos(2*pi*f1*n*T)+5*cos(2*pi*f2*n*T)+cos(2*pi*f3*n*T); x1=x(1:N); X1=fft(x1,N); w=((0:N-1)*ws/N)/(2*pi); %横坐标范围 subplot(2,2,1);stem(w,abs(X1)); ylabel('幅频特性曲线');xlabel(' 采样点为 16 的 16 点 FFT 频率 Hz'); x2=[x(1:N) zeros(1,L-N)]; X2=fft(x2,L); w=((0:L-1)*ws/L)/(2*pi); subplot(2,2,3);stem(w,abs(X2)); ylabel('幅频特性曲线');xlabel(' 采样点为 16 补零到 256 后的 256FFT'); n=0:L-1;
华 北 电 力 大 学实 验Βιβλιοθήκη 报 告实验 环境 实验 名称
MATLAB 6.5 实验一： FFT 的应用 1、熟悉 MATLAB 在数字信号处理中的应用。 2、掌握利用 FFT 计算序列线性卷积的基本原理及编程实现。 3、掌握对连续信号进行采样的基本原理和方法，并利用 FFT 对信号进行频谱分 析。 1.对于两个序列：x(n)=nR16(n),h(n)=R8(n) （1）在同一图形窗口中绘出两序列的时域图形。 （2）利用 FFT 编程计算两序列的线性卷积，绘出的时域图形。 2.利用 FFT 对信号进行谱分析 对于连续信号 xa(t)=cos(2π f1t) +5cos(2π f2t) +cos(2π f3t) ，其中 f1=6.5kHz, f2=7kHz, f3=9kHz, 以采样频率 fs=32 kHz 对其进行采样， （1）对 xa(t) 信号采集 16 点样本，分别作 16 点和补零到 256 点的 FFT，并 分别绘出对应的幅频特性曲线。 （2）对 xa(t)信号采集 256 点样本，分别作 256 点和 512 点的 FFT，并分别绘 出对应的幅频特性曲线。 （3）比较（1）和（2）中的结果，分析采样点数和傅里叶变换点数对 FFT 的影响，说明高密度频谱和高分辨率频谱的特点与区别。 题目分析： 1.利用 FFT 计算线性卷积 首先.周期卷积是线性卷积以 L 为周期的周期延拓序列的主值序列； 两个 长度为 M，N 的序列的线性卷积可用长度均为 L 的循环卷积来代替， 用循环卷积计算线性卷积的方法归纳如下: 将长为 N2 的序列 x(n)延长到 L,补 L-N2 个零 将长为 N1 的序列 h(n)延长到 L,补 L-N1 个零 如果 L≥N1+N2-1,则循环卷积与线性卷积相等,此时,可有 FFT 计算线性 卷积,方法如下: a.计算 X(k)=FFT[x(n)] b.求 H(k)=FFT[h(n)] c.求 Y(k)=H(k)Y(k) k=0～L-1 d.求 y(n)=IFFT[Y(k)] n=0～L-1 可见，只要进行二次 FFT,一次 IFFT 就可完成线性卷积计算。计算表 明,L>32 时,上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因 此,也称上述循环卷积方法为快速卷积法。 题目利用 FFT 计算线性卷积 以下为程序及结果：
题目分析： 利用 FFT 进行频谱观测分析 频率分辨率可以理解为在使用 DFT 时，在频率轴上的所能得到的最小频率 间隔 f0=fs/N=1/NTs=1/T,其中 N 为采样点数，fs 为采样频率，Ts 为采样间隔。 因为最小的频率差值为 7KHz—6.5KHz=0.5KHz， 采样频率为 32K,所以最小 的样本数目应当是 32/0.5=64 个，当采样点不足时，必然发生混叠失真，即不能 观测出原信号的频率分布。 不是采样点数越多， 频率分辨力就越高， 因为一段数据拿来就确定了时间 T， 注意：f0=1/T，而 T=NTs，增加 N 必然减小 Ts ，因此，增加 N 时 f0 是不变的。 只有增加点数的同时导致增加了数据长度 T 才能使分辨率越好。 还有容易搞混的 一点，我们在做 DFT 时，常常在有效数据后面补零达到对频谱做某种改善的目 的，我们常常认为这是增加了 N，从而使频率分辨率变好了，其实不是这样的， 补零并没有增加有效数据的长度，仍然为 T。但是补零其实有其他好处： 1.使数据 N 为 2 的整次幂，便于使用 FFT。 2.补零后，其实是对 DFT 结果做了插值，，使谱外观平滑化。 3.由于对时域数据的截短必然造成频谱泄露， 因此在频谱中可能出现难以 辨认的谱峰，补零在一定程度上能消除这种现象。 那么选择 DFT 时 N 参数要注意： 1.由采样定理：fs>=2fh； 2.频率分辨率.f0=fs/N，所以一般情况给定了 fh 和 f0 时也就限制了 N 范 围：N>=fs/f0。 当取样点 N=16<64 时，左上角的信号的 16 点 FFT 变换所得频谱图上几乎反 应出信号的任何频率信息，因为幅值的大小并不能代表对应的频率即为信号的频 率分布。 当采样点 N=16<64 时，左下角的信号的 256 点 FFT 变换所得频谱图也不能
实 验 结 果 及 分 析
第 3 页
华 北 电 力 大 学
实 验 报 告
实验 名称
设计性实验二：IIR 数字滤波器的设计 1、本实验为设计性实验。 2、掌握用双线性变换法设计 IIR 数字滤波器的基本原理和设计方法。 3、掌握用双线性变换法设计 IIR 数字 Butterworth 滤波器的原理和设计方法。 IIR 数字滤波器的设计 用双线性变换法设计一个 IIR 数字 Butterworth 低通滤波器。 技术指标为： 通带截止频率 fp=1kHz ，阻带截止频率 fs=1.5kHz ，通带衰减 Rp≤1dB，阻 带衰减 Rs ≥40dB ，采样频率 Fs=10kHz。绘出滤波器的幅频特性曲线和相频 特性曲线，判断设计是否符合要求。 题目分析： 双线性变换法，其主要的优点： 避免了频率响应的混叠现象。这是因为 S 平面与 Z 平面是单值的一一对 应关系。S 平面整个 jΩ轴单值地对应于 Z 平面单位圆一周，即频率轴是单值 变换关系。在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近 于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时， ω终止在折叠频率ω=π处， 因而双线性变换就不会出现由于高频部分超过折 叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。 但是双线性变换的这个特点是靠频率的严重非线性关系而得到的。 由于 这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位 的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有 的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分 段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型 的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字 滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变， 双线性变换法幅 度和相位特性的非线性映射 对于分段常数的滤波器，双线性变换后，仍得 到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸 变，这种频率的畸变，可以通过频率的预畸来加以校正。 预畸变就是：将临界模拟频率事先加以畸变，然后经变换后正好映射到 所需要的数字频率上。 由于用脉冲响应不变法设计 IIR 数字滤波器，存在从 S 平面到 Z 平面的 多值映射的缺点，而且还有频谱混叠现象，所以本次课程设计就采用双线性 变换法来设计 IIR 数字滤波器。 设计滤波器的步骤： 1.得到数字指标（Wn 等） 2.双线性变换为模拟低通指标 3.归一化模拟拟指标 4.利用通带衰减与阻带衰减的值求 Butterworth 的阶数 N 及归一 化的模拟系统函数
实 验 目 的
实 验 内 容
第 4 页
华 北 电 力 大 学
实 验 报 告
5.将这个再经过去归一化得到想要的滤波器类型 6.用双线性变换法变为数字滤波器 程序代码： fp=1000;%通带截止频率 fs=1500;%阻带截止频率 Rp=1;Rs=40;Fs=10000; wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; T=1;Fs=1/T;Wp=(2/T)*tan(wp/2);Ws=(2/T)*tan(ws/2);% 数 字 域 转 模 拟 域 %预畸变 [N,Wn]=buttord(Wp,Ws,Rp,Rs,'s');%求得滤波器的阶数 N 与 3DB 点对应的频 率 Wn [Z,P,K]=buttap(N);%建立传递函数 Ha(s)。 [bz,az]=zp2tf(Z,P,K); [b,a]=lp2lp(bz,az,Wn); [B,A]=bilinear(b,a,Fs); [H,w]=freqz(B,A,512);%从帮助里面找到返回的单位是 Hz，而不是 rad/s H=(H(1:501))';w=(w(1:501))'; mag=abs(H);pha=angle(H); subplot(2,1,1);plot(w/pi,mag);ylabel('幅度');xlabel('以\pi 为单位的 频率'); title('幅频特性曲线');grid; subplot(2,1,2);plot(w/pi,pha);title('相频特性曲线'); ylabel('相位');xlabel('以\pi 为单位的频率'); grid;