% ex9.3
fn=input('Enter WAV filename : ','s');
[x,fs,nb]=wavread(fn);
%
NF=128;         % number of filters
BW=300;         % filter bandwidth
OR=1000;        % output frame rate
%
% get output energy sample points
x=reshape(x,[length(x) 1]);         % force to be a column
eidx=floor(fs/OR):floor(fs/OR):length(x);         % sampling indices (OR frame/sec)
ne=length(eidx);                    % number of output energy frames
%
% calculate frequency cut offs
cuts=linspace(100,fs/2-100-BW,NF);         % from 100Hz to fs/2-100Hz
%
% initialise the output energy table
e=zeros(ne,NF);
%
% build a smoothing filter
[sb,sa]=butter(4,2*(OR/2)/fs);          % low-pass at half output frame rate
%
% for each channel in turn 
for i=1:NF
    fprintf('Calculating channel %d (%.1f-%.1fHz)\n',i,cuts(i),cuts(i)+BW);
    % build the band-pass filter
    [b,a]=butter(4,[2*cuts(i)/fs 2*(cuts(i)+BW)/fs]);
    % filter the signal
    y=filter(b,a,x);
    % rectify and smooth
    sy=filter(sb,sa,abs(y));
    % downsample 
    e(:,i)=sy(eidx);
end
%
% calculate amplitude 50dB down from maximum
emin=max(max(abs(e)))/300;
%
% plot top 50dB as image
t=(0:ne-1)/OR;                  % vector of frame times
imagesc(t,cuts,20*log10(max(abs(e'),emin)/emin));
%
% label plot
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
%
% build and use a grey scale
lgrays=zeros(100,3);
for i=1:100
    lgrays(i,:) = 1-i/100;
end
colormap(lgrays);