多采样应用——语音信号的降采样与重构

1语音信号的降采样处理和插值重构1.引言本文是根据网上找的“信号的分析与处理综合实验”的内容，通过学习和MATLAB实践后的学习总结。实验内容为真实语音信号的采样重构，具体要求如下：录制一段自己的语音信号，并对录制的信号进行采样；画出采样前后语音信号的时域波形和频谱图；对降采样后的信号进行插值重构，滤波，恢复原信号。2.降采样原理对采样数据每隔M-1个点（M为整数）取一个，结果使得在时间间隔里的采样数据被压缩M倍，如图所示。设数据采样率为fs，奈奎斯特频率为fs/2，则降采样后数据采样率为fs/M，其奈奎斯特频率为fs/2M。采样率降低后，它相对于原始输入信号的频带就成了欠采样，这样可能产生混迭，为避免混迭，输入的采样数据必须先进行反混叠滤波（这个反混叠滤波不是系统的模拟反混叠滤波器，是针对采样数据的数字化反混叠滤波器），使它所含最多频率分量低于降采样后的奈奎斯特频率fs/2M。反混叠滤波和降采样合称为抽取器。抽取器在时域里的输入-输出关系表示为：kkymWmMhkxmMkWnhkxnk3.语音信号的降采样处理H0.wav文件是用xp系统自带的录音机功能用22KHz采样率采集的男声“零”的发音，长度为两秒。调用Matlab的wavread函数读取波形文件，比较不同降采样率后的声音，从波形图、频谱和听觉上对比。用sound函数播放时，采样声音要与采样频率对应才能不失真。ttx(n)y(n)TT23.1Matlab中波形读取与降采样函数y=wavread(file)，%读取file所规定的wav文件，返回采样值放在向量y中。[y,fs,nbits]=wavread(file)，%采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。y=x((1:N:length(x)));%对原始信号每隔N个点取一位，即采样率变为原来的1/Ny=resample(yn,L,M);%采样率变为原来的L/M倍y=downsample(yn,N);%%采样率变为原来的1/N倍3.2语音信号读取与降采样处理原始语音信号的最高频率分量在4kHz左右，采样率为22kHz，故不用加抗混叠滤波器。通过观察原始信号的频谱发现直流分量特别大，影响了其他分量的观察，故滤除直流分量后作为参考信号。对比如图3-1：图3-1信号波形图和频谱图改变采样率为原来的1/2倍,1/10倍,1/100倍，分别得到降采样后的信号波形和频谱图，如图3-2。可以发现，降低采样率后，频谱出现了混叠和泄漏，尤其是到1/10倍采样后，低频分量丢失，混叠和泄漏严重，因此在降采样为1/M倍时，若信号最高频率超过fs/2M，则需要加抗混叠滤波器。通过对比1/2倍和原始信号频谱可以发现，虽然1/2倍采样信号的采样频谱为11kHz，接近原始语音信号的最高频率分量4kHz的两倍，满足Nyquist定理，但也有明显的频谱混叠和泄漏。所以采样率尽量取信号最高频率的4~10倍，这样得到的频谱更准确。用Matlab的sound函数播放声音，听觉上1/2倍采样信号和原始信号无太大差别，有意思的是原始信号有沙沙的噪音，而1/2倍采样信号却没有了。1/10倍,1/100倍降采样信号失3真就很严重了，已经听不清零了。需要注意的是调用sound函数时，格式为sound(y,fs,nbits)，对于1/M倍采样信号，fs应改为fs/M，这样才能等效还原语音。还有进行1/10倍,1/100倍降采样前，应该对原始信号进行低通滤波，以满足fs/2M大于信号最高频率4kHz。图3-2：1/2倍、1/10倍、1/100倍降采样信号波形和频谱图4.语音信号重构降采样后，信号的采样率和采样点数同时变化。如要恢复原始信号，信号长度和采样频率须要变为原来同样大小。因此，必须对降采样信号插值重构，即通过升采样恢复信号长度和采样频谱。4.1升采样对信号音质影响对原始语音信号抗混叠滤波后进行1/N降采样后再进行N倍升采样，波形和频谱如图，代码见附录2。调用sound函数感受插值后的声音，发现N越大，恢复后的声音尖锐噪声越4明显。本文对N=2,5,10进行了测试，随着N增大，高频声音越来越大。图4-1：N=2,5时升采样信号波形和频谱图4.2反镜像滤波作用升采样恢复的信号存在高频噪音，采用反镜像滤波器滤波后，高频噪音得到有效抑制。反镜像滤波器也是低通滤波器，本文采用的是用firls函数生成的fir滤波器，具体代码见附录2。对比原始语音和N=2和N=5的恢复信号，发现N=2时能有效恢复原始信号，而N=5时恢复效果不理想。通过分析发现，原始语音信号的最高频率成分是4kHz，而N=5降采样时的截止频率fs/(2*N)=2.2kHz，通过反混叠滤波器已经滤除了原始信号一些有用的频率成分，故恢复的效果不好。5图4-2：反镜像滤波后信号波形和频谱5.总结通过本实验的练习，加深了镜像频谱的理解，心得如下：正常以fs采样的信号，比如采了N个点，则在频谱的第N点处表示频率fs.所以增加采样点数只是提高了频率的分辨率，而没有拓宽频谱的观察范围。采用升采样后，比如以L倍插值，相当于拓宽频谱的观察范围到L*fs所以能看到原频谱以fs为周期拓展后的结果，也就是镜像频谱。FFT运算只是根据点数来运算，和采样频率并没有关系。原始的采样信号点数不管多少，是信号真实的反应。而采用升采样插值后，采样点所对应的信号其实变成了新的信号。如果插的是零点，新的信号的频谱相当于原始信号频谱观察范围拓宽L倍。新频谱的最大点数N对应新的采样频率L*fs。升采样的反镜像滤波在时域上相当于平滑滤波为原来的1/L,说的是信号的幅值。反镜像滤波后，内插的零点都有了值，相当于原来的模拟信号用L倍的采样频率进行采样，故单位时间的采样点数增加了，但信号的幅值降为原来的1/L，故要乘以系数L才和原来相同。滤波后的频谱也满足实际的采样情况。本实验的不足之处有以下几点：（1）滤波器的设计无详细考虑，滤波器的类型和阶数都影响最后结果；（2）升采样和降采样的算法没有深入研究6附录1：语音降采样和频谱分析Matlab程序clearall[y,fs,nbits]=wavread('H0.wav');N=2;y=y-mean(y);%滤除直流分量y1=y(1:N:length(y));figure,subplot(2,2,1);plot(y);title('原始信号波形');xlabel('采样点数');gridon;subplot(2,2,2);num=length(y);plot((-num/2:num/2-1)/num*fs,abs(fftshift(fft(y))));title('原始信号频谱');xlabel('frequenceHz');gridon;subplot(2,2,3);plot(y1);title('1/2采样信号波形');xlabel('采样点数');gridon;subplot(2,2,4);num=length(y1);plot((-num/2:num/2-1)/num*fs,abs(fftshift(fft(y1))));title('新信号频谱');xlabel('frequenceHz');gridon;附录2：原始信号降采样后用升采样用重构程序clear;closeall;[y,fs,nbits]=wavread('H0.wav');y=y-mean(y);%滤除直流分量%N=5的抗混叠滤波器也是反镜像滤波器%fs=22.05kHz，fs/（2*5）=2.2kHz,因此相对截止频率为0.1L=100;%滤波器阶数，即点数B5=firls(L,[00.090.111],[1100]);%N=10的抗混叠滤波器也是反镜像滤波器%fs=22.05kHz，fs/（2*10）=1.1kHz,因此相对截止频率为0.05L=100;%滤波器阶数，即点数B10=firls(L,[00.0490.0511],[1100]);%N=2的反镜像滤波器7%fs=22.05kHz，fs/（2*2）=5.5kHz,因此相对截止频率为0.25L=100;%滤波器阶数，即点数B2=firls(L,[00.240.261],[1100]);%滤波后结果y5=filter(B5,1,y);y10=filter(B10,1,y);%降采样输出yy2=downsample(y,2);%1/2倍降采样yy5=downsample(y5,5);%1/5倍降采样yy10=downsample(y10,10);%升采样恢复x2=zeros(1,2*length(yy2));x2(1:2:length(y))=yy2;x5=zeros(1,5*length(yy5));x5(1:5:length(y))=yy5;x10=zeros(1,10*length(yy10));x10(1:10:length(y))=yy10;%反镜像滤波输出xx2=filter(B2,1,x2);xx5=filter(B5,1,x5);xx10=filter(B10,1,x10);%作图figure,num=length(y);subplot(3,2,1);plot(y);title('原始信号波形');xlabel('采样点数');gridon;subplot(3,2,2);plot((-num/2:num/2-1)/num*fs,abs(fftshift(fft(y))));title('原始信号频谱');xlabel('frequenceHz');gridon;subplot(3,2,3);plot(x2);title('N=2信号波形');xlabel('采样点数');gridon;subplot(3,2,4);plot((-num/2:num/2-1)/num*fs,abs(fftshift(fft(x2))));title('N=2信号频谱');xlabel('frequenceHz');gridon;subplot(3,2,5);8plot(x5);title('N=5信号波形');xlabel('采样点数');gridon;subplot(3,2,6);plot((-num/2:num/2-1)/num*fs,abs(fftshift(fft(x5))));title('N=5信号频谱');xlabel('frequenceHz');gridon;