林数字信号处理课程设计

1 广西工学院 数字信号处理 课程设计课程设计 说说 明明 书书 设计题目设计题目 语音信号的处理与滤波语音信号的处理与滤波 系 别 计算机工程系 专业班级 通信 062 班 学生姓名 林芳慧 学 号 200600402059 指导教师 周坚和 日 期 2009 年 1 月 2 日 2 目录目录 一 设计目的 3 二 设计内容 3 三 设计原理 4 四 设计步骤与程序 5 五 设计思考 19 六 结论 19 七 参考文献 20 3 数字信号处理课程设计数字信号处理课程设计 一 设计目的一 设计目的 这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波 器设计 通过理论推导得出相应结论 并利用 MATLAB 或者 开发系统作为 工具进行实现 从而复习巩固课堂所学的理论知识 提高对所学知识的综合应 用能力 并从实践上初步实现对数字信号的处理 二 设计内容二 设计内容 语音信号的处理与滤波 1 熟悉并掌握 MATLAB 中有关声音 wave 录制 播放 存储和读取的函数 在 MATLAB 环境中 使用声音相关函数录制 2 秒左右自己的声音 抽样率是 8000Hz s 并考虑如何解决一个实际问题 录制刚开始时 常会出现实际发出 声音落后录制动作半拍的现象 如何排除对这些无效点的采样 2 语音信号的频谱分析 分别取 8000 个和 16000 个数据进行频谱分析 得到幅度和相位谱 比较二 者异同并分析原因 要求首先画出语音信号的时域波形 然后对语音信号进行 频谱分析 在 MATLAB 中 可以利用函数 fft 对信号进行快速傅立叶变换 得到 信号的频谱特性 从而加深对频谱特性的理解 3 设计数字滤波器和画出频率响应 针对电话信道 最高 3500Hz 设计一个 FIR 或 IIR 滤波器进行滤波 把 抽样率转变为 7000Hz s 并进行频谱分析 得到幅度和相位谱 根据语音信号 的特点给出有关滤波器的性能指标 1 低通滤波器性能指标 fp 1000Hz fc 1200 Hz As 80dB Ap 1dB 2 高通滤波器性能指标 fc 2800 Hz fp 3000 Hz As 80dB Ap 1dB 3 带通滤波器性能指标 fp1 1200 Hz fp2 3000 Hz fc1 1000 Hz fc2 3200 Hz As 80dB Ap 1dB 要求学生首先用窗函数法设计上面要求的三种滤波器 在 MATLAB 中 可以利用 函数 fir1 设计 FIR 滤波器 然后在用双线性变换法设计上面要求的三种滤波器 之后再利用函数 butter 和 cheby1 设计上面要求的三种 IIR 滤波器 最后 利 用 MATLAB 中的函数 freqz 画出各滤波器的频率响应 4 4 用滤波器对信号进行滤波 比较 FIR 和 IIR 两种滤波器的性能 然后用性能好的各滤波器分别对采集 的信号进行滤波 在 MATLAB 中 FIR 滤波器利用函数 fftfilt 对信号进行滤波 IIR 滤波器利用函数 filter 对信号进行滤波 5 比较滤波前后语音信号的波形及频谱 要求在一个窗口同时画出滤波前后的波形及频谱 6 回放语音信号 在 MATLAB 中 函数 sound 可以对声音进行回放 其调用格式 sound x fs bits 可以感觉滤波前后的声音有变化 三 设计原理三 设计原理 1 IIR 滤波器双线性变换法设计原理 双线性变换法中 S 平面与 z 平面之间 满足以下映射关系 s 平面的虚轴单值地映射于 z 平面的单位圆上 s 平面的左半平面完全映射到 z 平面的单位圆内 双线性变换不存在混叠问题 双线性变换时一种非线性变换 这种非线性引起的幅频特性畸变可通过预畸而得到校正 双线性法设计 IIR 数字滤波器的步骤 1 将数字滤波器的频率指标 由 k Wk 2 T tan 转换为模拟滤波器的频率指标 2 根据通带临 k w k 界频率和阻带临界频率计算模拟低通原型滤波器的阶数 N 并求得低通原型的 传递函数 H s 3 用上面的双线性变换公式代入 H s 求出所设计的传递 函数 H z 4 分析滤波器特性 检查其指标是否满足要求 2 窗函数法设计线性相位 FIR 滤波器步骤 1 确定数字滤波器的性能要求 临界频率 k 滤波器单位脉冲响应长度 N 根据性能要求 合理选择单位脉冲响应 h n 的奇偶对称性 从而确定理想频率 响应 Hd ej 的幅频特性和相频特性 2 求理想单位脉冲响应 hd n 在实际计算中 可对 Hd ej 按 M M 远大于 N 点等距离采样 并对其求 IDFT 得 hM n 用 hM n 代替 hd n 3 选择适当的窗函数 w n 根据过渡带宽及阻带衰减要求 选择窗函数的类 型并估计窗口长度 N 或阶数 M N 1 窗函数类型可根据最小阻带衰减 As 独立 选择 因为窗口长度 N 对最小阻带衰减 As 没有影响 在确定窗函数类型以后 可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度 N 设待求滤波 5 器的过渡带宽为 w 它与窗口长度 N 近似成反比 窗函数类型确定后 其计 算公式也确定了 不过这些公式是近似的 得出的窗口长度还要在计算中逐步 修正 原则是在保证阻带衰减满足要求的情况下 尽量选择较小的 N 在 N 和 窗函数类型确定后 即可调用 MATLAB 中的窗函数求出窗函数 wd n 根据 h n hd n w n 求所需设计的 FIR 滤波器单位脉冲响应 求 H ej 分析其幅频 特性 若不满足要求 可适当改变窗函数形式或长度 N 重复上述设计过程 以得到满意的结果 窗函数的傅式变换 W ej 的主瓣决定了 H ej 过渡带宽 W ej 的旁瓣大小和多少决定了 H ej 在通带和阻带范围内波动幅度 四 课程设计的步骤四 课程设计的步骤 1 语音信号的采集 1 在 MATLAB 软件平台下 利用函数 wavread 对语音信号进行采样 记住采 样频率和采样点数 采样频率即取样频率 指每秒钟取得声音样本的次数 它的采样频率越高 声音的质量也就越好 但人耳的分辨率很有限 所以太高的频率就分辨不出好坏 来 对于高于 48KHz 的采样频率人耳已无法辨别出来了 由采样定理知 采样 频率要大于信号最高频率的两倍 否则会产生混叠 人的声音频率在 300Hz 3400Hz 所以取采样频率为 8000Hz 这次实验是利用 MATLAB 来录制的 采样位数 可以理解为声卡处理声音的解析度 这个数值越大 解析度就越高 录制和回 放的声音就越真实 电脑中的声音文件是用数字 0 和 1 来表示的 所以在电脑 上录音的本质就是把模拟声音信号转换成数字信号 反之 在播放时则是把数 字信号还原成模拟声音信号输出 实验中 wavread 函数返回的 Nbits 为 16 位 声音效果还是很好的 我的录音程序 fs 8000 采样频率 fprintf 开始录音 y wavrecord 2 fs fs double 录制声音 2 秒 fprintf 录音结束 wavplay y fs 播放录音 wavwrite y fs d shengyin 声音的存储 wavread d shengyin 读取声音 figure 1 plot abs y 画出声音采样后的时域波形 grid on 添加网格 title 语音信号的时域波形 给图形加注标签 z fft y 8192 对信号做 8192 点 FFT 变换 f fs 0 4095 8192 采样点数 figure 2 plot f abs z 1 4096 绘制原始语音信号的幅频响应图 axis 0 4000 0 20 title FFT 后声音的频谱特性 给图形加注标签 grid on 添加网格 6 录制刚开始时 常会出现实际发出声音落后录制动作半拍的现象 现在排除 对这些无效点 fs 8000 size wavread d shengyin size 声音的数据长度 z wavread d shengyin 2000 16000 读取声音 wavplay z fs 播放声音 wavwrite z fs d shengyin1 声音的存储 wavread d shengyin1 读取声音 figure 3 plot abs z 画出声音采样后的时域波形 grid on 添加网格 title 排除无效点后语音信号的时域波形 给图形加注标签说明 z fft y 8192 对信号做 fs 点 FFT 变换 figure 4 f fs 0 4095 8192 采样点数 plot f abs z 1 4096 绘制语音信号的幅频响应图 axis 0 4000 0 20 title 排除无效点经 FFT 后声音的频谱特性 给图形加注标签说明 grid on 添加网格 0200040006000800010000120001400016000 0 0 01 0 02 0 03 0 04 0 05 0 06 个 个 个 个 个 个 个 个 个 050001000015000 0 0 01 0 02 0 03 0 04 0 05 0 06 个 个 个 个 个 个 个 个 个 个 个 个 个 个 个 7 05001000150020002500300035004000 0 2 4 6 8 10 12 14 16 18 20 FFT个 个 个 个 个 个 个 个 05001000150020002500300035004000 0 2 4 6 8 10 12 14 16 18 20 个 个 个 个 个 个 FFT个 个 个 个 个 个 个 个 2 取 8000 个和 16000 个数据进行频谱分析 得到幅度和相位谱 比较二者异 同并分析原因 程序如下 x1 wavread D shengyin 8000 读取语音信号的数据 赋给变量 x1 sound x1 8000 播放语音信号 x2 wavread D shengyin 16000 读取语音信号的数据 赋给变量 x2 sound x1 8000 播放语音信号 y1 fft x1 8192 对信号做 8192 点 FFT 变换 y2 fft x2 8192 f fs 0 4095 8192 figure 1 subplot 1 2 1 plot x1 title 取 8000 数据的语音信号 xlabel 时间 ylabel 幅值 subplot 1 2 2 plot x2 title 取 16000 数据的语音信号 xlabel 时间 ylabel 幅值 figure 2 8 subplot 1 2 1 plot angle y1 title 相位图 subplot 1 2 2 plot angle y2 title 相位图 figure 3 subplot 1 2 1 plot f abs y1 1 4096 axis 0 4000 0 20 title 语音信号频谱 xlabel 频率 ylabel 幅值 subplot 1 2 2 plot f abs y2 1 4096 axis 0 4000 0 20 title 语音信号频谱 xlabel 频率 ylabel 幅值 0500010000 0 0 005 0 01 0 015 0 02 0 025 0 03 0 035 0 04 0 045 0 05 个 8000个 个 个 个 个 个 个 个 个 个 个 0500010000 0 0 005 0 01 0 015 0 02 0 025 0 03 0 035 0 04 0 045 0 05 个 16000个 个 个 个 个 个 个 个 个 个 个 0500010000 4 3 2 1 0 1 2 3 4 个 个 个 0500010000 4 3 2 1 0 1 2 3 4 个 个 个 9 0100200300400 0 2 4 6 8 10 12 14 16 18 20 8000个 个 个 个 个 个 个 个 个 个 个 个 个 0100200300400 0 2 4 6 8 10 12 14 16 18 20 16000个 个 个 个 个 个 个 个 个 个 个 个 个 分析异同点 从两者的时域图看出取 8000 点数据的语音信号被截断了 而取 16000 点数据的 声音则完好地呈现出来 从频谱图看出幅值大小与 fft 选用点数 N 有关 N 越 大 幅值就相对大一点 3 设计数字滤波器和对信号滤波 窗函数法设计 FIR 滤波器 1 窗函数设计低通滤波器 程序设计如下 clear fs 8000 采样频率 fprintf 开始录音 y wavrecord 2 fs fs double 录制声音 2 秒 fprintf 录音结束 wavplay y fs 播放录音 wavwrite y fs d shengyin3 声音的存储 z1 fs bits wavread d shengyin3 Y1 fft z1 8192 fp 1000 fc 1200 As 100 Ap 1 Fs 7000 wc 2 pi fc Fs wp 2 pi fp Fs wdel wc wp beta 10 056 根据阻带衰减选凯泽窗 N0 ceil 6 4 2 pi wdel N N0 mod N0 1 2 确保 N 为奇数 wn kaiser N 1 beta ws wp wc 2 pi b fir1 N ws wn figure 1 freqz b 1 x fftfilt b z1 X fft x 8192 10 figure 2 f fs 0 4095 8192 subplot 3 2 1 plot f abs Y1 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波前信号频谱 subplot 3 2 2 plot f abs X 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波后信号频谱 subplot 3 2 3 plot angle Y1 title 滤波前相位谱 subplot 3 2 4 plot angle X title 滤波后相位谱 subplot 3 2 5 plot z1 title 滤波前信号波形 subplot 3 2 6 plot x title 滤波后信号波形 sound x fs bits wavwrite x fs d shengyin4 wav 声音的存储 11 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 00 511 52 x 10 4 0 1 0 0 1 个 个 个 个 个 个 个 00 511 52 x 10 4 0 1 0 0 1 个 个 个 个 个 个 个 可以看出 FIR 低通滤波器把 1000 赫兹以上的频率成分基本滤去 达到信号处理 效果 在实际系统中 一个输入信号可以分解为多个正弦信号的叠加 为了 使得输出信号不会产生相位失真 必须要求它所包含的这些正弦信号通过系 统的时间是一样的 从系统的频率响应来看 就是要求它的相频特性是一条 直线 在 FIR 滤波器的线性相位的性质 这就 避免相位失真 2 窗函数设计高通滤波器 程序设计如下 clear close all wavwrite y fs d shengyin3 Wav 声音的存储 z1 fs bits wavread d shengyin3 读取声音 Y1 fft z1 8192 fp 2800 fc 3000 As 100 Ap 1 Fs 7000 wc 2 pi fc Fs wp 2 pi fp Fs wdel wc wp beta 10 056 N1 ceil 6 4 2 pi wdel N N1 mod N1 1 2 高通滤波器要确保 N 为奇数 wn kaiser N beta ws wp wc 2 pi b fir1 N 1 ws high wn figure 1 12 freqz b 1 x fftfilt b z1 X fft x 8192 figure 2 f fs 0 4095 8192 subplot 3 2 1 plot f abs Y1 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波前信号频谱 subplot 3 2 2 plot f abs X 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波后信号频谱 subplot 3 2 3 plot angle Y1 title 滤波前相位谱 subplot 3 2 4 plot angle X title 滤波后相位谱 subplot 3 2 5 plot z1 title 滤波前信号波形 subplot 3 2 6 plot x title 滤波后信号波形 sound x fs bits wavwrite x fs d shengyin5 wav 声音的存储 图形分析如下 13 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 00 511 52 x 10 4 0 1 0 0 1 个 个 个 个 个 个 个 00 511 52 x 10 4 0 01 0 0 01 个 个 个 个 个 个 个 3 窗函数设计带通滤波器 z1 fs bits wavread d shengyin3 Y1 fft z1 8192 fp1 1200 fp2 3000 fc1 1000 fc2 3200 As 100 Ap 1 Fs 7000 wp1 2 pi fp1 Fs wc1 2 pi fc1 Fs wp2 2 pi fp2 Fs wc2 2 pi fc2 Fs wdel wp1 wc1 beta 10 056 N1 ceil 6 4 2 pi wdel N N1 mod N1 1 2 确保 N 为奇数 ws wp1 wc1 2 pi wp2 wc2 2 pi wn kaiser N 1 beta b fir1 N ws wn figure 1 freqz b 1 x fftfilt b z1 X fft x 8192 f fs 0 4095 8192 figure 2 f fs 0 4095 8192 subplot 3 2 1 plot f abs Y1 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 14 title 滤波前信号频谱 subplot 3 2 2 plot f abs X 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波后信号频谱 subplot 3 2 3 plot angle Y1 title 滤波前相位谱 subplot 3 2 4 plot angle X title 滤波后相位谱 subplot 3 2 5 plot z1 title 滤波前信号波形 subplot 3 2 6 plot x title 滤波后信号波形 sound x fs bits wavwrite x fs d shengyin8 wav 00 10 20 30 40 50 60 70 80 91 15000 10000 5000 0 5000 Normalized Frequency rad sample Phase degrees 00 10 20 30 40 50 60 70 80 91 200 100 0 100 Normalized Frequency rad sample Magnitude dB 15 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 00 511 52 x 104 0 1 0 0 1 个 个 个 个 个 个 个 00 511 52 x 104 0 1 0 0 1 个 个 个 个 个 个 个 双线性变换法设计 IIR 低通滤波器 选用 butter 程序设计如下 clear close all z1 fs bits wavread d shengyin3 Y1 fft z1 8192 fp 1000 fc 1200 As 100 Ap 1 Fs 7000 wc 2 fc Fs wp 2 fp Fs N ws buttord wc wp Ap As b a butter N ws figure 1 h w freqz b a plot w pi abs h grid xlabel omega pi ylabel 振幅 幅值 title Butterworth 型数字低通滤波器的幅频响应 x filter b a z1 X fft x 8192 figure 2 f fs 0 4095 8192 subplot 3 2 1 plot f abs Y1 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 16 title 滤波前信号频谱 subplot 3 2 2 plot f abs X 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波后信号频谱 subplot 3 2 3 plot angle Y1 title 滤波前相位谱 subplot 3 2 4 plot angle X title 滤波后相位谱 subplot 3 2 5 plot z1 title 滤波前信号波形 subplot 3 2 6 plot x title 滤波后信号波形 sound x fs bits wavwrite x fs d shengyin6 wav 图形分析 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 00 511 52 x 104 0 1 0 0 1 个 个 个 个 个 个 个 00 511 52 x 104 0 1 0 0 1 个 个 个 个 个 个 个 17 5 双线性变换法设计高通滤波器 选用 cheby1 程序设计如下 z1 fs bits wavread d shengyin3 Y1 fft z1 8192 fp 3000 fc 2800 As 100 Ap 1 Fs 7000 wc 2 fc Fs wp 2 fp Fs N ws cheb1ord wc wp Ap As b a cheby1 N Ap wp high figure 1 h w freqz b a plot w pi abs h grid xlabel omega pi ylabel 振幅 幅值 title 切比雪夫型数字高通滤波器的幅频响应 x filter b a z1 X fft x 8192 f fs 0 4095 8192 figure 2 f fs 0 4095 8192 subplot 3 2 1 plot f abs Y1 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波前信号频谱 subplot 3 2 2 plot f abs X 1 4096 axis 0 4000 0 10 xlabel 频率 ylabel 振幅 title 滤波后信号频谱 subplot 3 2 3 plot angle Y1 title 滤波前相位谱 subplot 3 2 4 plot angle X title 滤波后相位谱 subplot 3 2 5 plot z1 title 滤波前信号波形 subplot 3 2 6 plot x title 滤波后信号波形 sound x fs bits 18 wavwrite x fs d shengyin7 wav 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 01000200030004000 0 5 10 个 个 个 个 个 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 0500010000 5 0 5 个 个 个 个 个 个 00 511 52 x 104 0 1 0 0 1 个 个 个 个 个 个 个 00 511 52 x 104 0 01 0 0 01 个 个 个 个 个 个 个 以上几种滤波 我们都可以从信号滤波前后的波形图以及频谱图上看出变化 当然 也可以用 sound 函数来播放滤波后的语音 从听觉上直接感受语音信 号的变化 但由于人耳听力的限制 有些情况下我们很难听出来的 五 设计报告思考 1 在双线性变换法中 模拟频率与数字频率不再是线性关系 所以一个线性相 位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了 如 以上实验过程中 采用双线性变化法设计的 butter 和 cheby1 数字滤波器 从 图中可以看到这种非线性关系 而窗函数设计法则保持原有的线性相位 2 分析实际信号频谱过程中 分别用 N 8092 点 N 4096 来分析幅度谱 由于采 样频率 f 8000Hz 故 Nyquist 频率为 4000Hz 整个频谱图是以 Nyquist 频率 为轴对称的 因此利用 fft 对信号作频谱分析 只要考察 0 Nyquist 频率 采 样频率一半 范围的幅频特性 比较可见 幅值大小与 fft 选用点数