第一篇：数字信号处理实验-FFT的实现
实
验
报
告
学生姓名：
学 号：
指导教师：
一、实验室名称：数字信号处理实验室
二、实验项目名称：FFT的实现
三、实验原理：
一．FFT算法思想：
1．DFT的定义：
对于有限长离散数字信号{x[n]}，0 n N-1，其离散谱{x[k]}可以由离散付氏变换（DFT）求得。DFT的定义为：
N1X[k]通常令ej2Nx[n]en0j2Nnk，k=0,1,…N-1 WN，称为旋转因子。
2．直接计算DFT的问题及FFT的基本思想：
由DFT的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（N-1）次加法。因此，对于一些相当大的N值（如1024）来说，直接计算它的DFT所作的计算量是很大的。
FFT的基本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。例如，若N为偶数，将原有的N
22点序列分成两个（N/2）点序列，那么计算N点DFT将只需要约[(N/2)·2]=N/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。
3．基2按时间抽取（DIT）的FFT算法思想：
设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。
将长度为N2L的序列x[n](n0,1,...,N1)，先按n的奇偶分成两组：
x[2r]x1[r]x[2r1]x2[r]，r=0,1,…,N/2-1 DFT化为：
N1N/21N/21X[k]DFT{x[n]}N/21n0x[n]WnkN2rkr0x[2r]W2rkNr0x[2r1]WN(2r1)kN/21r0N/21x1[r]Wx1[r]W2rkNWWkNr0N/21x2[r]WN
r0rkN/2kNr0x2[r]WN/22rkrk上式中利用了旋转因子的可约性，即：WNN/21NrkN/21rkWN/2。又令
rkX1[k]r0x[1r]W,/X2[k]2r0x[r]WN2，则上式可以写成： /2X[k]X1[k]WNX2[k](k=0,1,…,N/2-1)
k可以看出，X1[k],X2[k]分别为从X[k]中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了X[k]前N/2点的DFT值，还需要继续利用X1[k],X2[k]表示X[k]的后半段本算法推导才完整。利用旋转因子的周期性，有：WN/2WN/2X1[N2N/21rkr(kN/2)，则后半段的DFT值表达式：
rkk]r0x1[r]W2N/2r(Nk)N/21r0x1[r]WN/2X1[k]，同样，X2[N2k]X2[k]
（k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表达式获得，中间还利用到WN(N2k)NWN2Wk得到后半段的X[k]值表达式W，k为：X[k]X1[k]WNkX2[k]（k=0,1,…,N/2-1）。
这样，通过计算两个N/2点序列x1[n],x2[n]的N/2点DFTX1[k],X2[k]，可以组合得到N点序列的DFT值X[k]，其组合过程如下图所示：
X1[k] X1[k]WNkX2[k]
X2[k] WNnk-1 X1[k]WNkX2[k]
比如，一个N = 8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：
x(0)W0x(1)W0x(2)W0x(3)W2W0W1W0x(5)W0x(6)W0x(7)W2X(7)W3X(6)W2X(5)X(3)X(2)X(1)X(0)x(4)X(4)
4．基2按频率抽取（DIF）的FFT算法思想：
设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。
在把X[k]按k的奇偶分组之前，把输入按n的顺序分成前后两半：
N1N/21nkNN1X[k]DFT{x[n]}N/21N/21x[n]Wn0(nn0N2)kx[n]WnkNnN/2x[n]WNnkn0N/21x[n]WnkNn0x[nNkN2]WNnk
Nn0[x[n]x[nN2NkN2]W2N]WN,k0,1,...,N1因为W2N1，则有WX[k](1)，所以：
kkN/21n0[x[n](1)x[nN2]]WN,k0,1,...,N1
nk按k的奇偶来讨论，k为偶数时：
N/21X[2r]n0[x[n]x[nN2]]WN,k0,1,...,N1 N22rnN/21k为奇数时：X[2r1]前面已经推导过WNN/21n0[x[n]x[n]]WN(2r1)n,k0,1,...,N1
2rkWN/2，所以上面的两个等式可以写为：
N2]]WN/2,r0,1,...,N/21 N2rnrkX[2r]n0[x[n]x[nN/21X[2r1]n0{[x[n]x[n]]WN}WN/2,r0,1,...,N/21
nnr通过上面的推导，X[k]的偶数点值X[2r]和奇数点值X[2r1]分别可以由组合而成的N/2点的序列来求得，其中偶数点值X[2r]为输入x[n]的前半段和后半段之和序列的N/2点DFT值，奇数点值X[2r1]为输入x[n]的前半段和后半段之差再与WN相乘序列的N/2点DFT值。
令x1[n]x[n]x[nN/21nN2]，x2[n][x[n]x[nN/21N2]]WN，则有：
nX[2r]n0x1[n]WrnN/2,X[2r1]n0x2[n]WrnN/2,r0,1,...,N21
这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示：
x[n] x[n]x[nN2]
x[nN2]-1 WNn [x[n]x[nN2]]WNn
二．在FFT计算中使用到的MATLAB命令：
函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R 四、实验目的： 离散傅氏变换（DFT）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。 本实验通过直接计算DFT，利用FFT算法思想计算DFT，以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。 五、实验内容： a)计算实数序列x(n)cos516n,0n256的256点DFT。 b)计算周期为1kHz的方波序列（占空比为50％，幅度取为＋/-512，采样频率为25kHz，取256点长度）256点DFT。 六、实验器材（设备、元器件）： 安装MATLAB软件的PC机一台，DSP实验演示系统一套。 七、实验步骤： （1）先利用DFT定义式，编程直接计算2个要求序列的DFT值。 （2）利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值。（3）(拓展要求)不改变序列的点数，仅改变DFT计算点数（如变为计算1024点DFT值），观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。 （4）利用FFT的基本思想（基2－DIT或基2－DIF），自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果进行对比。 （5）(拓展要求)尝试对其他快速傅立叶变换算法（如Goertzel算法）进行MATLAB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行对比。 （6）（拓展要求）在提供的DSP实验板上演示要求的2种序列的FFT算法（基2－DIT），用示波器观察实际计算出来的频谱结果，并与理论结果对比。 八、实验数据及结果分析： 程序：（1）对要求的2种序列直接进行DFT计算的程序 （2）对要求的2种序列进行基2－DIT和基2－DIF FFT算法程序（3）对要求的2种序列用MATLAB中提供的FFT函数进行计算的程序 结果：（1）对2种要求的序列直接进行DFT计算的频域波形 （2）对2种要求的序列进行基2－DIT和基2－DIF FFT算法频域波形（3）对2种要求的序列用MATLAB中提供的FFT函数计算的频域波形。（4）（拓展要求）分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。 （5）（拓展要求）保持序列点数不变，改变DFT计算点数（变为1024点），观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。 九、实验结论： 十、总结及心得体会： 十一、对本实验过程及方法、手段的改进建议： 数字信号处理 说 明 书 目录 一．摘要…………………………………3 二．课程设计目的………………………3 三．设计内容……………………………3 四．设计原理……………………………4 4.1．语音信号的采集…………………………….4 4.2．滤波器……………………………………….4 4．21．IIR滤波器原理…………………………………….4 4.22．FIR滤波器原理………………………………………5 五．设计步骤……………………………6 5.1录制女音………………………………………6 5.2采样语音信号并画出时域波形和频谱图……7 5.3采用双线性变换法设计IIR滤波器…………10 5.4窗函数法设计FFR滤波器………………......12 5.5用IIR滤波器对信号进行滤波………………14 5.6用FIR滤波器对信号进行滤波………………16 5.7男女声语音信号频谱特点分析………………19 5.8有背景噪声的信号分析………………………20 六．心得体会…………………………….22 七．参考文献…………………………….23 一．摘要： 这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。 MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。 在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。 二．课程设计目的： 综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。 三．设计内容： 内容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数 法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。 四．设计原理： 4.1．语音信号的采集 熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有： a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据； b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y； c：wavwrite((y,fs,wavfile);创建音频文件； d：y=wavread(file);读取音频文件； 关于声音的函数还有sound()；soundsc()；等。4.2滤波器： 4．21.IIR滤波器原理 冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。 双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也 就消除了频谱混叠现象。 双线性法设计IIR数字滤波器的步骤： 1）将数字滤波器的频率指标{ k}由Wk＝（2/T）*tan(wk),转换为模拟滤波器的频率指标{k}.2)由模拟滤波器的指标设计H(s).3)由H(s)转换为H(z)21z1H(z)H(s)sT1z1 4.22.FIR滤波器原理 FIR滤波器与IIR滤波器特点不同，IIR滤波器的相位是非线性的，若需线性相位则要采用全通网络进行相位校正。而有限长单位冲激响应（FIR）数字滤波器就可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。 由于FIR系统的冲激响应就是其系统函数各次项的系数，所以设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲激响应作为H（z）的系数，冲激响应长度N就是系统函数H（z）的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。这种时域设计、频域检验的方法一般要反复几个回合，不像IIR DF设计靠解析公式一次计算成功。给出的理想滤波器频率响应是，它是w的周期函数，周期 由傅立叶反变换导出，即 hd(n)1Hd(ejw)ejwndw2，再将hd(n)与窗函数，因此可展开成傅氏级数w(n)相乘就可以得到h(n)。、的计算可采用傅氏变换的现成公式和程序，窗函数也是现成的。但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算，整个设计可用计算机编程来做。 窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有： 矩形窗 w(n)=RN(n)； Hanning窗 ； Hamming窗 ； Blackmen窗 ； Kaiser窗。 式中Io(x)为零阶贝塞尔函数。 五．设计步骤： 5.1录制女音： 利用MATLAB中的函数录制声音。function nvyin()fs=11025; %采样频率 str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 开始录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double'); %录制声音3秒 str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 播放录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs); %播放录音 str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp(' 播放录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音'); %声音的存储 5.2采样语音信号并画出时域波形和频谱图 读取语音信号，画出其时域波形和频谱图，与截取后的语音信号的时域波形和频谱图比较，观察其变化。程序如下： [x,fs,bits]=wavread('女音.wav'); %读取声音 N=length(x); %计数读取信号的点数 t=(1:N)/fs; %信号的时域采样点 f0=fs/N; %采样间隔 n=1:N/2; %取信号的一半 figure(1);subplot(2,2,1); %把画图区域划分为2行2列，指定第一个图 plot(t, x); %画出声音采样后的时域波形 title('原女音信号的时域波形'); %给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid; %添加网格 y=fft(x); %对信号做N点FFT变换 k=(n-1)*f0; %频域采样点 subplot(2,2,3); %把画图区域划分为2行2列，指定第三个图 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid; %添加网格 subplot(2,2,4); %把画图区域划分为2行2列，指定第四个图 if y~=0 %判断指数是否为0 plot(k,20*log10(abs(y(n)))); %画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid; %添加网格 %实际发出声音落后录制动作半拍的现象的解决 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]); %截取语音信号 N=length(x1); %计数读取信号的点数 t=(1:N)/fs; %信号的时域采样点 f0=fs/N; %采样间隔 n=1:N/2; %取信号的一半 figure(2);subplot(2,2,1); %把画图区域划分为2行2列，指定第一个图 plot(t,x1); %画出声音采样后的时域波形 title('截取后女音信号的时域波形'); %给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid; %添加网格 y1=fft(x1); %对信号做N点FFT变换 subplot(2,2,3); %把画图区域划分为2行2列，指定第三个图 k=(n-1)*f0; %频域采样点 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid; %添加网格 subplot(2,2,4); %把画图区域划分为1行2列，指定第二个图 if y1~=0 %判断指数是否为0 plot(k,20*log10(abs(y1(n)))); %画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid; %添加网格 原女音信号的时域波形10.5A/幅0振-0.5-10123时间/tFFT变换后声音的频谱特性FFT变换后声音的频谱特性30050A200贝/值分/幅0幅100振00200040006000-***频率/HzHz 截取后女音信号的时域波形10.5振幅/A0-0.5-10123FFT变换后声音的频谱特性50时间/tFFT变换后声音的频谱特性300200振幅/分贝幅值/A01000020004000频率/Hz6000-5002000Hz40006000 结果分析： 由原女音信号的时域波形可知录取开始时实际发出声音大概落后3500个采样点，我们把前3500点去除即可解决实际发出声音落后录制动作半拍的现象。由原女音的的频谱图和截取后声音的频谱图可看出，对声音的截取并不会影响它们频谱分布。 5.3采用双线性变换法设计IIR滤波器： 人的声音频率一般在（1~~4）kHZ之间，则我们只需要设计一个带通滤波器即可滤去声音频带以外的无用噪声，得到比较清晰的声音。根据声音的频谱图分析，设计一个带通滤波器性能指标如下： fp1＝1000 Hz，fp2＝3000 Hz，fsc1＝500 Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000 程序如下： %iir带通的代码： %w=2*pi*f/fs Ap=1; %通带波纹系数 Az=100; %最小阻带衰减 wp=[0.2 0.6]; %归一化通带数字截止频率 wz=[0.1 0.7]; %归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az); %估计契比雪夫I型滤波器阶数 [b,a]=cheby1(N,Ap,wn); %N指定滤波器阶数，wn归一化 截 %止频率，Ap通带波动 [h,w]=freqz(b,a); %求数字滤波器的复频率响应 figure(1);subplot(2,1,1);plot(w/pi,abs(h)); %绘制数字滤波器的频谱图 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型带通滤波器的幅频响应');subplot(2,1,2);if abs(h)~=0 %判断指数是否为0 plot(w/pi,20*log10(abs(h))); %绘制数字滤波器频谱的分贝图 end grid;xlabel('omega/pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型带通滤波器的幅频响应');契比雪夫Ⅰ型带通滤波器的幅频响应1振幅（幅值）0.5000.10.20.50.60.70.8/契比雪夫Ⅰ型带通滤波器的幅频响应0.30.40.910振幅(分贝)-200-400-60000.10.20.30.40.5/0.60.70.80.91 5.4窗函数法设计FFR滤波器 线性相位FIR滤波器通常采用窗函数法设计。窗函数法设计FIR滤波器的基本思想是：根据给定的滤波器技术指标，选择滤波器长度N和窗函数ω（n），使其具有最窄宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h(n)。工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯泽（Kaiser）窗。 这次设计我采用的是布莱克曼来设计给定数字带通滤波器的参数如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下： Ap=1; %通带波纹系数 Az=100; %最小阻带衰减 fs=10000; %采样频率 wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2)); %---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW); %---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2); %---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想带通滤波器 wdWindow=blackman(N); %布拉克曼窗 hr=wdWindow.*hdWindow'; %点乘 n=0:N-1; %阶数 subplot(2,2,1);stem(n,wdWindow); %绘制布拉克曼窗时域波形 xlabel('时间');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1); %求滤波器频率响应 subplot(2,2,3);plot(W/pi,abs(H)) %绘制滤波器频域波形 xlabel('omega/pi');ylabel('振幅');title('FIR带通滤波器幅频特性');subplot(2,2,4); if abs(H)~=0 %判断指数是否为0 plot(W/pi,20*log10(abs(H))); %画滤波器频谱的分贝图 end xlabel('omega/');ylabel('振幅/分贝');title('FIR带通滤波器幅频特性');grid; %添加网格 %---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录ideallp.m function hd = ideallp(wc,N);% 理想低通滤波器的脉冲响应子程序 % hd = 点0 到 N-1之间的理想脉冲响应 % wc = 截止频率(弧度)% N = 理想滤波器的长度 tao =(N-1)/2; % 理想脉冲响应的对称中心位置 n = [0:(N-1)]; % 设定脉冲响应长度 m = n-tao + eps; % 加一个小数以避免零作除数 hd = sin(wc*m)./(pi*m); % 理想脉冲响应 布拉克曼窗1振幅0.500406080时间FIR带通滤波器幅频特性500振幅/分贝20FIR带通滤波器幅频特性1.51振幅-50-100-15000.5/10.5000.5/1 5.5用IIR滤波器对信号进行滤波 用自己设计的IIR滤波器分别对采集的信号进行滤波，在Matlab中，IIR滤波器利用函数filter对信号进行滤波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x); %计数读取信号的点数 t=(1:N)/fs; %信号的时域采样点 f0=fs/N; %采样间隔 n=1:N/2; %取信号的一半 y=fft(x); %对信号做N点FFT变换 k=(n-1)*f0; %频域采样点 subplot(2,1,1); %把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid; %iir带通的代码： Ap=1; %通带波纹系数 Az=100; %最小阻带衰减 wp=[0.2 0.6]; %归一化通带数字截止频率 wz=[0.1 0.7]; %归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az); %估计契比雪夫I型滤波器阶数 [b,a]=cheby1(N,Ap,wn); %N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x); %对声音滤波 wavplay(x1)wavwrite(x1,'IIR滤波后女音.wav');N=length(x1); %计数读取信号的点数 t=(1:N)/fs; %信号的时域采样点 f0=fs/N; %采样间隔 n=1:N/2; %取信号的一半 y=fft(x1); %对信号做N点FFT变换 k=(n-1)*f0; %频域采样点 subplot(2,1,2); %把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid; 滤波前女音的频谱特性300幅值/A***030004000频率/Hz滤波后女音的频谱特性500060006040幅值/A***0频率/Hz400050006000 结果分析： 由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3400Hz的频谱成分。回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。 5.6用FIR滤波器对信号进行滤波 用自己设计的FIR滤波器分别对采集的信号进行滤波，在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波 程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x); %计数读取信号的点数 t=(1:N)/fs; %信号的时域采样点 f0=fs/N; %采样间隔 n=1:N/2; %取信号的一半 y=fft(x); %对信号做N点FFT变换 k=(n-1)*f0; %频域采样点 subplot(2,1,1); %把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性'); %给图形加注标签说明 xlabel('频率/omega');ylabel('幅值/A');grid; %FIR带通滤波器代码 fs=10000;wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2)); %---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW); %---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2); %---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);wdWindow=blackman(N);hr=wdWindow.*hdWindow';x1=fftfilt(hr,x); %对声音滤波 wavplay(x1)wavwrite(x1,'FIR滤波后女音.wav');N=length(x1); %计数读取信号的点数 t=(1:N)/fs; %信号的时域采样点 f0=fs/N; %采样间隔 n=1:N/2; %取信号的一半 y=fft(x1); %对信号做N点FFT变换 k=(n-1)*f0; %频域采样点 subplot(2,1,2); %把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid; 滤波前女音的频谱特性300200幅值/A***004000频率/l滤波后女音的频谱特性500060006040幅值/A20005001000***03000频率/Hz***0 结果分析： 由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3500Hz的频谱成分。和用IIR滤波器滤波一样，回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。5.7男女声语音信号频谱特点分析 换一个男音录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点。程序如下： [x,fs,bits]=wavread('女音.wav'); %读取声音 N=length(x); %计数读取信号的点数 t=(1:N)/fs; f0=fs/N; n=1:N/2; y=fft(x); k=(n-1)*f0; subplot(2,1,1); plot(k,abs(y(n))); title('FFT变换后女音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid; [x,fs,bits]=wavread('明明.wav'); N=length(x); t=(1:N)/fs; f0=fs/N; n=1:N/2; y=fft(x); k=(n-1)*f0; subplot(2,1,2); plot(k,abs(y(n))); title('FFT变换后男音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid; %信号的时域采样点 %采样间隔 %取信号的一半 %对信号做N点FFT变换 %频域采样点 %把画图区域划分为2行1列，指定第一个图%绘制原始语音信号的幅频响应图 %给图形加注标签说明 %添加网格 %读取声音 %计数读取信号的点数 %信号的时域采样点 %采样间隔 %取信号的一半 %对信号做N点FFT变换 %频域采样点 %把画图区域划分为2行1列，指定第二个图%绘制原始语音信号的幅频响应图 %给图形加注标签说明 %添加网格 axis([0 6000 0 300])； %改变横纵坐标便于比较频谱图 FFT变换后女音的频谱特性300200幅值/A***00频率/FFT变换后男音的频谱特性***200幅值/A***00频率/400050006000 结果分析： 通过比较上面女音频谱图和男音频谱图可知，男音的频谱集中在低频部分，高频成分底，谱线较平滑，声音听起来低沉。5.8有背景噪声的信号分析 从硬盘中把一段噪声（频谱能量集中在某个小范围内）叠加到语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除； 程序如下： z=wavread('女音.wav',[1 24000]); %读取声音在1-24000之间 f=wavread('noise.wav',[1 24000]);x=z+f;wavplay(x);fs=11025;N=length(x);f0=fs/N; %采样间隔 n=1:N; %取信号的一半 y=fft(x,N);%对信号做N点FFT变换 k=(n-1)*f0; %频域采样点 subplot(2,1,1); %把画图区域划分为1行2列，指定第二个图 plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图 title('加噪声后声音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格 %iir带通滤波器的代码： Ap=1; %通带波纹系数 Az=70; %最小阻带衰减 wp=[0.2 0.7]; %归一化通带数字截止频率 wz=[0.1 0.8]; %归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az); %估计契比雪夫I型滤波器阶数 [b,a]=cheby1(N,Ap,wn); %N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x); %对声音滤波 wavplay(x1)； wavwrite(x1,'滤除噪音后女音.wav');N=length(x1);f0=fs/N; %采样间隔 n=1:N; %取信号的一半 y1=fft(x1,N); %对信号做fs点FFT变换 subplot(2,1,2); %把画图区域划分为1行2列，指定第二个图 k=(n-1)*f0; %频域采样点 plot(k,abs(y1(n))); %绘制原始语音信号的幅频响应图 title('滤除噪声后声音的频谱特性'); %给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格 加噪声后声音的频谱特性3000幅值/A***0008000频率/Hz滤除噪声后声音的频谱特性***030幅值/A***000频率/Hz80001000012000 结果分析 观察加噪声后声音的频谱图可知，噪音频率主要在4000Hz处，只要我们设计一个，滤波器滤除大概在4000Hz的频谱即可，回放滤波后的语音信号，可证噪音基本滤除。 六．心得体会： 通过这次课程设计，让我对MATLAB的基本应用有了更深的了解，还有数字信号处理在MATLAB中的一些函数的用法。通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。 在这次实验中，也遇到了很多问题，比如画信号频谱的分贝图时（20*log10(abs(y))）指数为零时的处理。滤波器的设计也花了好大的功夫，刚开始不会设计参数，一头雾水，通过同学的指导和讨论，得知通过观察信号的频谱图，看噪音频率集中在那一部分，设计滤波器把其滤除即可。可反复设置参数直到滤波后语音信号的效果好为止。 七．参考文献： （1）《MATLAB LabVIEW SystemView》翁剑枫 叶志前 编著, 机械工业出版社； （2）《MATLAB及在电子信息课程中的应用》陈怀琛 吴大正 高西全编著，电子工业出版社； （3）《MATLAB在数字信号处理中的应用》（弟2版）薛年喜 编著，清华大学出版社； （4）《MATLAB扩展编程》何强 何英 编著，清华大学出版社；（5）《MATLAB7简明教程》吴清 曹辉林 编著，清华大学出版社；（6）MATLAB5.3精要.编程及高级应用》程卫国 冯峰 王雪梅 刘艺 编著，机械工程出版社。 邯 郸 学 院 讲 稿 2010 ～2011 学年 第 一 学期 分院(系、部)： 信息工程学院 教 研 室： 电子信息工程 课 程 名 称： 数字信号处理 授 课 班 级： 07级电子信息工程 主 讲 教 师： 王苗苗 职 称： 助教（研究生） 使 用 教 材： 《数字信号处理》 制 作 系 统： Word2003 邯郸学院制 实验一..Matlab仿真软件介绍 一、实验目的 熟悉Matlab仿真软件 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、学习Matlab仿真软件的安装 2、熟悉Matlab仿真软件的操作环境 3、直接在Matlab仿真软件的命令窗口实现数值计算 4、编写M文件 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 1、熟悉Matlab仿真软件 2、参阅Matlab及在电子信息类课程中的应用（第2版）唐向宏 电子工业出版社 实验二 离散信号和系统分析的Matlab实现 一、实验目的 1、Matlab实现离散信号和系统分析 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、利用Matlab产生离散信号 2、利用Matlab计算离散卷积 3、利用Matlab求解离散LTI系统响应 4、利用Matlab计算DTFT 5、利用Matlab实现部分分式法 6、利用Matlab计算系统的零极点 7、利用Matlab进行简单数字滤波器设计 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验三 利用Matlab实现信号DFT的计算 一、实验目的 1、Matlab实现信号DFT的计算 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、利用Matlab计算信号的DFT 2、利用Matlab实现由DFT计算线性卷积 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验四 利用Matlab实现滤波器设计 一、实验目的 1、Matlab实现实现滤波器设计 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、利用Matlab实现模拟低通滤波器的设计 2、利用Matlab实现模拟域频率变换 3、利用Matlab实现脉冲响应不变法 4、利用Matlab实现双线性变换法 5、利用Matlab实现数字滤波器设计 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验五 利用Matlab实现FIR滤波器设计 一、实验目的 1、Matlab实现实现滤波器设计 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、利用Matlab实现窗函数法 2、利用Matlab实现频率取样法 3、利用Matlab实现优化设计 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验六..随机信号功率谱估计的Matlab实现 一、实验目的 1、Matlab实现实现滤波器设计 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、利用Matlab实现随机序列 2、利用Matlab计算相关函数的估计 3、利用Matlab进行非参数功率谱估计 4、利用Matlab进行AR模型功率谱估计 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验七..数字滤波器结构的Matlab实现 一、实验目的 1、Matlab实现实现滤波器设计 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、利用Matlab实现数字滤波器直接型设计 2、利用Matlab实现数字滤波器级联设计 3、利用Matlab实现数字滤波器并联型设计 4、利用Matlab实现数字滤波器格型设计 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验八....利用Matlab实现信号小波分析 一、实验目的 1、Matlab实现实现滤波器设计 2、进一步熟悉Matlab软件操作 二、实验设备和元器件 含Matlab仿真软件的计算机 三、实验内容和步骤 1、小波测试信号 2、分解与重构滤波器组 3、离散小波变换 4、离散小波反变换 5、基于小波的信号去噪 6、基于小波的信号压缩 四、实验报告要求 按照《Matlab程序设计》模板提交实验报告 五、预习要求 预习课本上的相关内容 实验五 FFT 算法的应用 1.进一步加深对离散信号 DFT 的理解。 2.运用其 FFT 算法解决实际问题。 1.微机。 2.Matlab 编程环境。 1.熟悉 Matlab 的编程环境和编程语言。 2.学习教材 P213-227，P249-263，掌握快速傅里叶变换（FFT）的原理。 1.实验重点、难点、特点 快速傅里叶变换（FFT）的原理及应用。难点在 FFT 的应用以及 Matlab 编程中矩阵乘和数乘的应用。 2.实验原理 一、实验目的 二、实验仪器设备 三、实验学时 学时 四、预习要求 五、实验特点及实验原理简介 1.已知 2N 点实数序列 N=64。理论计算 X(k)=DFT[x(n)]2N。 六、实验内容及步骤 N=64;n=0:1:N-1;n2=0:1:2*N-1;k=0:1:2*N-1;xn1=cos(2*pi/N*7*2*n)+1/2*cos(2*pi/N*19*2*n);xn2=cos(2*pi/N*7*(2*n+1))+1/2*cos(2*pi/N*19*(2*n+1));wn=xn1+i*xn2;[wk,kk]=dft_my(wn,N);xk1=1/2*(wk+conj(wk(mod(-kk,N)+1)));xk2=1/2*(wk-conj(wk(mod(-kk,N)+1)))/i;XK1=fft(xn1,N);XK2=fft(xn2,N);%subplot(2,2,1);%stem(kk,abs(xk1));grid on;%subplot(2,2,3);%stem(kk,abs(XK1));grid on;%subplot(2,2,2);%stem(kk,abs(xk2));grid on;%subplot(2,2,4);%stem(kk,abs(XK2));grid on; xk=[xk1+(exp(-1i*pi/N).^kk).*xk2,xk1-(exp(-1i*pi/N).^kk).*xk2];subplot(2,1,1);stem(k,abs(xk),'r');grid on;xlabel('k');ylabel('2*NFFT');title('n点DFT完成2N点FFT');xn=cos(2*pi/N*7*n2)+1/2*cos(2*pi/N*19*n2);xkk=fft(xn,2*N);subplot(2,1,2);stem(k,abs(xkk),'g');grid on;xlabel('k');ylabel('2*NFFT');title('Matlab2N点FFT'); 七、问题思考 1.两次编程计算与理论计算相比较，结果一致吗？说明原因。 2.两次编程计算在编程方面、运算量上的比较。 八、心得体会 北京信息科技大学 数字信号处理实验报告 题 目： 数字信号处理课程设计实验 学 院： 信息与通信工程学院 专 业： 通信工程专业 姓 名： 班 级： 学 号： 指导老师： 实验目的 1、熟悉IIR数字滤波器的设计原理与方法。 2、掌握数字滤波器的计算机软件实现方法。 3、通过观察对实际心电图信号的滤波作用，学习数字滤波器在实际中的应用。 实验仪器及材料 计算机，MATLAB软件 实验内容及要求 1．设计巴特沃斯低通数字滤波器对人体心电信号进行滤波 (1)人体心电图信号在测量过程中会受到工业高频干扰，所以必须经过低通滤波处理，才能作为判断心脏功能的有用信息。以下为一个实际心电图信号采样序列x(n)，其中存在高频干扰，抽样周期Ts=1秒。在实验中，以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] 对序列x(n)用FFT做频谱分析，生成x(n)的频谱图。 (2)设计一个巴特沃斯低通IIR数字滤波器H(z)。 设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB； 在阻带内 [0.3π, π]频率区间上，最小衰减大于15dB。 j|H(e)|。写出数字滤波器H(z)的表达式，画出滤波器的幅频响应曲线 (3)用所设计的滤波器对实际心电图信号采样序列x(n)进行滤波处理，编写程序，求滤波后的序列y(n)，并分别画出滤波前后的心电图信号波形图和频谱图。 y(n)= [0,0,0,0, 0,0,0,0,-0.14025,0.40279,-0.56085 ,0.33328,0.023981,-0.18809,0.11843,-0.1038,0.11576,-0.1225,0.099815 ,-0.13769 ,0.095249,-0.0070273,0.018867,0.090543,-0.11257,-0.070884 ,0.17676,-0.55407,0.24813,-0.34732,-0.30428,0.59426,-0.29574,-0.063869,0.34018,-0.73334,1.0293,-0.57107,-0.2461,0.83605,-0.83026,0.45459,0.011551,-0.25667,0.23896,-0.17361,0.20829,-0.28417,0.28765 ,-0.2035,0.02865,0.066164,0.077916,-0.36052, 0.53517,-0.5571] 源程序 clear all，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];%未经滤波的心电图信号 L=length(x);l=0:L-1; y=fft(x,L);Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=15; [N,Wn] = buttord(Wp,Ws,Rp,Rs,'s');[b,a] = butter(N,Wn,'s');[numa,dena]=impinvar(b,a,1);w=linspace(0,pi,1024);h=freqz(numa,dena,w);norm=max(abs(h));numa=numa/norm;[z,p]=tf2zp(b,a);figure(1) plot(w,20*log10(abs(h)/norm));grid; xlabel('数字频率');ylabel('幅度响应dB');figure(2)plot(w,abs(h));grid; xlabel('数字频率'); ylabel('幅度响应|H(e^(jw))|');figure(3)zplane(z,p);xx=filter(b,a,x);yy=fft(xx,L);figure(4)subplot(2,1,1)stem(l,x); title('未经滤波的心电图信号');xlabel('n');subplot(2,1,2)stem(l,xx); title('经滤波之后的心电图信号');xlabel('n');figure(5)subplot(2,1,1)plot(l,abs(y)); title('未经滤波的心电图信号的频谱');subplot(2,1,2)plot(l,abs(yy)); title('经滤波处理的心电图信号的频谱'); 2．用help查看内部函数cheb1ord.m及cheby1.m，了解调用格式。 编程设计教材习题6-2，求模拟滤波器Ha(s)的表达式。 源程序 close all clear all clc Wp=2*pi*3000;Rp=2;Ws=2*pi*12000;Rs=50;Fs=24000; w=linspace(0,pi,1024); [N,Wn]=cheb1ord(Wp,Ws,Rp,Rs,'s');e=sqrt(10^(Rp/10)-1);[b,a]=cheby1(N,e,Wn,'s')[numa,dena]=impinvar(b,a,Fs);h=freqz(numa,dena,w);norm=max(abs(h)); plot(w*Fs/pi,20*log10(abs(h)/norm))title('幅度响应')xlabel('频率(Hz)')ylabel('幅度(dB)')grid 3．模拟滤波器的数字化 用内部函数impinvar及bilinear实现教材习题6-5，求数字滤波器H(z)的表达式。 源程序 close all clear all clc b1=[0 0 1];a1=[1 1 1];b2=[0 0 1];a2=[2 3 1]; [numa1,dena1]=impinvar(b1,a1,0.5)[numa2,dena2]=bilinear(b1,a1,0.5)[numa3,dena3]=impinvar(b2,a2,0.5) [numa4,dena4]=bilinear(b2,a2,0.5) 本实验所用的部分MATLAB函数 L=length(x)：求序列x长度。 y=fft(x,L)：将序列x(n)做L点快速傅立叶变换，结果赋给序列y(n)。 [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s')：计算模拟Butterworth滤波器的最小阶次n和截止频率为Wn。 [b,a] = butter(n,Wn,'s')：设计模拟截止频率为Wn（rad/s）的n阶 Butterworth低通滤波器，返回值为模拟滤波器的系数。 y=filter(b,a,x): 将序列x(n)通过滤波器滤波后生成序列y(n)，滤波器的分母多项式系数构成a向量，分子多项式系数构成b向量。 [BZ,AZ] = impinvar(B,A,Fs)：冲激响应不变法，返回值为数字滤波器的系数。 [BZ,AZ] = bilinear(B,A,fs)：双线性变换，返回值为数字滤波器的系数。 [H w]=freqz(b,a)：由滤波器分母多项式系数构成的a向量和分子多项式系数构成的b向量求系统频响。 截图 实验体会 这次的实验让这让我看到了理论与实践相结合的优势与用处，让我受益匪浅。我认识到了自己理论知识的不足，也认识到了我们学习的基础知识究竟能运用于什么领域，如何运用。我们在老师的耐心指导下调试电路，直到得到要求的效果。让我们在学习电路、信号等理论知识的同时，明白如何把这些应用于实际。第二篇：MATLAB实现数字信号处理
第三篇：数字信号处理实验讲稿
第四篇：数字信号处理实验5
第五篇：数字信号处理实验二
点击咨询 