DAY4.5音频信号处理

1、 滤波器 去除信号中不要的成分或者增强所需成分 模拟滤波器 电阻、电容、晶体管等组成的物理网络 数字滤波器 处理由模拟信号采样得到的数字信号 理论上可以实现任何滤波效果 寄存器、延时器、加法器和乘法器等基本 数字电路实现 September 21, 20122 数字滤波器的限制 成本 集成电路成本不断降低，数字滤波器越来越常见， 收音机、蜂窝电话等日常用品的重要组成部分 速度 不可能比滤波器内部的数字电路的运算速度更快 September 21, 20123 数字滤波器的特性-I 精度更高 例：很容易做到一个1000Hz的低通滤波器，允许 999Hz信号通过，并完全阻止1001Hz的信号，模

2、拟滤波器无法区分如此接近的信号 更高的信噪比 避免了模拟电路中噪声(如电阻热噪声)的影响 主要噪声源 在数字系统之前的模拟电路引入的电路噪声 数字系统输入端模数转换过程中产生的量化噪声 September 21, 20124 数字滤波器的特性-II 模拟滤波器不能比拟的可靠性 组成模拟滤波器的电子元件的电路特性会随着时 间、温度、电压的变化而漂移 数字电路没有这种问题。只要在数字电路的工作 环境下，数字滤波器就能够稳定可靠的工作 线性时不变 叠加原理(参见下页) 系统参数不随时间而变化 不管输入信号作用的时间先后，输出信号响应的形状均 相同，仅是从出现的时间不同 September 21, 2

3、0125 数字滤波器分类-I 线性与非线性 是否满足叠加原理 September 21, 20126 数字滤波器分类-II 因果系统 一个实际的物理系统，其当前时刻的输出只能和 当前时刻的输入、过去时刻的输入和输出有关， 而不能和将来时刻的输入输出有关 y(n)=fx(n),x(n-k),y(n-m) k0,m0 非因果系统 y(n)=fx(n+1),x(n+2),. September 21, 20127 数字滤波器分类-III 传递函数的Z变换表示 无限脉冲响应(IIR) M1，称为阶数，系统中反馈环的个数。反馈的 存在使IIR滤波器的脉冲响应为无限长 有限脉冲响应(FIR) A(z)=1

4、，系统的脉冲响应的长度为N+1 September 21, 20128 IIR滤波器的设计过程 根据参数要求设计对应的模拟滤波器(如巴 特沃斯滤波器、切比雪夫滤波器等) 通过映射(如脉冲响应不变法、双线性映射 等等)将模拟滤波器变换为数字滤波器 IIR滤波器的优点 直接利用模拟滤波器设计的成果 IIR滤波器的缺点 稳定性无法保证，数字运算可能溢出 September 21, 20129 FIR滤波器的优点 不存在系统极点，绝对稳定 确保了线性相位 不同频率信号分量经过FIR滤波器后时间差不变 不需要反馈，实现比IIR滤波器简单 FIR滤波器的缺点 性能不如同样阶数的IIR滤波器 优点远大于缺点

5、，比IIR滤波器的应用更广 September 21, 201210 FIR与IIR滤波器的比较 September 21, 201211 低通滤波器 容许低频信号通过，但减少频率高于截止 频率的信号的通过 绘制长期走势或均化 高通滤波器 容许高频信号通过，但减少频率低于截止 频率信号的通过 强调细节 September 21, 201212 带通滤波器 通过某一频率范围内的频率分量，但将其 他范围的频率分量衰减到极低水平 带阻滤波器 通过大多数频率分量，但将某些范围的频 率分量衰减到极低水平的滤波器 September 21, 201213 September 21, 201214 理想的滤

6、波器 通带完全平坦，没有增益/衰减 通带外所有频率被完全衰减掉 通带外的转换在极小的频率范围完成 理论上在频域中用信号乘以矩形函数得到， 等价于在时域与sinc(x)=sin(x)/x函数作 卷积 September 21, 201215 实际的滤波器-I sinc函数延伸到无穷远 为执行卷积需预测未来并记录过去的数据 预先录制好的数字信号可在信号后边补零 实时应用中的实际滤波器 将信号延时一小段时间让它们能够看到未来的一小部分 来近似地实现理想滤波器 近似精度越高所需要的延时越长 September 21, 201216 实际的滤波器-II 滚降 实际的滤波器并不能将期望频率范围外的所有频

