spectrogram_v6.m

Spectrogram function in matlab [For speech]

function  []=Spectrogram_V6(inp,Fs,segsize,overlap,PREEMPH,SCF,DR,CLR,XL,YL,labell)

%DEFAULT:
%segsize=128; overlap=segsize/2; PREEMPH=1; SCF=0; DR=-80; CLR=0; XL=1; YL=1;
%Spectrogram_V6(inp,Fs,segsize,overlap,PREEMPH,SCF,DR,CLR,XL,YL,'(a)');


%   SEGSIZE        => SEG SIZE
%   OVERLAP        => OVERLAP SIZE 
%   PREEMPH        => 1 FOR PREEMPHASIS
%   SCF            => Smoothing Coefficient (-0.2)
%   DR             => Dynamic Range (-80 dB); 
%   CLR            => 1 FOR COLOR SPECTROGRAM   
%   XL, YL         => 1 TO DISPLAY LABEL 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Spectrogram Final: Default Parameters

% (i)Seg size: Highly change the nature of the spectrogram
% Low (64)     ->   High Contrast
% High (256)   ->  Low Contrast
 
% (ii)	Overlap size: It will not affect much
% Default= segsize/2
% 
% (iii)	 Preemphasis:  Increase the high frequency component value
% Default =1
 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
A=0.99;   % PREEMPHASIS FILTER COEFFICIENT  
SPEC_EXP=0.1;

% (iv)	Spec_Exp  :   Higly the affect the nature of the spectrogram
% Default:  0.1  (<0.1)  Low  ->  High Contrast (0.1)  High  -> Low Contrast  (0.5)



disp('------------------- SPECTROGRAM DISPLAY -------------------');

if(size(inp,1)==1)
    inp=inp';
end


upFreq=1.0;   % Upper frequency (percentage of Srate/2) in spectrogram
Srate=Fs;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

yLab = [0.5*Srate*upFreq/1000 0];
ex=(length(inp))/Srate;
xLab = [0 ex];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    meen=mean(inp);
    inp=inp-meen;
    if norm(inp,2) < 1.0e-7
      if meen < 1.0e-3
        meen=1.0e-3;
      end
      inp=inp+meen;
    end
    
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    nSamples=length(inp);
    
    if nSamples<1000
      ns=4;
    elseif nSamples>1000 & nSamples<4000
      ns=3;
    else
      ns=3;
    end
    
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if CLR == 0
		nMap=length(colormap);
		x1=gray(nMap);
		x4=flipud(x1);
		colormap(x4);
	else
	   colormap('jet');
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   
 image(xLab,yLab,spectrg_V2(inp,segsize,overlap,ns,upFreq,PREEMPH,A,SPEC_EXP,SCF,DR));
 axis xy
 text(max(xLab)+0.05,range(yLab)/2,labell)

 
 if(XL==1)
 xlabel('Time (secs)');
 end
 
 if(YL==1)
 ylabel('Freq. (KHz)');
 end
 
 if(CLR==0)
 set(gca,'Xcolor','k'); set(gca,'Ycolor','k');
 else
 set(gca,'Xcolor','k'); set(gca,'Ycolor','k');
 end     
   
%**************************************************************************

function spec = spectrg_V2(wave,segsize,overlap,ntrans,upfr,PREEMPH,A,SPEC_EXP,SCF,DR)
 
%`````````````````````````````````````````````````````````````````````````````````````````````````````````````` 
%PREEMPHASIS

if PREEMPH>0
   wave = filter([1 -A],1,wave);
end
%`````````````````````````````````````````````````````````````````````````````````````````````````````````````` 
%COMPUTATION OF STFT

slap=overlap;

s = length(wave);
 
nsegs = floor(s/slap)-floor(segsize/slap)+1;
 
slap=floor(slap);
 
ysize=floor(ntrans/2*segsize*upfr);
spec = zeros(ysize,nsegs);
 
seg2=zeros(ysize,1);
 
window = hanning(segsize);

%`````````````````````````````````````````````````````````````````````````````````````````````````````````````` 
%---take the FFT of  the waveform ----

k=1;

for i = 1:nsegs
 
 seg = zeros(ntrans*segsize,1); % leave half full of zeroes
 seg(1:segsize) = window.*wave(k:k+segsize-1);
 
 seg = abs(fft(seg));
 seg2=seg(1:ysize);
 
 spec(:,i)=flipud(filter(1,[1 SCF],seg2.*seg2));
 
 k=k+slap;
 
end 

%`````````````````````````````````````````````````````````````````````````````````````````````````````````````` 
% smooth the channels

for i=1:ysize     
    spec(i,:) = filter(1,[1 SCF],spec(i,:));
end

%`````````````````````````````````````````````````````````````````````````````````````````````````````````````` 

% =================== compress with root of amplitude  ================

    offsetvalue=10^(DR/20);
    off = offsetvalue*max(max(spec));
    spec = (off+spec).^SPEC_EXP-off^SPEC_EXP;   
    spec = 255/max(max(spec))*spec;

%``````````````````````````````````````````````````````````````````````````````````````````````````````````````