matlab声音处理



由于最近的Project要做声音分析,需要用到MATLAB,之前一直没怎么接触过,所以乘着做Project学习下。真的用了才知道MATLAB真是神器啊,呵呵~~~其强大的函数库和数学运算能力彻底让我折服了。言归正传,我们来讨论下用MATLAB做声音文件处理。

1. 读取WAV声音文件

  1. % wavread(filename) 读取一个WAVE文件,并返回采样数据到向量y中,Fs表示采样频率, bits表示采样位数  
  2. [y, Fs, bits] = wavread('drum.wav');  
  3.   
  4. %假设声音文件有两个声道,我们只分析第一个声道,如果要分析第二个声道可以改为:ft=y(:,2);  
  5. ft = y(:,1);  
  6. sigLength = length(ft); %获取声音长度  
  7.   
  8. %可以使用sound函数来播放声音  
  9. sound(y, Fs, bits)  
% wavread(filename) 读取一个WAVE文件,并返回采样数据到向量y中,Fs表示采样频率, bits表示采样位数
[y, Fs, bits] = wavread('drum.wav');

%假设声音文件有两个声道,我们只分析第一个声道,如果要分析第二个声道可以改为:ft=y(:,2);
ft = y(:,1);
sigLength = length(ft); %获取声音长度

%可以使用sound函数来播放声音
sound(y, Fs, bits)


2. 绘制波形图
  1. t=(0:sigLength-1)/Fs;   
  2. figure;  
  3. subplot(2,1,1);  
  4. plot(t, ft), title('Plot of the Tone'),grid;  
  5. xlabel('Time(s)');  
  6. ylabel('Amplitude');  
t=(0:sigLength-1)/Fs; 
figure;
subplot(2,1,1);
plot(t, ft), title('Plot of the Tone'),grid;
xlabel('Time(s)');
ylabel('Amplitude');

3. 绘制振幅频谱图
  1. %Y = fft(X) 使用快速傅里叶变换算法返回向量X的离散型傅里叶变换  
  2. %Y = fft(X,n) 返回n点的离散傅里叶变换,如果向量X的长度小于n,函数要将向量X补零到长度n;如果向量X的长度大于n, 则函数阶段X使之长度为n。若X是矩阵,按相同方法对X进行处理。  
  3.   
  4. Y = fft(ft,sigLength);   
  5. halfLength = floor(sigLength/2);  
  6. Pyy =Y(1:halfLength + 1); % 只选取前半截部分  
%Y = fft(X) 使用快速傅里叶变换算法返回向量X的离散型傅里叶变换
%Y = fft(X,n) 返回n点的离散傅里叶变换,如果向量X的长度小于n,函数要将向量X补零到长度n;如果向量X的长度大于n, 则函数阶段X使之长度为n。若X是矩阵,按相同方法对X进行处理。

Y = fft(ft,sigLength); 
halfLength = floor(sigLength/2);
Pyy =Y(1:halfLength + 1); % 只选取前半截部分

波形的傅里叶变换返回震级和相位信息,并且用复数的形式表达,通过计算绝对值来获取其频率的振幅
  1. Pyy = abs(Pyy);%用于计算复向量的Y的振幅  
  2. f = ((0:halfLength)+1)* Fs/sigLength;   
  3. subplot(2,1,2);  
  4. plot(f,Pyy), title('Frequency spectrum'),grid;  
  5. xlabel('Frequency(Hz)');  
  6. ylabel('Amplitude');  
Pyy = abs(Pyy);%用于计算复向量的Y的振幅
f = ((0:halfLength)+1)* Fs/sigLength; 
subplot(2,1,2);
plot(f,Pyy), title('Frequency spectrum'),grid;
xlabel('Frequency(Hz)');
ylabel('Amplitude');

4. 绘制能量频谱图
  1. Y = fft(ft,sigLength);   
  2. halfLength = floor(sigLength/2);  
  3. Pyy =Y(1:halfLength + 1); % 只选取前半截部分  
  4. Pyy = abs(Pyy);%用于计算复向量的Y的振幅  
  5. f = ((0:halfLength)+1)* Fs/sigLength;   
  6. %通过点数调整比例,从而使振幅不依赖于信号长度或采样频率,说实话这部分还不是很明白  
  7. Pyy = Pyy/sigLength;   
  8. Pyy = Pyy.^2; %求平方得到能量  
  9. % 乘以2, 请参照 http://www.mathworks.com/support/tech-notes/1700/1702.html  
  10. if rem(sigLength, 2)  %奇数的 n fft 不包含奈奎斯特(Nyquist)点  
  11.    Pyy(2:end) = Pyy(2:end)*2;  
  12. else  
  13.    Pyy(2:end -1) = Pyy(2:end-1)*2;  
  14. end  
  15.   
  16. plot(f/1000, 10*log10(Pyy), 'k')   
  17. xlabel('Frequency (kHz)')   
  18. ylabel('Power (dB)')   
Y = fft(ft,sigLength); 
halfLength = floor(sigLength/2);
Pyy =Y(1:halfLength + 1); % 只选取前半截部分
Pyy = abs(Pyy);%用于计算复向量的Y的振幅
f = ((0:halfLength)+1)* Fs/sigLength; 
%通过点数调整比例,从而使振幅不依赖于信号长度或采样频率,说实话这部分还不是很明白
Pyy = Pyy/sigLength; 
Pyy = Pyy.^2; %求平方得到能量
% 乘以2, 请参照 http://www.mathworks.com/support/tech-notes/1700/1702.html
if rem(sigLength, 2)  %奇数的 n fft 不包含奈奎斯特(Nyquist)点
   Pyy(2:end) = Pyy(2:end)*2;
else
   Pyy(2:end -1) = Pyy(2:end-1)*2;
end

plot(f/1000, 10*log10(Pyy), 'k') 
xlabel('Frequency (kHz)') 
ylabel('Power (dB)') 



引用:
The MethWorks Support. Technique notes 1702,  http://www.mathworks.com/support/tech-notes/1700/1702.html
  • 14
    点赞
  • 62
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值