Matlab音频信号读取函数audioread

昨天有人问我怎么用Matlab读取WAV或MP3等音频信号，然后提取音频，并作傅里叶变换，找到音频的振频和谐波。然后，他给了我一段m文件代码，说是会报错，让我改一下看看。然后，我就一顿操作，使用audioread函数读取，稍微修改了两行代码，就帮他解决了问题。

今天，把这个过程及最后的代码发出来，供各位参考。虽然，我也不知道到底啥意思，反正就是解决那些报错信息就可以了。

原版m文件

接下来，先把他给的m文件代码发出来：

clc; 
clear all;
[x,FS,NBITS]=wavread('1.WAV'); 
N=length(x)
if mod(N,2)==0;N=N;else x(N)=[];N=N-1;end;
tx=(0:N-1)/FS;
subplot(3,2,1);plot(tx,x);
xf=fft(x);
fx=(0:N/2)*FS/N;
subplot(3,2,2);plot(fx,abs(xf(1:N/2+1)));

运行后，报错信息如下：

未定义函数或变量 'wavread'。

出错 dai (line 3)
[x,FS,NBITS]=wavread('1.WAV');

问题判断：不存在wavread这个函数，估计是用法问题。经过一番查询，觉得可以使用audioread这个函数，我们先看看这个函数的功能吧：

>> help audioread
 audioread Read audio files
    [Y, FS]=audioread(FILENAME) reads an audio file specified by the 
    string FILE, returning the sampled data in Y and the sample rate 
    FS, in Hertz. 
    
    [Y, FS]=audioread(FILENAME, [START END]) returns only samples START 
    through END from each channel in the file.
    
    [Y, FS]=audioread(FILENAME, DATATYPE) specifies the data type format of
    Y used to represent samples read from the file.
    If DATATYPE='double', Y contains double-precision normalized samples.
    If DATATYPE='native', Y contains samples in the native data type
    found in the file.  Interpretation of DATATYPE is case-insensitive and
    partial matching is supported.
    If omitted, DATATYPE='double'.  
    
    [Y, FS]=audioread(FILENAME, [START END], DATATYPE);
 
    Output Data Ranges
    Y is returned as an m-by-n matrix, where m is the number of audio 
    samples read and n is the number of audio channels in the file.
 
    If you do not specify DATATYPE, or dataType is 'double', 
    then Y is of type double, and matrix elements are normalized values 
    between -1.0 and 1.0.
 
    If DATATYPE is 'native', then Y may be one of several MATLAB 
    data types, depending on the file format and the BitsPerSample 
    of the input file:
    
     File Format      BitsPerSample  Data Type of Y     Data Range of Y
     ----------------------------------------------------------------------
     WAVE (.wav)            8           uint8             0 <= Y <= 255
                           16           int16        -32768 <= Y <= 32767
                           24           int32         -2^32 <= Y <= 2^32-1
                           32           int32         -2^32 <= Y <= 2^32-1
                           32           single         -1.0 <= Y <= +1.0
     ----------------------------------------------------------------------
     FLAC (.flac)           8           uint8             0 <= Y <= 255
                           16           int16        -32768 <= Y <= 32767
                           24           int32         -2^32 <= Y <= 2^32-1
     ----------------------------------------------------------------------
     MP3 (.mp3)            N/A          single         -1.0 <= Y <= +1.0
     MPEG-4(.m4a,.mp4)
     OGG (.ogg)
     ----------------------------------------------------------------------
 
    Call audioinfo to learn the BitsPerSample of the file.
 
    Note that where Y is single or double and the BitsPerSample is 
    32 or 64, values in Y might exceed +1.0 or -1.0.

我们发现，读取后的音频会有两个参数：一个Y，是一个矩阵，一般是一个两列多行矩阵，甚至多行多列矩阵的采样数据；FS是这个音频的采样率，即一秒钟的音频所包含的数据组数。

因此，我们把对方发来的m文件，作一部分修改，得到如下结果：

clc; 
clear all;
[Y,FS]=audioread('1.WAV'); 
N=length(Y)
if mod(N,2)==0;N=N;else Y(N)=[];N=N-1;end;
tx=(0:N-1)/FS;
subplot(3,2,1);
plot(tx,Y);
xf=fft(Y);
fx=(0:N/2)*FS/N;
subplot(3,2,2);
plot(fx,abs(xf(1:N/2+1)));

运行后，得到如下图形：

对方，觉得振频那部分，即右边处理的不是太好。经过观察，我发现需要调整下振频部分的横坐标范围就可以了。顺便，我重新调整了图片放置的位置，通过修改subplot这个函数可以实现。最后，修改后的代码如下：

clc; 
clear all;
[Y,FS]=audioread('1.WAV'); 
N=length(Y)
if mod(N,2)==0;N=N;else Y(N)=[];N=N-1;end;
tx=(0:N-1)/FS;
subplot(2,1,1);
plot(tx,Y);
xf=fft(Y);
fx=(0:N/2)*FS/N;
subplot(2,1,2);
plot(fx,abs(xf(1:N/2+1)));
axis([100 300,-inf,inf])

其运行结果如下：

好了，到此完美的解决了对方的问题。这里说一下，axis，即最后一行，是修改plot画图结果的坐标范围：其中，100 300表示横坐标范围为100到300，-inf inf表示纵坐标范围不做处理。

当然，也可以读取mp3文件，比如我这里再处理一个mp3文件，如下结果：

开头，我们说Y是一个两列矩阵，其实我们也可以通过如下代码分别展开这两列，来绘图，结果跟文中的m文件结果是一致的。算了，最后这部分可以忽略。其实plot自己会识别多列数组的，它会自动在一个图中绘制结果。

a=Y(:,1);
b=Y(:,2);
plot(tx,a,tx,b)
plot(tx,Y)