www.pudn.com > LPCToolbox.rar > SPECGRM.M

function array = specgrm(wave,segsize,nlap,ntrans,pr); 
% array = specgrm(wave,segsize,nlap,ntrans[,pr]) 
%  
% Defaults specgrm(wave,128,8,4,0.95)  
% For broadband, try specgrm(wave,64,8,4) 
% For narrowband,try specgrm(wave,512,16,2) 
% 
%   wave   : waveform to analyze 
% 
%   segsize: size of analysis frame 
% 
%   nlap   : number of hamming windows overlapping a point, i.e. the 
%            analysis frame moves in increments of segsize/nlap. So, 
%            total # of frames = (length(wave)/(segsize/nlap))-nlap+1 
% 
%   ntrans : factor by which transform is bigger than segment, 
%            i.e. for the FFT, the analysis frame is padded with  
%            segsize*(ntrans-1) zeros. Note:  
%                frequency resolution = Fs/(ntrans*segsize),  
%             => length(FFT(...)) = ntrans*segsize 
%             => # of positive frequencies = segsize*ntrans/2 
% 
%   pr     : amount of preemphasis (0 <= pr <= 1.0, def: 0.95) 
% 
%   array  : spectrogram, compressed with 4th root of power, filter- 
%            smoothed, and formatted for 'imagesc' display (array(1,:)  
%            is the HIGHEST frequency). Each column is one time 
%            slice, each row is one frequency. 
% 
% N.B: The signal is preemphasized 95%. Each segment is  
% hamming-windowed prior to the FFT.  
%  
% This function is based on the SPECTROGRAM.M function 
% in Malcolm Slaney's Auditory Toolbox: 
% http://rvl4.ecn.purdue.edu/~malcolm/interval/1998-010/ 
% 
 
if nargin < 5, pr=.95; end 
if nargin < 4; ntrans=4; end 
if nargin < 3; nlap=8; end 
if nargin < 2; segsize=128; end 
 
wave = filter([1 -pr],[1],wave(:)); % preemphasis 
 
s = length(wave); 
updatelen = segsize/nlap; 
fftlen = ntrans*segsize; 
nsegs = floor(s/updatelen)-nlap+1; 
array = zeros(ceil(fftlen/2),nsegs); 
positivefreq = [1:ceil(fftlen/2)]; 
window = hamming(segsize); 
 
for k = 1:nsegs, 
 segstart = ((k-1) * updatelen) + 1; 
 ii = [segstart : segstart+segsize-1]; 
 seg = abs( fft(window.*wave(ii),fftlen) ); 
 array(:,k) = seg(positivefreq); 
end 
 
% get the power spectrum & invert for display 
array = flipud(array.*array);   
 
array = filter([.2 1 .2], [1], array, [], 1); % smooth columns 
array = filter([.2 1 .2], [1], array, [], 2); % smooth rows  
 
% compress with square root of amplitude (fourth root of power) 
off = 0.0001*max(max(array)); 
array = (off+array).^0.25-off^0.25; 	 
array = 255/max(max(array))*array;