% CumulativeSpectralDecay(fs,imp)
%
% fs   - samplingsfrekvens
% imp  - impulssvaret (rad- eller kolonnvektor)
%
% Plottar impulssvaret imp som vattenfall (Cumulative Spectral Decay).
%
% L-J Brännmark, Signals & Systems, 2008

function CumulativeSpectralDecay(fs,imp);

imp=imp(:).';

f_lo=20; f_hi=fs/2;    % Frekvensgränser för vad som visas i mesh-plotten.
numoffreqs=320;        % Mesh-plottens upplösning i frekvens.
numoftimesteps=120;    % Mesh-plottens upplösning i tid.
impsize=size(imp);
implen=impsize(2);     % Impulslängden i samples.

displaylen=round(min([0.25*fs implen])); % Visa högst 0.25 sekunder.

fftlength=max([16384 implen]);           % Antal fft-punkter.
if mod(fftlength,2)~=0
    fftlength=fftlength+1;
end

timestep=floor(displaylen/numoftimesteps);       % Tiddsteg i plotten
t=(1/fs)*linspace(0,displaylen-1,numoftimesteps);% Tidsvektor till plotten.

freqs=logspace(log10(f_lo),log10(f_hi),numoffreqs); % Skapa logaritmiskt fördelade
                                                 % frekvenser, f_lo-f_hi [Hz].
nyqvistindex=floor(fftlength/2+1);               % Index som motsvarar fs/2 i fft:n
freqindex=ceil( (freqs/(fs/2) )*nyqvistindex);   % Index motsvarande logaritmiska
                                                 % frekvensskalan.

M=zeros(numoftimesteps,fftlength);
for k = 1:numoftimesteps
    M(k,:) = abs( fft( imp( (k-1)*timestep+1 : end ),fftlength ) ).^2;
end

M_plot=10*log10(abs( M(:,freqindex) ) );         % Behåll bara de tids- och frekvensindex
                                                 % som skall plottas
M_plot(find(M_plot<-30))=-30;                    % Nivågolv =-30 dB för vad som
                                                 % skall visas.

mesh(freqs,t,M_plot);                            % Plotta "vattenfallet"
shading interp;
set(gca,'xscale','log','xdir','reverse');        % Sätt logaritmisk frekvensskala. 
                                                 % och vänd på frekvensaxeln.
axis([f_lo f_hi 0 1/fs*(displaylen+20) -30 10]);view(-159,68);
xlabel('Frequency [Hz)]');ylabel('Time [s]');zlabel('Amplitude [dB]');
title('Cumulative Spectral Decay');