%
% Beispiel zur Nutzung von ntelOktave
%
%
%   Dipl.-Ing. Armin Rohnen LbA
%   Labor für Maschinendynamik
%   Hochschule für angewandte Wissenschaften München
%   University of applied sciences Munich
%
%   05. Januar 2016
%
%
%%
% Parameter die als Übergabeparameter erforderlich sind

clear;
datetime

%filename = '15ccw_01.wav';
%[data, fs] = audioread(filename);
load('spieluhr.mat');
ntel = 12;                       % ntelOktav
filterord = 1000;               % Steilheit der ntelOktavfilter 1000 = 6000 dB / Oktave
tdata = OData(:,1);             % Zeitsignal für die ntelOktav-Analyse
zeitk = 25;                     % Zeitkonstante in ms für die RMS-Pegelberechnung
overlap = 50;                   % Overlap für die RMS-Pegelberechnung 0 = kein Overlap

[dataOut, timeOut, mittenfrequenzen, frequenzen] = ntelOktave(ntel, filterord, fs, tdata, zeitk, overlap);

datetime

%%
%
%   Darstellung als "Treppenstufen" 
%   Diese entstehen, wenn für fu und fo des jeweiligen Frequenzbereichs
%   der Pegelwert der Mittenfrequenz abgelegt wird.
%   fu des "neuen" Frequenzbereichs ist 0.00000001 Hz größer als fo des
%   vorherigen.
%   Die Grenzfrequenz für den Treppensprung wird in die Mitte der beiden
%   Frequenzbereiche gelegt
%

anz_f = length(mittenfrequenzen);
anz_time = length(timeOut);

%
% neue Frequenzachse erstellen
%

new_f_point = 0;
for k=1:anz_f
    
    if k==1
        fu = frequenzen(k,1);
    else
        fu = fo + 0.00000001;
    end;
    
    if k<anz_f
        fo = (frequenzen(k,2)+frequenzen(k+1,1))/2;
    else
        fo = frequenzen(k,2);
    end;
    
    new_f_point = new_f_point +1;
    fAxis(1,new_f_point) = fu;
    new_f_point = new_f_point +1;
    fAxis(1,new_f_point) = mittenfrequenzen(1,k);
    new_f_point = new_f_point +1;
    fAxis(1,new_f_point) = fo;
end;
    
for n=1:anz_time
    new_f_point = 0;
    for k=1:anz_f
        new_f_point = new_f_point +1;
        cornerData(n,new_f_point) = dataOut(n,k);
        new_f_point = new_f_point +1;
        cornerData(n,new_f_point) = dataOut(n,k);
        new_f_point = new_f_point +1;
        cornerData(n,new_f_point) = dataOut(n,k);
    end;
end;

mesh(fAxis, timeOut, cornerData);
set(gca,'XScale','log');
axis([100 10000 0 25 0.1 1])
colorbar;
view(60,30);