function mfcc=frame2mfcc(frame, fs, filterNum, mfccNum, plotOpt)
% frame2mfcc: Frame to MFCC conversion.
%    Usage: mfcc=frame2mfcc(frame, fs, filterNum, mfccNum, plotOpt)
%
%    For example:
%        waveFile='what_movies_have_you_seen_recently.wav';
%        [y, fs, nbits]=wavReadInt(waveFile);
%        startIndex=12000;
%        frameSize=512;
%        frame=y(startIndex:startIndex+frameSize-1);
%        frame2mfcc(frame, fs, 20, 12, 1);%    Roger Jang 20060417if nargin<1, selfdemo; return; end
if nargin<2, fs=16000; end
if nargin<3, filterNum=20; end
if nargin<4, mfccNum=12; end
if nargin<5, plotOpt=0; endframeSize=length(frame);
% ====== Preemphasis should be done at wave level
%a=0.95;
%frame2 = filter([1, -a], 1, frame);
frame2=frame;
% ====== Hamming windowing
frame3=frame2.*hamming(frameSize);
% ====== FFT
[fftMag, fftPhase, fftFreq, fftPowerDb]=fftOneSide(frame3, fs);
% ====== Triangular band-pass filter bank
triFilterBankPrm=getTriFilterBankPrm(fs, filterNum);    % Get parameters for triangular band-pass filter bank
% Triangular bandpass filter.
for i=1:filterNumtbfCoef(i)=dot(fftPowerDb, trimf(fftFreq, triFilterBankPrm(:,i)));%得到filterNum个滤波系数
end
% ====== DCT
mfcc=zeros(mfccNum, 1); %DCT变换的前后个数也没有变
for i=1:mfccNumcoef = cos((pi/filterNum)*i*((1:filterNum)-0.5))'; %mfcc中的前mfccNum个系数mfcc(i) = sum(coef.*tbfCoef');%直接按照DCT公式
end
% ====== Log energy
%logEnergy=10*log10(sum(frame.*frame));
%mfcc=[logEnergy; mfcc];if plotOptsubplot(2,1,1);plot(frame, '.-');set(gca, 'xlim', [-inf inf]);title('Input frame');subplot(2,1,2);plot(mfcc, '.-');set(gca, 'xlim', [-inf inf]);title('MFCC vector');
end% ====== trimf.m (from fuzzy toolbox)
function y = trimf(x, prm) %由频率的横坐标算出三角形内的纵坐标,0~1
a = prm(1); b = prm(2); c = prm(3);
y = zeros(size(x));
% Left and right shoulders (y = 0)
index = find(x <= a | c <= x);
y(index) = zeros(size(index)); %只考虑三角波内的量
% Left slope
if (a ~= b)index = find(a < x & x < b);y(index) = (x(index)-a)/(b-a);
end
% right slope
if (b ~= c)index = find(b < x & x < c);y(index) = (c-x(index))/(c-b);
end
% Center (y = 1)
index = find(x == b);
y(index) = ones(size(index));% ====== Self demo
function selfdemo
waveFile='what_movies_have_you_seen_recently.wav';
[y, fs, nbits]=wavReadInt(waveFile);
startIndex=12000;
frameSize=512;
frame=y(startIndex:startIndex+frameSize-1);
feval(mfilename, frame, fs, 20, 12, 1);

waveFile='csNthu.wav';
frameSize=256;
overlap=0;
[y, fs, nbits]=wavread(waveFile);
frameMat=enframe(y, frameSize, overlap);
frameNum=size(frameMat, 2);
for i=1:frameNumframeMat(:,i)=frameMat(:,i)-mean(frameMat(:,i));    % mean justification
end
zcr=sum(frameMat(1:end-1, :).*frameMat(2:end, :)<0);
sampleTime=(1:length(y))/fs;
frameTime=((0:frameNum-1)*(frameSize-overlap)+0.5*frameSize)/fs;
subplot(2,1,1); plot(sampleTime, y); ylabel('Amplitude'); title(waveFile);
subplot(2,1,2); plot(frameTime, zcr, '.-');
xlabel('Time (sec)'); ylabel('Count'); title('ZCR');

waveFile='sunday.wav';
[wave, fs, nbits] = wavread(waveFile);
frameSize = 256;
overlap = 128;wave=wave-mean(wave);                % zero-mean substraction
frameMat=buffer2(wave, frameSize, overlap);    % frame blocking,每一列代表一帧
frameNum=size(frameMat, 2);            % no. of frames
volume=frame2volume(frameMat);        % volume,求每一帧的能量，绝对值或者平方和,volume为行向量
volumeTh1=max(volume)*0.1;            % volume threshold 1
volumeTh2=median(volume)*0.1;            % volume threshold 2
volumeTh3=min(volume)*10;            % volume threshold 3
volumeTh4=volume(1)*5;                % volume threshold 4
index1 = find(volume>volumeTh1); %找出volume大于阈值的那些帧序号
index2 = find(volume>volumeTh2);
index3 = find(volume>volumeTh3);
index4 = find(volume>volumeTh4);
%frame2sampleIndex()为从帧序号找到样本点的序号(即每一个采样点的序号)
%endPointX长度为2,包含了起点和终点的样本点序号
endPoint1=frame2sampleIndex([index1(1), index1(end)], frameSize, overlap);
endPoint2=frame2sampleIndex([index2(1), index2(end)], frameSize, overlap);
endPoint3=frame2sampleIndex([index3(1), index3(end)], frameSize, overlap);
endPoint4=frame2sampleIndex([index4(1), index4(end)], frameSize, overlap);subplot(2,1,1);
time=(1:length(wave))/fs;
plot(time, wave);
ylabel('Amplitude'); title('Waveform');
axis([-inf inf -1 1]);
line(time(endPoint1(  1))*[1 1], [-1, 1], 'color', 'm');%标起点终点线
line(time(endPoint2(  1))*[1 1], [-1, 1], 'color', 'g');
line(time(endPoint3(  1))*[1 1], [-1, 1], 'color', 'k');
line(time(endPoint4(  1))*[1 1], [-1, 1], 'color', 'r');
line(time(endPoint1(end))*[1 1], [-1, 1], 'color', 'm');
line(time(endPoint2(end))*[1 1], [-1, 1], 'color', 'g');
line(time(endPoint3(end))*[1 1], [-1, 1], 'color', 'k');
line(time(endPoint4(end))*[1 1], [-1, 1], 'color', 'r');
legend('Waveform', 'Boundaries by threshold 1', 'Boundaries by threshold 2', 'Boundaries by threshold 3', 'Boundaries by threshold 4');subplot(2,1,2);
frameTime=frame2sampleIndex(1:frameNum, frameSize, overlap);
plot(frameTime, volume, '.-');
ylabel('Sum of Abs.'); title('Volume');
axis tight;
line([min(frameTime), max(frameTime)], volumeTh1*[1 1], 'color', 'm');
line([min(frameTime), max(frameTime)], volumeTh2*[1 1], 'color', 'g');
line([min(frameTime), max(frameTime)], volumeTh3*[1 1], 'color', 'k');
line([min(frameTime), max(frameTime)], volumeTh4*[1 1], 'color', 'r');
legend('Volume', 'Threshold 1', 'Threshold 2', 'Threshold 3', 'Threshold 4');

