A 3D alternative to spectrogram

There have occasionally been threads here on visual presentation of 
displaying time sequential spectra. Spectrograms (gray scale or color) 
are very popular. [e.g. 
http://en.wikipedia.org/wiki/Image:Praat-spectrogram-tatata.png ]

I sat for many hours in front of RF spectrum analyzers, so I'm not 
oriented to that format. Scilab is free and has reasonable graphics, so 
I wrote my own.

I post it (as it is relatively short) as some might find it useful and 
with hope of constructive comments. Enjoy.



// A 3D alternative to spectrogram popular in speech recognition and 
related fields
// (q.v.
//      http://en.wikipedia.org/wiki/Image:Praat-spectrogram-tatata.png
//      http://en.wikipedia.org/wiki/Spectrogram )
// Richard Owlett (rowlett@atlascomm.net) Aug 6,2008

meg = 1000000;          // I like to look at ridiculously large data sets
stacksize(10 * meg);
set("old_style","off");  // force new graphics mode as default

datafile = xgetfile('*.wav','C:\ATEST\',title = 'Chose an input WAV file')
[data,y] = loadwave(datafile);

maxpoints = y(7)/2;
rate = y(3);

time_plot = scf();
plot2d((1:maxpoints)/rate,data)


//  menu  for choosing fft parameters

labels=["First Window";"Last Window";"Window Offset";"Window Width"];
[ok,t0,tlast,toffset,twidth]=getvalue("Define FFT parameters (all times 
in seconds)", ...
   labels, list("vec",1,"vec",1,"vec",1,"vec",1),["0.";"1.";".01";".05";]);

// t0 -- start time of first fft window
// tlast -- start time of last fft window
// toffset -- offset between each window
// twidth -- fft window width

n0 = max(int(t0*rate),1);          // start point
offset = int(toffset*rate);  // offset in points
width = int(twidth*rate);    // window width 1024 points
nw = int((tlast-t0)/twidth); // max number of fft's
nc = int(width/2);           // other half of fft is redundant

zz = 1;
deltapoints = width-1;
while zz <= nw;
     stpoint = (zz-1) * deltapoints + n0;
     endpoint = stpoint + deltapoints;
     if (maxpoints >= endpoint) then
         v = data(stpoint:endpoint) .* window('kr',width,8.6) ;
         aa = abs(fft(v,1));
         f(zz,1:nc) = aa(1:nc);
         xx(zz,1:nc)=zz;
     end
     zz = zz+1;
end;

maxf = rate * .0005;        // plot spectrum in kHz
nr=size(f,1);
bb = [1:nc];
yy = bb(ones([1:nr]),:);

spectra_plot =scf();

drawlater           // postpone drawing to avoid annoying redraws

param3d1(xx', yy' * maxf/nc , 10 * f'/max(f))
a=get("current_axes");     //get the handle of the newly created axes
     a.axes_reverse = ["on","off","off"];
     a.data_bounds = [t0,0,0;nw,int(maxf),10];
//  a.x_ticks.locations   // reminder to improve visuals
     a.x_ticks.labels = string(t0+toffset*(a.x_ticks.locations-1));

h=a.children;               //get the handle of the param3d entity: an 
Compound composed of nr curves
     h.children(:).mark_mode="on";
     h.children(:).line_mode="off";
     h.children(:).mark_style = 0;
     h.children(:).mark_size=0;

xset("fpf"," ")             // null labels for each contour
contour(xx(:,1) ,yy(1,:)* maxf/nc ,10 * f/max(f),4,flag=[0 0 4])

drawnow