% function [c1,a]=melcepst(s,fs,w,nc,p,n,inc,fl,fh)

[fname,pname] = uigetfile('*.wav')
s = readwav([pname fname],'r');
s = s(1:256)
n =256;
inc = 64;
fs = 8000;
w ='Mpt';
nc =16;
p = 24;
fl = 0;
fh = 0.5;
%MELCEPST Calculate the mel cepstrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH)
%
%
% Simple use: c=melcepst(s,fs)	% calculate mel cepstrum with 12 coefs, 256 sample frames
%				  c=melcepst(s,fs,'e0dD') % include log energy, 0th cepstral coef, delta and delta-delta coefs
%
% Inputs:
%     s	 speech signal
%     fs  sample rate in Hz (default 11025)
%     nc  number of cepstral coefficients excluding 0'th coefficient (default 12)
%     n   length of frame (default power of 2 <30 ms))
%     p   number of filters in filterbank (default floor(3*log(fs)) )
%     inc frame increment (default n/2)
%     fl  low end of the lowest filter as a fraction of fs (default = 0)
%     fh  high end of highest filter as a fraction of fs (default = 0.5)
%
%		w   any sensible combination of the following:
%
%				'R'  rectangular window in time domain
%				'N'	Hanning window in time domain
%				'M'	Hamming window in time domain (default)
%
%		      't'  triangular shaped filters in mel domain (default)
%		      'n'  hanning shaped filters in mel domain
%		      'm'  hamming shaped filters in mel domain
%
%				'p'	filters act in the power domain
%				'a'	filters act in the absolute magnitude domain (default)
%
%			   '0'  include 0'th order cepstral coefficient
%				'e'  include log energy
%				'd'	include delta coefficients (dc/dt)
%				'D'	include delta-delta coefficients (d^2c/dt^2)
%
%		      'z'  highest and lowest filters taper down to zero (default)
%		      'y'  lowest filter remains at 1 down to 0 frequency and
%			   	  highest filter remains at 1 up to nyquist freqency
%
%		       If 'ty' or 'ny' is specified, the total power in the fft is preserved.
%
% Outputs:	c     mel cepstrum output: one frame per row
%


%      Copyright (C) Mike Brookes 1997
%
%      Last modified Thu Jun 15 09:14:48 2000
%
%   VOICEBOX is a MATLAB toolbox for speech processing. Home page is at
%   http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if nargin<2 fs=11025; end
% if nargin<3 w='M'; end
% if nargin<4 nc=12; end
% if nargin<5 p=floor(3*log(fs)); end   %  if fs =8000 then p = 26
% if nargin<6 n=pow2(floor(log2(0.03*fs))); end   % if fs = 8000 then n = 128
% if nargin<9
%    fh=0.5;   
%    if nargin<8
%      fl=0;
%     if nargin<7
%        inc=floor(n/2);
%     end
%  end
% end

% 預加重濾波器
s=filter([1 -0.9375],1,s)
z=enframe(s,hamming(n),inc)
f=rfft(z.')
pw= f.*conj(f)
m=melbankm(p,n,fs,fl,fh)
%%m = m/max(m(:));
%figure(1)
%[a1,b1] = size(pw)
%[x1,y1] = meshgrid(1:1:a1, 1:1:b1);
%plot3(x1,y1,pw')
%title('pw' )


pth=max(pw(:))*1E-6

log_pw = log(max(pw,pth))
y=log(max(m*pw,pth))
%figure(3)
%[a1,b1] = size(pw)
%[x1,y1] = meshgrid(1:1:a1, 1:1:b1);
%plot3(x1,y1,log_pw')
%title('log_pw' )

%figure(6)
%[x1,y1] = meshgrid(1:1:b2, 1:1:a2);
%plot3(x1,y1,y)
%title('y ')

c=rdct(y).'

%plot(c(1:100))
nf=size(c,1)
nc=nc+1
if p>nc
   c(:,nc+1:end)=[]
elseif p<nc
   c=[c zeros(nf,nc-p)]
end
if ~any(w=='0')
   c(:,1)=[];
end
%c(1:100)

%figure(9)
%[a3,b3] = size(c)
%[x1,y1] = meshgrid(1:1:a3, 1:1:b3);
%plot3(x1,y1,c')
%title('c ' );

[a,b] = size(c);
%[fname1,pname1] = uiputfile('*.hdg')
 %  filename1 = [pname1, fname1] 
 
 c(:,nc) = [1:1:a]';
 %fprintf('size(c) = %f,%f',size(c));
 
 c1(1:2) =[a,b];
 %-----c1(3:a*(b+1)+2) = reshape(c',1,a*(b+1));
 
   %count  = fwrite(fid,16,'float32')
   % c(1:10,1)'
   % c(:,end+1) = [1:1:a]';
   % c(1:10,1)'
    
  %  c1 = reshape(c', a*(b+1),1);
   %  c1(1:10,1)'
  %  c(1:10,1)'
    
    
    
    
%figure(10)
%[x1,y1] = meshgrid(1:1:b3, 1:1:a3);
%plot3(x1,y1,c)
%title('c ')


%% calculate derivative

%if any(w=='D')
%  vf=(2:-1:-2);
%  af=(1:-1:-1);
%  ww=ones(2,1)
  
%  cx=[c(ww,:); c;c(nf*ww,:)];
%  vx=reshape(filter(vf,1,cx(:)),nf+4,nc-1);
%  vx = vx/3;
    
%  ax=reshape(filter(af,1,vx(:)),nf+4,nc-1)
%  ax([1 2 nf+1 nf+2],:)=[]
%  ax = ax/2;
  
%  vx([1 2 nf+1 nf+2],:) = [];
  
%  figure(11)
%  [x7,y7] = meshgrid(1:1:nf, 1:1:nc-1);
%  plot3(x7,y7,vx')
%  title('dc');
  
%  figure(12)
%  [x7,y7] = meshgrid(1:1:nc-1, 1:1:nf);
%  plot3(x7,y7,vx)
% title('dc');
  
%  figure(13)
%  [x7,y7] = meshgrid(1:1:nf, 1:1:nc-1);
%  plot3(x7,y7,ax')
%  title('ddc');
  
% figure(14)
 %[x7,y7] = meshgrid(1:1:nc-1, 1:1:nf);
 %plot3(x7,y7,ax)
% title('ddc');
  
 % if any(w=='d')
 %    c=[c vx ax];
 %else
 %    c=[c ax];
 % end
 %elseif any(w=='d')
 % vf=(2:-1:-2);
 % ww=ones(2,1);
 % cx=[c(ww,:); c; c(nf*ww,:)];
 % vx=reshape(filter(vf,1,cx(:)),nf+4,nc-1);
 % vx([1 2 nf+1 nf+2],:) = [];
 % vx = vx/3;
 % c=[c vx];
 %end
 
%%soundview(sunion,8000)

%% if nargout<1
%%   [nf,nc]=size(c)
%%   t=((0:nf-1)*inc+(n-1)/2)/fs
%%   ci=(1:nc)-any(w=='0')-any(w=='e')
%%   imh = imagesc(t,ci,c.')
%%   axis('xy')
%%   xlabel('Time (s)')
%%   ylabel('Mel-cepstrum coefficient')
%%	map = (0:63)'/63
%%	colormap([map map map])
%%	colorbar
    % end