function [M, V, W, logProb] = gmmTrain(data, gaussianNum, dispOpt)
% gmmTrain: Parameter training for gaussian mixture model (GMM)
%    Usage: function [M, V, W, logProb] = gmm(data, gaussianNum, dispOpt)
%        data: dim x dataNum matrix where each column is a data point
%        gaussianNum: No. of Gaussians or initial centers
%        dispOpt: Option for displaying info during training
%        M: dim x meanNum matrix where each column is a mean vector
%        V: 1 x gaussianNum vector where each element is a variance for a Gaussian
%        W: 1 x gaussianNum vector where each element is a weighting factor for a Gaussian% Roger Jang 20000610if nargin==0, selfdemo; return; end
if nargin<3, dispOpt=0; endmaxLoopCount = 50;    % Max. iteration
minImprove = 1e-6;    % Min. improvement
minVariance = 1e-6;    % Min. variance
logProb = zeros(maxLoopCount, 1);   % Array for objective function
[dim, dataNum] = size(data);% Set initial parameters
% Set initial M
%M = data(1+floor(rand(gaussianNum,1)*dataNum),:);    % Randomly select several data points as the centers
if length(gaussianNum)==1,% Using vqKmeans to find initial centersfprintf('Start KMEANS to find the initial mu...\n');
%    M = vqKmeansMex(data, gaussianNum, 0);M = vqKmeans(data, gaussianNum, 0); %利用聚类的方法求均值,聚成gaussianNum类
%    M = vqLBG(data, gaussianNum, 0);fprintf('Start GMM training...\n');if any(any(~isfinite(M))); keyboard; end
else% gaussianNum is in fact the initial centersM = gaussianNum;gaussianNum = size(M, 2);
end
% Set initial V as the distance to the nearest center
if gaussianNum==1V=1;
elsedistance=pairwiseSqrDist(M);%pairwiseSqrDist是dll%distance=pairwiseSqrDist2(M);distance(1:(gaussianNum+1):gaussianNum^2)=inf;    % Diagonal elements are inf[V, index]=min(distance);    % Initial variance for each Gaussian
end
% Set initial W
W = ones(1, gaussianNum)/gaussianNum;    % Weight for each Gaussian,初始化时是均分权值if dispOpt & dim==2, displayGmm(M, V, data); end
for i = 1:maxLoopCount  %开始迭代训练参数,EM算法% Expectation step:% P(i,j) is the probability of data(:,j) to the i-th Gaussian% Prob为每个样本在GMM下的概率[prob, P]=gmmEval(data, M, V, W);logProb(i)=sum(log(prob)); %所有样本的联合概率if dispOptfprintf('i = %d, log prob. = %f\n',i-1, logProb(i));endPW = diag(W)*P;BETA=PW./(ones(gaussianNum,1)*sum(PW));    % BETA(i,j) is beta_i(x_j)sumBETA=sum(BETA,2);% Maximization step:  eqns (2.96) to (2.98) from Bishop p.67:M = (data*BETA')./(ones(dim,1)*sumBETA');DISTSQ = pairwiseSqrDist(M, data);                    % Distance of M to data%DISTSQ = pairwiseSqrDist2(M, data);                    % Distance of M to dataV = max((sum(BETA.*DISTSQ, 2)./sumBETA)/dim, minVariance);    % (2.97)W = (1/dataNum)*sumBETA;                    % (2.98)if dispOpt & dim==2, displayGmm(M, V, data); endif i>1, if logProb(i)-logProb(i-1)<minImprove, break; end; end
end
[prob, P]=gmmEval(data, M, V, W);
logProb(i)=sum(log(prob));
fprintf('Iteration count = %d, log prob. = %f\n',i, logProb(i));
logProb(i+1:maxLoopCount) = [];% ====== Self Demo ======
function selfdemo
%[data, gaussianNum] = dcdata(2);
data = rand(1000,2);
gaussianNum = 8;
data=data';
plotOpt=1;
[M, V, W, lp] = feval(mfilename, data, gaussianNum, plotOpt);pointNum = 40;
x = linspace(min(data(1,:)), max(data(1,:)), pointNum);
y = linspace(min(data(2,:)), max(data(2,:)), pointNum);
[xx, yy] = meshgrid(x, y);
data = [xx(:) yy(:)]';
z = gmmEval(data, M, V, W);
zz = reshape(z, pointNum, pointNum);
figure; mesh(xx, yy, zz); axis tight; box on; rotate3d on
figure; contour(xx, yy, zz, 30); axis image% ====== Other subfunctions ======
function displayGmm(M, V, data)
% Display function for EM algorithm
figureH=findobj(0, 'tag', mfilename);
if isempty(figureH)figureH=figure;set(figureH, 'tag', mfilename);colordef blackplot(data(1,:), data(2,:),'.r'); axis imagetheta=linspace(-pi, pi, 21);x=cos(theta); y=sin(theta);sigma=sqrt(V);for i=1:length(sigma)circleH(i)=line(x*sigma(i)+M(1,i), y*sigma(i)+M(2,i), 'color', 'y');endset(circleH, 'tag', 'circleH', 'erasemode', 'xor');
elsecircleH=findobj(figureH, 'tag', 'circleH');theta=linspace(-pi, pi, 21);x=cos(theta); y=sin(theta);sigma=sqrt(V);for i=1:length(sigma)set(circleH(i), 'xdata', x*sigma(i)+M(1,i), 'ydata', y*sigma(i)+M(2,i));enddrawnow
end

