搜索
MATLAB在音频信号处理中的应用1.引言MATLAB是美国MathWorks公司推出的一种面向工程和科学计算的交互式计算软件、它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程实用软件。MATLAB软件除可以对文本数据进行操作外，还可以对图像文件、WAV类型的音频文件和AVI类型的视频文件进行读写.借助windows系统自带的音频录制播放设备和MATLAB软件，可以完成音频信号的综合分析处理。音频是信号的一种，处理数字音频信号也是一种数字信号分析与处理。人能听到的声音频率范围为20Hz到20kHz，语音信号频率为300Hz到3.4kH。木文首先对音频信号的采集与频谱分析作介绍，然后介绍了用MATLAB处理音频信号的整体流程，并且对我们平常所听的音乐文件进行一系列的处理，然后用我们的耳朵亲自去听，去感受不同处理方法的效果。2.音频信号的采集与频谱分析进行音频的频谱分析时，首先要对声音信号进行采。MATLAB的数据采集工具箱提供了一整套命令和函数,通过调用这些函数和命令,可直接控制声卡进行数据采集。Windows自带的录音机程序也可驱动声卡来采集语音信号，并能保存为WAV格式文件供MATLAB相关函数直接读取、写入或播放。函数[y,fs,bits]=wavread(‘Blip’,N)。用于读取音频，对音频信号进行采集，采样值放在参数y中，fs表示每秒采样点数，即采样频率，bits表示每个采样点在编码时所占位数。N表示采样点总数。参数’Blip’为音频所在地址，如:'C:\yinpinl'。调用函数fft可对己采集音频信号进行时频转换，通过函数abs()和angle()可分别得到信号频谱的幅频图和相频图。对放在C盘目录下手机短信铃声”test.wav”进行采样与频谱分析。其体代码如下:clearall;[y0,Fs0,nbits0]=wavread('test.wav');a=round(length(y0));y01=fft(y0);y011=abs(y01);t0=linspace(0,a/Fs0,a);figure(1);subplot(2,1,1);plot(t0,y0);title('(a)','fonts',10.5,'position',[25,-1.5]);xlabel('时间/s','fonts',10.5,'position',[45,-1.4]);ylabel('幅值n','fonts',10.5,'position',[-3,0.7]);f0=linspace(0,Fs0,a);figure(2);subplot(2,1,1);plot(f0,y011);title('(a)','fonts',10.5,'position',[8000,-780]);xlabel('频率f/hz','fonts',10.5,'position',[14000,-600]);ylabel('幅值n','fonts',10.5,'position',[-1500,2000]);figure(3);subplot(2,1,1);plot(f0(1:round(a/16)),y011(1:round(a/16)));title('(a)','fonts',10.5,'position',[500,-780]);xlabel('频率f/hz','fonts',10.5,'position',[850,-600]);ylabel('幅值n','fonts',10.5,'position',[-90,2000]);运行程序，得到音频信y0，采样率Fs0=44100/s;采样总点数a=441000;音频总时长t0=10s。通过图1音频信号时域波形图，可以看出音频频率的大小，但不够定量，频率分布定量计算需要通过频谱分析。通过MATLAB对音频进行处理，得到图1的音频信号时域信号波形图，和图2、图3的幅频图。图1音频信号时域信号波形图图2音频信号幅频图图3音频信号频谱密集部分幅频图由图2的音频信号幅频图可以看图，此音频信号的频率主要集中在低频区和高频区，此外我们还可以用MATLAB计算频率在1000hz以下的频谱部分占0-45000hz频谱部分的比例。具体实现语句为：w0=(sum(y011(1:10000)))/sum(y011(1:45000));运行结果w0=0.5324。3.MATLAB处理音频信号的流程我们听到的声音一般都属于复音,其声音信号由不同的振幅与频率的波合成而得。WAV格式音频信号作为分析处理的输入数据。用MATLAB处理音频信号的基本流程是先将WAV格式音频信号经wavread函数转换成MATLAB列数组变量，再用MATLAB对音频信号进行音量标准化，声道分离合并与组合，数字滤波，数据转换。处理后的数据如是音频数据则可用wavwrite转换成WAV格式文件或用soundwavplay等函数直接回放。3.1音量标准化录制声音过程中需对声音电平进行量化处理，最理想的量化是最大电平对应最高量化比特。但实际却很难做到，常有音轻问题。利用MATLAB很容易实现音量标准化。即最大电平对应最高量化比特。基本步骤是：先用wavread函数将WAV文件转换成列数组变量，再求出数组变量的极值并对所有元素作归一化处理，最后用wavwrite函数还原成音量标准化的WAV文件。本文现对手机短信铃声的音频信号进行处理。先将其复制另存到MATLAB当前目录中，文件名为test.wav。再通过音量标准化处理后保存为outtest.wav文件。试听可知标准化处理后音量稍大,在MATLAB运行结果中采样频率FS=44100,量化比特NBITS=16，Ym=1.3.2声道分离合并与组合立体声或双声道音频信号有左右两个声道。利用MATLAB实现双声道分离，两路声道合并和两个单声道组合成一个双声道等效果。实际上是利用了MATLAB的矩阵抽取，矩阵相加和矩阵重组运算现以音量标准化处理生成的outtest.wav实现分离，合并和组合处理。可以试听声道分离，合并与组合的效果。也可对各文件大小进行比较。3.3数字滤波数字滤波是常用的音频处理技术，可根据技术指标先利用FDATool工具设计一个数字滤波器，再用Filter或Filter2函数即可实现滤波处理。调用的Filter函数格式是Y=filte(B,A,X)。其中B和A是滤波器传输函数的分子和分母系数,X是输入变量，Y是实现滤波后的输出变量。如果处理立体声音频信号，可分开处理，但用FIR滤波器时调用Filter2函数更方便。现以声道分离合并与组合生成的outtest12.wav实现数字滤波。程序运行结果如图3所示，由图可知滤波对波形影响不大，但对高频有较大衰减，试听会感觉到处理后的声音比较沉闷。图3数字滤波前后的波形对比3.4数据转换数据转换是指改变音频格式中的采样频率或量化位数。转换原理是先用矩阵插值或抽取技术实现变量变换，如果是抽取数据还需在变换前作滤波处理使之满足采样定理。变量变换完成后再用Wavwrite函数重新定义量化位数和采样频率即可实现数据转换。数据转换过程中要注意采样频率与原始采样频率及插值或抽取系数的关系。MATLAB实现插值或抽取的函数有decimate，interp和resample。这里以2倍抽取为例，将经过数字滤波后产生的outtest12Filter.WAV文件进行数据转换处理。运行数据转换程序，在得到的图形窗口中执行Edit/AxesProperties命令，再把各分图下X标签中的Limits设为0、0.01和0、1000，得到0--0.01秒的波形和0--1000Hz的频谱如图4所示。由图可知在满足采样定律条件下实现数据抽取。在原采样率下，波形变密、频谱变宽且幅度减半，但在新采样率下波形和频谱都很好。通过试听输出文件还可感受处理效果。图4数据转换波形对比图音量标准化程序：clear;closeall;clc;[Y,FS,NBITS]=wavread('test.WAV');%将WAV文件转换成变量FS,NBITS,%显示采样频率和量化比特Ym=max(max(max(Y)),max(abs(min(Y)))),%找出双声道极值X=Y/Ym;%归一化处理wavwrite(X,FS,NBITS,'outtest.wav')%将变量转换成WAV文件声道分离合并与组合程序：clear;closeall;clc;[x,FS,NBITS]=wavread('outtest.WAV');%将WAV文件转换成变量x1=x(:,1);%抽取第1声道x2=x(:,2);%抽取第2声道wavwrite(x1,FS,NBITS,'outtest1.WAV');%实现1声道分离wavwrite(x2,FS,NBITS,'outtest2.WAV');%实现2声道分离%如果合并位置不对前面补0%声道长度不对后面补0x12=x1+x2;%两路单声道列向量矩阵变量合并x12m=max(max(x12),abs(min(x12))),%找出极值y12=x12./x12m;%归一化处理wavwrite(y12,FS,NBITS,'outtest12.WAV');%实现两路声道合并%如果组合位置不对前面补0--声道长度不对后面补0x3=[x1,x2];%两路单声道变量组合wavwrite(x3,FS,NBITS,'outtest13.WAV');%实现两路声道组合数字滤波程序：clear;closeall;clc;[X,FS,NBITS]=wavread('outtest12.WAV');%将WAV文件转换成变量%利用FDATool设计一个LowpassButterworth滤波器%指标FS=22050HzFp=1000HzAp=1dBFs=3000HzAs=20dBB=[0.0062,0.0187,0.0187,0.0062];%分子系数A=[1,-2.1706,1.6517,-0.4312];%分母系数Y=filter(B,A,X);%实现数字滤波t=(0:length(X)-1)/FS;%计算数据时刻subplot(2,2,1);plot(t,X);%绘制原波形图title('原信号波形图');%加标题subplot(2,2,3);plot(t,Y);%绘制滤波波形图title('滤波后波形图');%加标题xf=fft(X);%作傅里叶变换求原频谱yf=fft(Y);%作傅里叶变换求滤波后频谱fm=3000*length(xf)/FS;%确定绘频谱图的上限频率f=(0:fm)*FS/length(xf);%确定绘频谱图的频率刻度subplot(2,2,2);plot(f,abs(xf(1:length(f))));%绘制原波形频谱图title('原信号频谱图');%加标题subplot(2,2,4);plot(f,abs(yf(1:length(f))));%绘制滤波后频谱图title('滤波后信号频谱图');%加标题wavwrite(Y,FS,NBITS,'outtest12Filter.WAV');%写成WAV文件数据转换程序：clear;closeall;clc;[x,FS,NBITS]=wavread('outtestFilter.WAV');%将WAV文件转换成变量N=length(x);%计算数据点数%不是偶数点化成偶数点ifmod(N,2)==0;N=N;elsex(N)=[];N=N-1;end;%原信号波形频谱分析tx=(0:N-1)/FS;%计算原信号数据点时刻subplot(3,2,1);plot(tx,x);%绘制原信号波形title('原信号波形图');%加标题xf=fft(x);%求原信号频谱fx=(0:N/2)*FS/N;%确定频谱图频率刻度subplot(3,2,2);plot(fx,abs(xf(1:N/2+1)));%绘制原信号频谱title('原信号频谱图');%加标题%实现数据抽取k=[1