the code for mfcc feature extraction is giiven below, and the code showing error at the time of running and the error also given below. how to rectify the error??

clear all;
close all;
% Step 0: Reading the File & initializing the Time and Freq.

    [x,fs]=audioread('cryrumble.wav');
    ts=1/fs;
    N=length(x);
    Tmax=(N-1)*ts;
    fsu=fs/(N-1);
    t=(0:ts:Tmax);
    f=(-fs/2:fsu:fs/2);
    figure, subplot(411),plot(t,x),xlabel('Time'),title('Original audio');
    subplot(412),plot(f,fftshift(abs(fft(x)))),xlabel('Freq (Hz)'),title('Frequency Spectrum');
% Step 1: Pre-Emphasis
% 
% a=[1];
%     b=[1 -0.95];
%     y=filter(b,a,x);
%     subplot(413),plot(t,y),xlabel('Time'),title('Signal After High Pass Filter - Time Domain');
%     subplot(414),plot(f,fftshift(abs(fft(y)))),xlabel('Freq (Hz)'),title('Signal After High Pass Filter - Frequency Spectrum');
% Step 2: Frame Blocking

     frameSize=1000;
%     frameOverlap=128;
%     frames=enframe(y,frameSize,frameOverlap);
%     NumFrames=size(frames,1);
frame_duration=0.06;
frame_len = frame_duration*fs;
framestep=0.01;
framestep_len=framestep*fs;
% N = length (x);
num_frames =floor(N/frame_len);
% new_sig =zeros(N,1);
% count=0;
% frame1 =x(1:frame_len);
% frame2 =x(frame_len+1:frame_len*2);
% frame3 =x(frame_len*2+1:frame_len*3);
frames=[];
for j=1:num_frames
     frame=x((j-1)*framestep_len + 1: ((j-1)*framestep_len)+frame_len);
%     frame=x((j-1)*frame_len +1 :frame_len*j);
%     identify the silence by finding frames with max amplitude less than
%     0.025
max_val=max(frame);
   if (max_val>0.025)
%     count = count+1;
%     new_sig((count-1)*frame_len+1:frame_len*count)=frames;
        frames=[frames;frame];
   end
end
% Step 3: Hamming Windowing

NumFrames=size(frames,1);
hamm=hamming(1000)';
    for i=1:NumFrames
    windowed(i,:)=frames(i,:).*hamm;
    end
    % Step 4: FFT 
% Taking only the positive values in the FFT that is the first half of the frame after being computed. 


    for i=1:NumFrames
    ft(i,:)=abs(fft((windowed(i,:)),480));     
    plot(ft(i,:))
    end
    
% Step 5: Mel Filterbanks
Lower_Frequency = 100;
Upper_Frequency = fs/2;
% With a total of 22 points we can create 20 filters.
    Nofilters=20;
    lowhigh=[300 fs/2];
    %Here logarithm is of base 'e'
    lh_mel=1125*(log(1+lowhigh/700));
    mel=linspace(lh_mel(1),lh_mel(2),Nofilters+2);
    melinhz=700*(exp(mel/1125)-1);
    %Converting to frequency resolution
    fres=floor(((frameSize)+1)*melinhz/fs); 
    %Creating the filters
    for m =2:length(mel)-1
        for k=1:frameSize/2
     if k<fres(m-1)
        H(m-1,k) = 0;
    elseif (k>=fres(m-1)&&k<=fres(m))
        H(m-1,k)= (k-fres(m-1))/(fres(m)-fres(m-1));
    elseif (k>=fres(m)&&k<=fres(m+1))
       H(m-1,k)= (fres(m+1)-k)/(fres(m+1)-fres(m));
    elseif k>fres(m+1)
        H(m-1,k) = 0;    
     end 
        end
    end
        %H contains the 20 filterbanks, we now apply it to the
    %processed signal.
    for i=1:NumFrames
    for j=1:Nofilters
        bankans(i,j)=sum((ft(i,:).*H(j,:)).^2);
    end
    end
    % Step 6: Nautral Log and DCT

%     pkg load signal
    %Here logarithm is of base '10'
    logged=log10(bankans);
    for i=1:NumFrames

        mfcc(i,:)=dct2(logged(i,:));
    end
    %plotting the MFCC
    figure 
    hold on
    for i=1:NumFrames
        plot(mfcc(i,1:13));
    end
    hold off
% save c5 mfcc