function [recogRate, confusionMatrix, speakerData]=speakerIdentify(speakerData, speakerGmm, useIntGmm)
% speakerIdentify: speaker identification using GMM parameters
%    Usage: [recogRate, confusionMatrix, speakerData]=speakerIdentify(speakerData, speakerGmm, useIntGmm)
%        speakerData: structure array generated by speakerDataRead.m
%        speakerGmm: speakerGmm(i).gmmPrm is the GMM parameters for speaker i.
%        useIntGmm: use fixed-point GMM%    Roger Jang, 20070517, 20080726if nargin<3, useIntGmm=0; end% ====== Speaker identification using GMM parameters
speakerNum=length(speakerData);
for i=1:speakerNum
%    fprintf('%d/%d: Recognizing wave files by %s\n', i, speakerNum, speakerData(i).name);for j=1:length(speakerData(i).sentence)
%        fprintf('\tSentece %d...\n', j);frameNum=size(speakerData(i).sentence(j).fea, 2);logProb=zeros(speakerNum, frameNum); %logProb(i,m)表示第i个人第j个句子中第m帧在GMM模型下的log概率%找出一个句子，看它属于哪个speakerfor k=1:speakerNum,
%            fprintf('\t\tSpeaker %d...\n', k);%    logProb(k, :)=gmmEval(speakerData(i).sentence(j).fea, speakerGmm(k).gmmPrm);if ~useIntGmm%    logProb(k, :)=gmmEvalMex(speakerData(i).sentence(j).fea, gmm(k).mean, gmm(k).covariance, gmm(k).weight);logProb(k, :)=gmmEval(speakerData(i).sentence(j).fea, speakerGmm(k).gmmPrm);else%    logProb(k, :)=gmmEvalIntMex(speakerData(i).sentence(j).fea, gmm(k).mean, gmm(k).covariance, gmm(k).weight);logProb(k, :)=gmmEvalIntMex(speakerData(i).sentence(j).fea, speakerGmm(i).gmmPrm);endendcumLogProb=sum(logProb, 2);[maxProb, index]=max(cumLogProb);speakerData(i).sentence(j).predictedSpeaker=index; %找出身份speakerData(i).sentence(j).logProb=logProb;end
end% ====== Compute confusion matrix and recognition rate
confusionMatrix=zeros(speakerNum);
for i=1:speakerNum,predictedSpeaker=[speakerData(i).sentence.predictedSpeaker];[index, count]=elementCount(predictedSpeaker);confusionMatrix(i, index)=count;
end
recogRate=sum(diag(confusionMatrix))/sum(sum(confusionMatrix));