% mv240_to.mlx
% Frequenzanalyse von Luftschallsignalen mittels Oktav-Bandpassfilter und Teiloktav-Bandpassfilter
%
% Dipl.-Ing. Armin Rohnen LbA
% Hochschule München - Munich University of Applied Sciences
% Faktultät 03 - Labor für Schwingungstechnik und Maschinendanmik
% Dachauer Straße 98b, 80335 München
% eMail: rohnen@hm.edu oder armin@rohnen.net
%
% München, 10. März 2022
%
% Lizenzbedingungen:
% Dieses Programm ist freie Software. Sie können es unter den Bedingungen der GNU General Public License, wie von der Free Software Foundation veröffentlicht, weitergeben und/oder modifizieren, entweder gemäß Version 3 der Lizenz oder (nach Ihrer Option) jeder späteren Version.
% Die Veröffentlichung dieses Programms erfolgt in der Hoffnung, daß es Ihnen von Nutzen sein wird, aber OHNE IRGENDEINE GARANTIE, sogar ohne die implizite Garantie der MARKTREIFE oder der VERWENDBARKEIT FÜR EINEN BESTIMMTEN ZWECK. Details finden Sie in der GNU General Public License.
% GNU General Public License http://www.gnu.org/licenses/    http://www.gnu.de/documents/gpl.de.html

clear;
load('20220225_Casablanca.mat');

plot(time, data)
title('Zeitsignal der Luftschallmessung')
ylabel('Schalldruck in Pa')
xlabel('Zeit in s')

% A-Bewertung des Signals

data_A = afilter(data, fs);

SPL = splMeter('FrequencyRange', [22 22050], ...
    'FrequencyWeighting', 'Z-weighting', ...
    'TimeWeighting', 'Fast', ...
    'PressureReference', 20e-6, ...
    'CalibrationFactor', 1, ...
    'SampleRate', fs);

ntel = 6;
filterdefs = fdesign.octave(ntel,'Class 0','N,F0',8, 1000,fs);

tAchse = time(1,1):1/1000:time(end,1);

frequenzen = validfrequencies(filterdefs);
Nfc = length(frequenzen);

for xi=1:Nfc
    filterdefs.F0 = frequenzen(xi);    

% der Pegelverlauf der einzelnen Teiloktaven wird auf eine Zeitauflösung von 1 ms reduziert. Das führt zu einer erheblichen Datenreduktion, ohne dabei das Pegelergebnis zu verändern

    to(:,xi) = interp1(time, SPL(filter(design(filterdefs, 'butter'),data_A)), tAchse);
end

mesh(frequenzen,tAchse, to);
axis([20 10000 0 20 30 90]);
colormap(jet);
colorbar;
caxis([30 90])
title('Teiloktavanalyse - 1/6-Oktaven')
xlabel('Frequenz in Hz')
ylabel('Zeit in s')
zlabel('Schalldruck in dB(A)')
view(15,45)