i= mfcc;
save i i
load i.mat
X=i;
k=1;
[IDXi,ci] = kmeans(X,k);
save c41i ci
mfcccombined.m
Displaying mfcccombined.m.





when I running this code  it showing error like thi" Reference to a cleared variable frames.
Error in mfccaudio (line 57)
    windowed(i,:)=frames(i,:).*hamm;" 
    how to rectify the error??

[Walter Roberson]

function hamming probably contains

clear all

[Suchithra K S]

sir, i am not getting can you please clear that.

Answer 1

[Walter Roberson]

I copied your code in and applied it to a random .wav that I had around.

Your code is not efficient. You have

hamm=hamming(1000)';
for i=1:NumFrames
windowed(i,:)=frames(i,:).*hamm;
end

You can prove that at that point frames is a single column vector, so you are effectively multiplying each element of that column vector by the hamm matrix, in a loop. You can be much much more efficient by coding

windowed = bsxfun(@times, frames, hamm);

with no loop.

You also have

    for i=1:NumFrames
    ft(i,:)=abs(fft((windowed(i,:)),480));     
    plot(ft(i,:))
    end

NumFrames is easily in the tens of thousands. "hold" is not on in the axis, so every iteration you are plotting something that will be ovewritten by the next iteration, which is pointless work: just plot the final iteration after the loop. You can also be much more efficient on doing the fft in a loop: just fft everything at once with no loop:

    ft = abs( fft(windowed, 480, 2) );
    

The ,2 is because you are fft'ing row by row and the default for fft() is column by column.

Eventually you get to

    for i=1:NumFrames
        for j=1:Nofilters
            bankans(i,j)=sum((ft(i,:).*H(j,:)).^2);
        end
    end

Because of the ,480 for the fft, we know that ft has 480 columns. The 480 columns are being multiplied by a column of H. How many columns does H have?

   for k=1:frameSize/2
     if k<fres(m-1)
        H(m-1,k) = 0;
    elseif (k>=fres(m-1)&&k<=fres(m))
        H(m-1,k)= (k-fres(m-1))/(fres(m)-fres(m-1));
    elseif (k>=fres(m)&&k<=fres(m+1))
       H(m-1,k)= (fres(m+1)-k)/(fres(m+1)-fres(m));
    elseif k>fres(m+1)
        H(m-1,k) = 0;    
     end
   end

So we see that it has frameSize/2 columns, and frameSize=1000 so half of that is 500. So you are trying to multiply 480 columns by 500 columns which is going to fail.

If the number of columns did match then you could no doubt improve efficiency by vectorizing, removing the inner loop.

bankans(i,:) = sum(bsxfun(@times, ft(i,:), H).^2,2);

Answer 2

[Suchithra K S]

Sir, run the above code with these chnges then also the code showing error. The error is like this "checking

Warning: Ignoring rows of X with missing data.

> In kmeans (line 164)

In checking (line 127)

Undefined variable "mffccombined" or class "mffccombined.m".

Error in checking (line 129)

mffccombined.m"

Sir, can you please rectify this error.

[Walter Roberson]

I think the last two liness are intended to convey that the author of the code had stored the code in a file named mfcccombined.

[Walter Roberson]

what is your current code version ?

[Suchithra K S]

matlab 2016 version.

[Walter Roberson]

what is your code after you modified it? II cannot assume that you fixed it the same way that II would have fixed it.

[Suchithra K S]

my modified code is given below,

[audio, fs1] = audioread('cryrumble.wav');
%sound(x,fs1);
   ts1=1/fs1;
    N1=length(audio);
    Tmax1=(N1-1)*ts1;
    t1=(0:ts1:Tmax1);
    figure;
  plot(t1,audio),xlabel('Time'),title('Original audio');
  fs2 = (20/441)*fs1;
  y=resample(audio,2000,44100);
%sound(y,fs2);
   ts2=1/fs2;
   N2=length(y);
   Tmax2=(N2-1)*ts2;
   t2=(0:ts2:Tmax2);
   figure;
   plot(t2,y),xlabel('Time'),title('resampled audio');
%   % spectrogram(y);
  nchan = size(y,2);
for chan = 1 : nchan
  %subplot(1, nchan, chan)
  spectrogram(y(:,chan), 256, [], 25, 2000, 'yaxis');
  title( sprintf('spectrogram of resampled audio ' ) );
end

% Step 2: Frame Blocking

     frameSize=1000;
