采用bic准则的音频流分割程序

happy

采用bic准则的音频流分割程序，能够准确找出不同音频种类的分界点位置，如音乐和语音。

给程序有两个函数组成，详细如下：

audioSeg.m

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. function [y,fs,position] = audioSeg(filenm,clipLen)
   
3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
4. %
   
5. %This is main audio stream segment function,it can classfication the audio file to different segments.
   
6. %
   
7. %input:
   
8. %     filenm = the flie which need to be classification;
   
9. %     clipLen = segment clip length,unit is frame,default value is 100;
   
10. %output:
    
11. %     y = the points vector of the file;
    
12. %     fs = sample rate;
    
13. %     position = type change points;
    
14. %
    
15. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
16. %% Written by Qian Yong Gao                   %%
    
17. %% Data:June 23,2006                          %%
    
18. %% Update:June 23,2006 by Qian Yong Gao       %%
    
19. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
20. 
    
21. if nargin < 2
    
22.     clipLen = 100;
    
23. end
    
24.          
    
25. [fid,msg] = fopen(filenm,'rb');%open file
    
26. error(msg);
    
27. fseek(fid,8,-1);%file point rewind
    
28. header = fread(fid,4,'uchar');
    
29. 
    
30. % makesure '*wav' file format
    
31. if header' ~= 'WAVE'
    
32.     error('This file is not a wavfile');
    
33. else
    
34.     fseek(fid,8,0);
    
35.     iPCM = fread(fid,1,'ushort');%PCM code format
    
36.     iMono = fread(fid,1,'ushort');%Mono channel
    
37.     if (iPCM ~= 1) & (iMono ~= 1)
    
38.         error('This wavefile is not Mono PCM type');
    
39.     end
    
40. end
    
41. fs = fread(fid,1,'ushort');%sample rate
    
42. fseek(fid,8,0);
    
43. nBit = fread(fid,1,'ushort');%each sample bits
    
44. nByte = nBit/8;%each sample bytes;
    
45. header = fread(fid,4,'uchar');
    
46. if header' ~= 'data'
    
47.     error('The file is corrupt');
    
48. end   
    
49. if nByte == 1
    
50.     nData = fread(fid,1,'ulong'); %file samples number
    
51.     rid = 'uchar';
    
52. else
    
53.     nData = fread(fid,1,'ulong') / 2;
    
54.     rid = 'short';
    
55. end   
    
56. 
    
57. position = 0;%beginning position is '0'
    
58. nRead = 0; %n samples read
    
59. nFrame = 0; %n frame read
    
60. segLen = 0; %n clip read
    
61. y = [];%signal read
    
62. frameLen = floor(0.02 * fs); %frame length
    
63. frameMov = ceil(0.01 * fs); %frame step
    
64. overlap = frameLen - frameMov;%frame overlap
    
65. 
    
66. while nRead <= nData - frameLen
    
67.     y1 = fread(fid,frameLen,rid);%read a frame
    
68.     nRead = nRead + frameMov;%read samples number increase frameMov
    
69.     if fid == 'uchar'
    
70.         y1 = (y1 - 128) / 128 ;
    
71.         fseek(fid,-overlap,0);%8 bit,rewind file pointer overlap bites
    
72.     else
    
73.         y1 = y1 / 32768;
    
74.         fseek(fid,-2 * overlap,0);%16 bit,rewind file pointer 2*overlap bites
    
75.     end   
    
76.     nFrame = nFrame + 1;%increase read frame number
    
77.     y = [y;y1'];%read signal     
    
78.     %when read frame is clipLen's integer times,nClip increase 1
    
79.     segLen = segLen + 1;
    
80.     if segLen == 3*clipLen
    
81.         clipMat = y([nFrame-3*clipLen+1:nFrame],:);%clip matrix which wait for classifcation
    
82.         f = fft(clipMat'.*(hamming(frameLen)*ones(1,3*clipLen)),1024);
    
83.         f = abs(f);
    
84.         mfc = melfc(f,30,30);
    
85.         segN = BeyasInfoCrit(mfc,frameLen,frameMov,0.6);
    
86.         if (segN>=clipLen) & (segN<=2*clipLen)
    
87.             position = [position;nFrame-3*clipLen+segN];
    
88.             segLen = 3 * clipLen - segN;
    
89.         else
    
90.             segLen = 2 * clipLen;
    
91.         end
    
92.     end
    
93. end   
    
94. fclose(fid);
    
95. position(2:end) = ((position(2:end)-1)*frameMov+frameLen) / fs;
    
96. y = y';
    
97. y1 = y(:,end);
    
98. y(:,end) = [];
    
99. y(frameMov+1:end,:) = [];
    
100. y = y(:);
     
101. y = [y;y1];
     
102. t = [1:length(y)] / fs;
     
103. plot(t,y,'k');
     
104. hold on
     
105. for l = 1:length(position)
     
106.     line([position(l) position(l)], [min(y) max(y)], 'Color', 'r');
     
107. end   
     
108. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

复制代码

happy

BeyasInfoCrit.m

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. function segN = BeyasInfoCrit(c,framelen,framemov,labeta)
   
3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
4. %
   
5. % This function use BIC to find the change point of c
   
6. %input:
   
7. %     c = frame feature matrix;
   
8. %     framelen = frame length
   
9. %     framemov = frame move
   
10. %output:
    
11. %     segN = segment change point
    
12. %
    
13. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
14. %% Written by Qian Yong Gao                   %%
    
15. %% Data:October 25,2005                       %%
    
16. %% Update:May 30,2006 by Qian Yong Gao        %%
    
17. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
18. if nargin < 4
    
19.     labeta = 0.8;
    
20. end   
    
21. [N,fd] = size(c);
    
22. u0 = mean(c);
    
23. sigma0 = sum((c-u0(ones(1,N),:)).^2)/N;
    
24. det0 = max(prod(max(sigma0,eps)),realmin);%行列式的值
    
25. detBIC = zeros(1,N-1);
    
26. for m = 1:N-1  %计算假定断点在m(1~N-1)处的BIC值
    
27.     %计算前一部分的参数
    
28.     c1 = c(1:m,:);
    
29.     u1 = mean(c1);
    
30.     sigma1 = sum((c1-u1(ones(1,m),:)).^2)/m;;
    
31.     %以对角线乘积代替行列式的值，避免数据不够带来的估计不正确
    
32.     det1 = max(prod(max(sigma1,eps)),realmin);
    
33.    
    
34.     c2 = c(m+1:N,:);
    
35.     u2 = mean(c2);
    
36.     sigma2 = sum((c2-u2(ones(1,N-m),:)).^2)/(N-m);
    
37.     det2 = max(prod(max(sigma2,eps)),realmin);
    
38.     %计算两种假设带来的差值
    
39.     detBIC(m) = 0.5*(N*log(det0)-m*log(det1)-(N-m)*log(det2))...
    
40.         - labeta*0.25*fd*(fd+3)*log(N);
    
41. end
    
42. [DB,index] = max(detBIC);
    
43. if DB > 0
    
44.     %segN = framemov*(index-1) + framelen;
    
45.     segN = index;
    
46. else
    
47.     segN = N;%framemov*(N-1) + framelen;
    
48. end

复制代码

890516k

请问mfc = melfc(f,30,30);中的melfc是什么函数，如有源程序能否发给我890516k@163.com

jingxinshulin

我也想要一份melfc 函数的代码，能否发给我：jingxinshulin@163.com
谢谢！

MISKERY

我也想要一份melfc 函数的代码，能否发给我：z.heng@foxmail.com,谢谢了