DAY4.5音频信号处理

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

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

3、非线性 是否满足叠加原理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，系统的脉冲响应的长度为N+1September

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

5、IR与IIR滤波器的比较September 21, 201211低通滤波器 容许低频信号通过，但减少频率高于截止频率的信号的通过 绘制长期走势或均化高通滤波器 容许高频信号通过，但减少频率低于截止频率信号的通过 强调细节September 21, 201212带通滤波器 通过某一频率范围内的频率分量，但将其他范围的频率分量衰减到极低水平带阻滤波器 通过大多数频率分量，但将某些范围的频率分量衰减到极低水平的滤波器September 21, 201213September 21, 201214理想的滤波器 通带完全平坦，没有增益/衰减 通带外所有频率被完全衰减掉 通带外的转换在极小的频率范围完成

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

7、，性能就与设计更加接近 吉布斯现象 随着滚降范围越来越小，通带就变得不再平坦开始出现波纹。在通带的边缘处尤其明显September 21, 201217September 21, 201218波特图September 21, 201219频率响应 一阶滤波器 频率增加一倍，信号强度减弱一半(约-6dB) 截止频率之下是一条水平线，截止频率之上则是一条斜线。边界为knee curve平缓转换 二阶滤波器 滚降速率更快，二阶的巴特沃斯滤波器频率增加一倍，信号强度衰减到最初的四分之一(-12dB) n阶滤波器 滚降速率是每倍频6ndBSeptember 21, 201220低通滤波器实例 固体屏障

8、当房间中播放音乐时，外边很容易听到音乐的低音，但是高音部分大部分被过滤掉了 封闭的汽车中非常大的音乐声在另外一个车中的人听来却是低音节拍低通滤波器应用 音频设备中的hiss滤波器、平滑数据、图像模糊处理 剔除短期波动、保留平滑的长期发展趋势September 21, 201221说明 低、高通滤波器仅是指滤波器响应的形状 高通滤波器的截止频率可能比任何低通滤波器截止频率更低 任何所期望的频率范围可以到微波 微波的频率范围 米波：波长10m1m。频率30MHz300MHz 分米波：波长10dm1dm。频率300MHz3000MHz 厘米波：波长10cm1cm。频率3000MHz30000MHz

9、毫米波：波长10mm1mm。频率3ooooMHz300000MHzSeptember 21, 201222常用滤波器 巴特沃斯响应 通带最平坦，接近DC信号 慢慢衰减至截止频率点为-3dB 逼近-20ndB/decade的衰减率，n为滤波器阶数 特别适用于低频应用，维护增益的平坦性 切比雪夫响应 通带具有一些纹波，更陡的衰减特性 对于元件的变化最不敏感，推荐使用September 21, 201223滤波器阶数 滤波器传递函数中的极点个数 决定了转折区的下降速度 一般每增加一阶(一个极点)，就会增加一20dBDec(一20dB每十倍频程) 滤波器技术指标 中心频率f0，即工作频带的中心 带宽B

10、W 通带衰减，即通带内的最大衰减 阻带衰减September 21, 201224滤波器常用参数 截止频率Fc 响应曲线在通带内下降到误差带以外的频率点 在巴特沃斯滤波器中被称作3dB点 通带纹波Amax 通带内响应的起伏 最小阻带衰减Amin 阻带内最小信号衰减September 21, 201225Matlab滤波器设计举例 % Fdatool；% GUI Filter design tool % Chebyshev Type II Lowpass filter，Fs = 44100，Fpass = 1000，Fstop = 1100，Apass = 1，Astop = 80，match

11、= 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 Lowpass filter designed usingSeptember 21, 201226% FDESIGN.LOWPASS.% All frequency values

12、 are in Hz.Fs = 48000; % Sampling FrequencyFpass = 3500; % Passband FrequencyFstop = 4000; % Stopband FrequencyApass = 1; % Passband Ripple (dB)Astop = 80; % Stopband Attenuation (dB)match = both; % Band to match exactly% Construct an FDESIGN object and call its ELLIP method.h = fdesign.lowpass(Fpas

13、s, Fstop, Apass, Astop, Fs);September 21, 201227Hd = design(h, ellip, MatchExactly, match);% EOF% x=wavread(RockOriginal.wav);x=wavread(PianoOriginal.wav);y = filter(Hd,x);wavplay(x,44100);wavplay(y,44100);September 21, 201228傅里叶变换-I 历史 傅立叶：法国数学家和物理学家 对热传递很感兴趣，1807年在法国科学学会上发表一篇论文，运用正弦曲线描述温度分布 当时具有争议

14、性的决断 任何连续周期信号可以由一组正弦曲线组合而成，简单 正弦曲线保真度性质 输入的正弦曲线信号变换后输出仍是正弦曲线 只有幅度和相位可能发生变化 频率和波的形状是一样的September 21, 201229傅里叶变换-II 文章发表的故事 审查者：数学家拉格朗日(1736-1813)，拉普拉斯(1749-1827). 拉格朗日坚决反对 认为傅立叶的方法无法表示带有棱角的信号，如在方波中出现非连续变化斜率 法国科学学会屈服于拉格朗日的威望，拒绝了傅立叶的工作，直到拉格朗日死后15年论文才发表 参加了政治运动，随拿破仑远征埃及，法国大革命后因会被推上断头台而一直在逃避 不能精确表示，但无限接

15、近，都对September 21, 201230傅里叶变换-III 线性积分变换 将满足一定条件的某个函数表示成三角函数(正弦和/或余弦函数)或者它们的积分的线性组合 应用广泛 物理学、声学、光学、结构动力学、数论、组合数学、概率论、统计学、信号处理、密码学、海洋学、通讯 将信号分解成振幅分量和频率分量September 21, 201231连续傅里叶变换离散傅里叶变换September 21, 201232傅里叶变换的实部与虚部September 21, 201233傅里叶变换的Z变换表示September 21, 201234傅里叶变换的性质 线性性质September 21, 20123

16、5傅里叶变换的性质 反转性September 21, 201236傅里叶变换的性质 尺度变换性质September 21, 201237September 21, 201238September 21, 201239傅里叶变换的性质 时移特性September 21, 201240卷积定理September 21, 201241帕塞瓦尔定理September 21, 201242Introduction DFTs with a million points are commonDFT (Discrete Fourier Transform) computational complexity O

17、(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 DFTSeptember 21, 201243Window length the length of the input data vector Transform length the length of the output, the computed DFT FFT pads or

18、chops the input to achieve the desired transform lengthSeptember 21, 201244Execution time of FFT depends on the transform length fastest when the transform 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 p

19、rime or have large prime factors in practice time differences，insignificant in modern FFT algorithms in Matlab， unnecessary to adjust the transform length for efficiencySeptember 21, 201245y = fft(x) window length m = length(x) transform length n = length(y) m=ny = fft(x,n) transform length is n If

20、length(x)n, only the first n elements of x are used to compute the transformSeptember 21, 201246September 21, 201247Example-1 x: with two component 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

21、*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 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); % Freque

22、ncy 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(bf Periodogram)September 21, 201248The first half of the frequency range (from 0 to the Nyquist frequency fs/2) is sufficient to identify the component fre

23、quencies the second half is just a reflection of the first halfSeptember 21, 201249Example-II Spectral Analysis of a Whale Call Because blue whale 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

24、 clearly audible cd C:Program FilesMATLABR2009bhelptechdocmathexamples; x,fs = auread(bluewhale.au) % fs=4000Hz plot(x); xlabel(Sample Number); 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); % amplitude tb = 1