7、率完全衰减掉，在所要的通带外还有一个被衰减 但是没有被隔离的范围 使用每十倍频的衰减幅度dB来表示 滚降范围越窄越好，性能就与设计更加接近 吉布斯现象 随着滚降范围越来越小，通带就变得不再平坦 开始出现波纹。在通带的边缘处尤其明显 September 21, 201217 September 21, 201218 波特图 September 21, 201219 频率响应 一阶滤波器 频率增加一倍，信号强度减弱一半(约-6dB) 截止频率之下是一条水平线，截止频率之上则是 一条斜线。边界为knee curve平缓转换 二阶滤波器 滚降速率更快，二阶的巴特沃斯滤波器频率增加 一倍，信号强度衰减到

8、最初的四分之一(-12dB) n阶滤波器 滚降速率是每倍频6ndB September 21, 201220 低通滤波器实例 固体屏障 当房间中播放音乐时，外边很容易听到音乐的低 音，但是高音部分大部分被过滤掉了 封闭的汽车中非常大的音乐声在另外一个车中的 人听来却是低音节拍 低通滤波器应用 音频设备中的hiss滤波器、平滑数据、图 像模糊处理 剔除短期波动、保留平滑的长期发展趋势 September 21, 201221 说明 低、高通滤波器仅是指滤波器响应的形状 高通滤波器的截止频率可能比任何低通滤 波器截止频率更低 任何所期望的频率范围可以到微波 微波的频率范围 米波：波长10m1m。频

9、率30MHz300MHz 分米波：波长10dm1dm。频率300MHz3000MHz 厘米波：波长10cm1cm。频率3000MHz30000MHz 毫米波：波长10mm1mm。频率3ooooMHz300000MHz September 21, 201222 常用滤波器 巴特沃斯响应 通带最平坦，接近DC信号 慢慢衰减至截止频率点为-3dB 逼近-20ndB/decade的衰减率，n为滤波器阶数 特别适用于低频应用，维护增益的平坦性 切比雪夫响应 通带具有一些纹波，更陡的衰减特性 对于元件的变化最不敏感，推荐使用 September 21, 201223 滤波器阶数 滤波器传递函数中的极点个数

10、 决定了转折区的下降速度 一般每增加一阶(一个极点)，就会增加一 20dBDec(一20dB每十倍频程) 滤波器技术指标 中心频率f0，即工作频带的中心 带宽BW 通带衰减，即通带内的最大衰减 阻带衰减 September 21, 201224 滤波器常用参数 截止频率Fc 响应曲线在通带内下降到误差带以外的频率点 在巴特沃斯滤波器中被称作3dB点 通带纹波Amax 通带内响应的起伏 最小阻带衰减Amin 阻带内最小信号衰减 September 21, 201225 Matlab滤波器设计举例 % Fdatool；% GUI Filter design tool % Chebyshev Typ

11、e II Lowpass filter，Fs = 44100， Fpass = 1000，Fstop = 1100，Apass = 1，Astop = 80， match = stopband % function Hd = untitled %UNTITLED Returns a discrete-time filter object. % M-File generated by MATLAB(R) 7.9 and the Signal % Processing Toolbox 6.12. % Generated on: 20-Apr-2013 16:49:10 % Elliptic Low

12、pass filter designed using September 21, 201226 % FDESIGN.LOWPASS. % All frequency values are in Hz. Fs = 48000; % Sampling Frequency Fpass = 3500; % Passband Frequency Fstop = 4000; % Stopband Frequency Apass = 1; % Passband Ripple (dB) Astop = 80; % Stopband Attenuation (dB) match = both; % Band t

13、o match exactly % Construct an FDESIGN object and call its ELLIP method. h = fdesign.lowpass(Fpass, Fstop, Apass, Astop, Fs); September 21, 201227 Hd = design(h, ellip, MatchExactly, match); % EOF % x=wavread(RockOriginal.wav); x=wavread(PianoOriginal.wav); y = filter(Hd,x); wavplay(x,44100); wavpla

14、y(y,44100); September 21, 201228 傅里叶变换-I 历史 傅立叶：法国数学家和物理学家 对热传递很感兴趣，1807年在法国科学学会上发 表一篇论文，运用正弦曲线描述温度分布 当时具有争议性的决断 任何连续周期信号可以由一组正弦曲线组合而成，简单 正弦曲线保真度性质 输入的正弦曲线信号变换后输出仍是正弦曲线 只有幅度和相位可能发生变化 频率和波的形状是一样的 September 21, 201229 傅里叶变换-II 文章发表的故事 审查者：数学家拉格朗日(1736-1813)，拉普拉 斯(1749-1827). 拉格朗日坚决反对 认为傅立叶的方法无法表示带有棱角的

