wolffd@0
|
1 function [array,raw] = spectrogram(wave,segsize,nlap,ntrans);
|
wolffd@0
|
2 %function array = spectrogram(wave,segsize,nlap,ntrans);
|
wolffd@0
|
3 % defaults spectrogram(wave,128,8,4)
|
wolffd@0
|
4 % nlap is number of hamming windows overlapping a point;
|
wolffd@0
|
5 % ntrans is factor by which transform is bigger than segment;
|
wolffd@0
|
6 % returns a spectrogram 'array' with fourth root of power,
|
wolffd@0
|
7 % filter smoothed and formatted for display.
|
wolffd@0
|
8
|
wolffd@0
|
9 % Added option to return raw spectrogram.... Malcolm 5/26/95
|
wolffd@0
|
10 % Added code so that input could be any direction ... Malcolm 5/26/95
|
wolffd@0
|
11 % (c) 1998 Interval Research Corporation
|
wolffd@0
|
12
|
wolffd@0
|
13 if nargin < 4; ntrans=4; end
|
wolffd@0
|
14 if nargin < 3; nlap=8; end
|
wolffd@0
|
15 if nargin < 2; segsize=128; end
|
wolffd@0
|
16
|
wolffd@0
|
17 [r c] = size(wave);
|
wolffd@0
|
18 if (r < c)
|
wolffd@0
|
19 wave = filter([1 -0.95],[1],wave');
|
wolffd@0
|
20 else
|
wolffd@0
|
21 wave = filter([1 -0.95],[1],wave);
|
wolffd@0
|
22 end
|
wolffd@0
|
23
|
wolffd@0
|
24 s = length(wave);
|
wolffd@0
|
25 nsegs = floor(s/(segsize/nlap))-nlap+1;
|
wolffd@0
|
26 array = zeros(ntrans/2*segsize,nsegs);
|
wolffd@0
|
27 window = 0.54-0.46*cos(2*pi/(segsize+1)*(1:segsize)');
|
wolffd@0
|
28 for i = 1:nsegs
|
wolffd@0
|
29 seg = zeros(ntrans*segsize,1); % leave half full of zeroes
|
wolffd@0
|
30 seg(1:segsize) = ...
|
wolffd@0
|
31 window.*wave(((i-1)*segsize/nlap+1):((i+nlap-1)*segsize/nlap));
|
wolffd@0
|
32 seg = abs(fft(seg));
|
wolffd@0
|
33 % reverse for image display
|
wolffd@0
|
34 array(:,i) = seg(((ntrans/2*segsize)+1):(ntrans*segsize));
|
wolffd@0
|
35 end
|
wolffd@0
|
36
|
wolffd@0
|
37 if nargout > 1
|
wolffd@0
|
38 raw = array;
|
wolffd@0
|
39 end
|
wolffd@0
|
40
|
wolffd@0
|
41 array = array .* array; % back into power domain for smoothing
|
wolffd@0
|
42
|
wolffd@0
|
43 for i=1:nsegs % smooth the spectral slices
|
wolffd@0
|
44 array(:,i) = filter([.2 1 .2],[1],array(:,i));
|
wolffd@0
|
45 end
|
wolffd@0
|
46
|
wolffd@0
|
47 for i=1:ntrans/2*segsize % smooth the channels
|
wolffd@0
|
48 array(i,:) = filter([.2 1 .2],[1],array(i,:));
|
wolffd@0
|
49 end
|
wolffd@0
|
50
|
wolffd@0
|
51 % compress with square root of amplitude (fourth root of power)
|
wolffd@0
|
52 off = 0.0001*max(max(array)); % low end stabilization offset,
|
wolffd@0
|
53 array = (off+array).^0.25-off^0.25; % better than a threshold hack!
|
wolffd@0
|
54 array = 255/max(max(array))*array;
|
wolffd@0
|
55
|
wolffd@0
|
56
|
wolffd@0
|
57
|