%     frameOverlap=128;
%     frames=enframe(y,frameSize,frameOverlap);
%     NumFrames=size(frames,1);
frame_duration=0.06;
frame_len = frame_duration*fs2;
framestep=0.01;
framestep_len=framestep*fs2;
% N = length (x);
num_frames =floor(N2/frame_len);
% new_sig =zeros(N,1);
% count=0;
% frame1 =x(1:frame_len);
% frame2 =x(frame_len+1:frame_len*2);
% frame3 =x(frame_len*2+1:frame_len*3);
frames=[];
for j=1:num_frames
     frame=y((j-1)*framestep_len + 1: ((j-1)*framestep_len)+frame_len);
%     frame=x((j-1)*frame_len +1 :frame_len*j);
%     identify the silence by finding frames with max amplitude less than
%     0.025
max_val=max(frame);
   if (max_val>0.025)
%     count = count+1;
%     new_sig((count-1)*frame_len+1:frame_len*count)=frames;
    frames=[frames;frame];
   end
   end
   % Step 3: Hamming Windowing

NumFrames=size(frames,1);
hamm=hamming(1000)';
    for i=1:NumFrames
    windowed(i,:)=frames(i,:).*hamm;
    end
    
      
     % Step 4: FFT 
% Taking only the positive values in the FFT that is the first half of the frame after being computed. 


    for i=1:NumFrames
    ft(i,:)=abs(fft((windowed(i,:)),480));     
    plot(ft(i,:))
    end
    
% Step 5: Mel Filterbanks
Lower_Frequency = 100;
Upper_Frequency = fs2/2;
% With a total of 22 points we can create 20 filters.
    Nofilters=20;
    lowhigh=[300 fs2/2];
    %Here logarithm is of base 'e'
    lh_mel=1125*(log(1+lowhigh/700));
    mel=linspace(lh_mel(1),lh_mel(2),Nofilters+2);
    melinhz=700*(exp(mel/1125)-1);
    %Converting to frequency resolution
    fres=floor(((frameSize)+1)*melinhz/fs2); 
    %Creating the filters
    for m =2:length(mel)-1
        for k=1:frameSize/2
     if k<fres(m-1)
        H(m-1,k) = 0;
    elseif (k>=fres(m-1)&&k<=fres(m))
        H(m-1,k)= (k-fres(m-1))/(fres(m)-fres(m-1));
    elseif (k>=fres(m)&&k<=fres(m+1))
       H(m-1,k)= (fres(m+1)-k)/(fres(m+1)-fres(m));
    elseif k>fres(m+1)
        H(m-1,k) = 0;    
     end 
        end
    end
        %H contains the 20 filterbanks, we now apply it to the
    %processed signal.
    for i=1:NumFrames
    for j=1:Nofilters
        bankans(i,j)=sum((ft(i,:).*H(j,:)).^2);
    end
    end
    % Step 6: Nautral Log and DCT

%     pkg load signal
    %Here logarithm is of base '10'
    logged=log10(bankans);
    for i=1:NumFrames

        mfcc(i,:)=dct2(logged(i,:));
    end
    %plotting the MFCC
    figure 
    hold on
    for i=1:NumFrames
        plot(mfcc(i,1:13));
    end
    hold off
% save c5 mfcc

i= mfcc;
save i i
load i.mat
X=i;
k=1;
[IDXi,ci] = kmeans(X,k);
save c41i ci
mfcccombined.m
Displaying mfcccombined.m.

[Walter Roberson]

where did you make even one of the changes I discussed??

[Suchithra K S]

sorry Sir, while I uploading the code is changed. This is my code.

[audio, fs1] = audioread('cryrumble.wav');
%sound(x,fs1);
   ts1=1/fs1;
    N1=length(audio);
    Tmax1=(N1-1)*ts1;
    t1=(0:ts1:Tmax1);
    figure;
  plot(t1,audio),xlabel('Time'),title('Original audio');
  fs2 = (20/441)*fs1;
  y=resample(audio,2000,44100);
%sound(y,fs2);
   ts2=1/fs2;
   N2=length(y);
   Tmax2=(N2-1)*ts2;
   t2=(0:ts2:Tmax2);
   figure;
   plot(t2,y),xlabel('Time'),title('resampled audio');
%   % spectrogram(y);
  nchan = size(y,2);
for chan = 1 : nchan
  %subplot(1, nchan, chan)
  spectrogram(y(:,chan), 256, [], 25, 2000, 'yaxis');
  title( sprintf('spectrogram of resampled audio ' ) );
end

% Step 2: Frame Blocking

     frameSize=1000;
%     frameOverlap=128;
%     frames=enframe(y,frameSize,frameOverlap);
%     NumFrames=size(frames,1);
frame_duration=0.06;
frame_len = frame_duration*fs2;
framestep=0.01;
framestep_len=framestep*fs2;
% N = length (x);
num_frames =floor(N2/frame_len);
% new_sig =zeros(N,1);
% count=0;
% frame1 =x(1:frame_len);
% frame2 =x(frame_len+1:frame_len*2);
% frame3 =x(frame_len*2+1:frame_len*3);
frames=[];
for j=1:num_frames
     frame=y((j-1)*framestep_len + 1: ((j-1)*framestep_len)+frame_len);
%     frame=x((j-1)*frame_len +1 :frame_len*j);
%     identify the silence by finding frames with max amplitude less than
%     0.025
max_val=max(frame);
   if (max_val>0.025)
%     count = count+1;
%     new_sig((count-1)*frame_len+1:frame_len*count)=frames;
    frames=[frames;frame];
   end
   end
   % Step 3: Hamming Windowing

NumFrames=size(frames,1);
hamm=hamming(1000)';
windowed = bsxfun(@times, frames, hamm);

      
     % Step 4: FFT 
% Taking only the positive values in the FFT that is the first half of the frame after being computed. 


       ft = abs( fft(windowed,500, 2) );
       plot(ft);
    
% Step 5: Mel Filterbanks
Lower_Frequency = 100;
Upper_Frequency = fs2/2;
% With a total of 22 points we can create 20 filters.
    Nofilters=20;
    lowhigh=[300 fs2/2];
    %Here logarithm is of base 'e'
    lh_mel=1125*(log(1+lowhigh/700));
    mel=linspace(lh_mel(1),lh_mel(2),Nofilters+2);
    melinhz=700*(exp(mel/1125)-1);
    %Converting to frequency resolution
    fres=floor(((frameSize)+1)*melinhz/fs2); 
    %Creating the filters
    for m =2:length(mel)-1
        for k=1:frameSize/2
     if k<fres(m-1)
        H(m-1,k) = 0;
    elseif (k>=fres(m-1)&&k<=fres(m))
        H(m-1,k)= (k-fres(m-1))/(fres(m)-fres(m-1));
    elseif (k>=fres(m)&&k<=fres(m+1))
       H(m-1,k)= (fres(m+1)-k)/(fres(m+1)-fres(m));
    elseif k>fres(m+1)
        H(m-1,k) = 0;    
     end 
        end
    end
        %H contains the 20 filterbanks, we now apply it to the
    %processed signal.
    for i=1:NumFrames
    for j=1:Nofilters
        bankans(i,j)=sum((ft(i,:).*H(j,:)).^2);
    end
    end
    % Step 6: Nautral Log and DCT

%     pkg load signal
    %Here logarithm is of base '10'
    logged=log10(bankans);
    for i=1:NumFrames

        mfcc(i,:)=dct2(logged(i,:));
    end
    %plotting the MFCC
    figure 
    hold on
    for i=1:NumFrames
        plot(mfcc(i,1:13));
    end
    hold off
% save c5 mfcc

i= mfcc;
save i i
load i.mat
X=i;
k=1;
[IDXi,ci] = kmeans(X,k);
save c41i ci
mffccombined.m
Displaying mfcccombined.m

[Walter Roberson]

delete the last two lines.

[Suchithra K S]

ok thank you sir.

[Suchithra K S]

Sir in the above code I used mel frequency filterbank instead of that I have to use greenwood filter bank. Can you tell me the code.

[Walter Roberson]

Sorry, I do not know anything about greenwood filter banks.

[Suchithra K S]

The given below attachments shows the green wood filter and i dont know how to write the matlab code for this .Sir, can you please give me the code for that.

[Walter Roberson]

No, I cannot give you the code for that. Perhaps it would be sufficient for you to use https://www.mathworks.com/matlabcentral/fileexchange/32212-gammatone-filterbank

[Suchithra K S]

ok thank you sir.

[Adarsh Sharma]

Hi, can you please highlight where did you do the pre-emphasis, the first step of getting MFCC features?

[hanif omar]

did you get for the first step ? pre emphasis ?