Daniel@0: function [array,raw] = spectrogram(wave,segsize,nlap,ntrans);
Daniel@0: %function array = spectrogram(wave,segsize,nlap,ntrans);
Daniel@0: % defaults spectrogram(wave,128,8,4)
Daniel@0: % nlap is number of hamming windows overlapping a point;
Daniel@0: % ntrans is factor by which transform is bigger than segment;
Daniel@0: % returns a spectrogram 'array' with fourth root of power,
Daniel@0: % filter smoothed and formatted for display.
Daniel@0: 
Daniel@0: % Added option to return raw spectrogram.... Malcolm 5/26/95
Daniel@0: % Added code so that input could be any direction ... Malcolm 5/26/95
Daniel@0: % (c) 1998 Interval Research Corporation  
Daniel@0: 
Daniel@0: if nargin < 4; ntrans=4; end
Daniel@0: if nargin < 3; nlap=8; end
Daniel@0: if nargin < 2; segsize=128; end
Daniel@0: 
Daniel@0: [r c] = size(wave);
Daniel@0: if (r < c)
Daniel@0: 	wave = filter([1 -0.95],[1],wave');
Daniel@0: else
Daniel@0: 	wave = filter([1 -0.95],[1],wave);
Daniel@0: end
Daniel@0: 
Daniel@0: s = length(wave);
Daniel@0: nsegs = floor(s/(segsize/nlap))-nlap+1;
Daniel@0: array = zeros(ntrans/2*segsize,nsegs);
Daniel@0: window = 0.54-0.46*cos(2*pi/(segsize+1)*(1:segsize)');
Daniel@0: for i = 1:nsegs
Daniel@0:  seg = zeros(ntrans*segsize,1); % leave half full of zeroes
Daniel@0:  seg(1:segsize) = ...
Daniel@0: 	 window.*wave(((i-1)*segsize/nlap+1):((i+nlap-1)*segsize/nlap));
Daniel@0:  seg = abs(fft(seg));
Daniel@0: 						% reverse for image display
Daniel@0:  array(:,i) = seg(((ntrans/2*segsize)+1):(ntrans*segsize));
Daniel@0: end
Daniel@0: 
Daniel@0: if nargout > 1
Daniel@0: 	raw = array;
Daniel@0: end
Daniel@0: 
Daniel@0: array = array .* array;			% back into power domain for smoothing
Daniel@0: 
Daniel@0: for i=1:nsegs    % smooth the spectral slices
Daniel@0:  array(:,i) = filter([.2 1 .2],[1],array(:,i));
Daniel@0: end
Daniel@0: 
Daniel@0: for i=1:ntrans/2*segsize    % smooth the channels
Daniel@0:  array(i,:) = filter([.2 1 .2],[1],array(i,:));
Daniel@0: end
Daniel@0: 
Daniel@0: % compress with square root of amplitude (fourth root of power)
Daniel@0: off = 0.0001*max(max(array));      	% low end stabilization offset,
Daniel@0: array = (off+array).^0.25-off^0.25; 	% better than a threshold hack!
Daniel@0: array = 255/max(max(array))*array;
Daniel@0: 
Daniel@0: 
Daniel@0: