% edrAndPlot.m  
% calculate and plot the Energy Decay Relief of a signal
%
% Music 421, Spring Quarter 2003-2004

function B_EDRdb = edrAndPlot(s,mode)

fileName = 'SLucia_22050mono.wav'; % name of soundfile to be analyzed
% addpath ~/mus421/hw3/; % location of soundfile
minPlotDB = -60; % minimum value in DB to be plotted
frameSizeMS = 30; % minimum frame length, in ms
overlap = 0.5; % fraction of frame overlapping
windowType = 'hann'; % type of windowing used for each frame

[signal,fs,bits] = wavread(fileName);
if (mode == 1)
signal = s;
end
% calculate STFT frames
minFrameLen = fs*frameSizeMS/1000; 
frameLenPow = nextpow2(minFrameLen);
frameLen = 2^frameLenPow; % frame length = fft size
eval(['frameWindow = ' windowType '(frameLen);']);
[B,F,T] = specgram(signal,frameLen,fs,frameWindow);

%%% WRITE CODE HERE %%%
     % Complete the rest of the code in order to compute the signal
     % spectrogram. Use the Matlab function "specgram", and
     % save the result into vectors [B,F,T] (this will be used
     % later for plotting the time and frequency axes).

B_EDRdb = zeros(size(B));

% calculate the EDR
     %%% WRITE CODE HERE %%%
     % The variable "B_EDRdb" is the EDR in units of decibels.
     % Keep the matrix orientation of the EDR the same as the spectrogram. 
     % Implementation hint: It is faster to calculate the EDR going
     % backward in time rather than forward.

B_EDRdb = 20*log10(B_EDRdb);

% normalize EDR to 0 dB and truncate the plot below a given dB threshold
offset = max(max(B_EDRdb));
B_EDRdbN = B_EDRdb - offset;
[nBins,nFrames] = size(B_EDRdbN);
for i=1:nFrames
  I = find(B_EDRdbN(:,i) < minPlotDB);
  if (I)
    B_EDRdbN(I,i) = minPlotDB;
  end
end

% plot the energy decay relief
figure(gcf);clf;
mesh(T,F/1000,B_EDRdbN);
view(130,30);
title('Normalized Energy Decay Relief (EDR)');
xlabel('Time (s)');ylabel('Frequency (kHz)');zlabel('Magnitude (dB)');
axis tight;zoom on;