15、信号，如在方波 中出现非连续变化斜率 法国科学学会屈服于拉格朗日的威望，拒绝了傅 立叶的工作，直到拉格朗日死后15年论文才发表 参加了政治运动，随拿破仑远征埃及，法国大革命后因 会被推上断头台而一直在逃避 不能精确表示，但无限接近，都对 September 21, 201230 傅里叶变换-III 线性积分变换 将满足一定条件的某个函数表示成三角函数(正 弦和/或余弦函数)或者它们的积分的线性组合 应用广泛 物理学、声学、光学、结构动力学、数论、组合 数学、概率论、统计学、信号处理、密码学、海 洋学、通讯 将信号分解成振幅分量和频率分量 September 21, 201231 连续傅里叶变换

16、 离散傅里叶变换 September 21, 201232 傅里叶变换的实部与虚部 September 21, 201233 傅里叶变换的Z变换表示 September 21, 201234 傅里叶变换的性质 线性性质 September 21, 201235 傅里叶变换的性质 反转性 September 21, 201236 傅里叶变换的性质 尺度变换性质 September 21, 201237 September 21, 201238 September 21, 201239 傅里叶变换的性质 时移特性 September 21, 201240 卷积定理 September 21, 20

17、1241 帕塞瓦尔定理 September 21, 201242 Introduction DFTs with a million points are common DFT (Discrete Fourier Transform) computational complexity O(n2) FFT (Fast Fourier Transform) computational complexity O(n log n) efficient method for computing the DFT Example For n = 220, almost 35,000 faster than D

18、FT September 21, 201243 Window length the length of the input data vector Transform length the length of the output, the computed DFT FFT pads or chops the input to achieve the desired transform length September 21, 201244 Execution time of FFT depends on the transform length fastest when the transf

19、orm length is a power of two almost as fast when the transform length has only small prime factors typically slower for transform lengths that are prime or have large prime factors in practice time differences，insignificant in modern FFT algorithms in Matlab， unnecessary to adjust the transform leng

20、th for efficiency September 21, 201245 y = fft(x) window length m = length(x) transform length n = length(y) m=n y = fft(x,n) transform length is n If length(x)n, only the first n elements of x are used to compute the transform September 21, 201246 September 21, 201247 Example-1 x: with two componen

21、t frequencies of differing amplitude and phase buried in noise fs = 100; % Sample frequency (Hz) t = 0:1/fs:10-1/fs; % 10 sec sample x = (1.3)*sin(2*pi*15*t) . % 15 Hz component + (1.7)*sin(2*pi*40*(t-2) . % 40 Hz component + (2.5)*randn(size(t); % Gaussian noise m = length(x); % Window length, 1000

22、 n = pow2(nextpow2(m); % Transform length, 1024 % p=nextpow2(m); 2p =abs(m) % n=pow2(p)=2p % n = m; y = fft(x,n); % DFT, 1024 f = (0:n-1)*(fs/n); % Frequency range, 1024 power = y.*conj(y)/n; % Power of the DFT % conj(Z) = real(Z) - i*imag(Z) plot(f,power) xlabel(Frequency (Hz) ylabel(Power) title(b

23、f Periodogram) September 21, 201248 The first half of the frequency range (from 0 to the Nyquist frequency fs/2) is sufficient to identify the component frequencies the second half is just a reflection of the first half September 21, 201249 Example-II Spectral Analysis of a Whale Call Because blue w

24、hale calls are so low, they are barely audible to humans. The time scale in the data is compressed by a factor of 10 to raise the pitch and make the call more clearly audible cd C:Program FilesMATLABR2009bhelptechdocmathexamples; x,fs = auread(bluewhale.au) % fs=4000Hz plot(x); xlabel(Sample Number)

25、; ylabel(Amplitude); title(bf Blue Whale Call); sound(x,fs); September 21, 201250 An A trill is followed by a series of B moans. September 21, 201251 B call is simpler and easier to analyze. Use the previous plot to determine approximate indices for the beginning and end of the first B call. Correct the time base for the factor of 10 speed-up in the data: bCall = x(2.45e4:3.10e4); % amp