clear all
close all
clc

%% Create filters
fs=48e3;
N = 16; % filter order
ThirdOct = octaveFilter(1000,'1/3 octave','SampleRate',fs);
F0 = getANSICenterFrequencies(ThirdOct);
F0(F0<100) = [];
F0(F0>16e3) = [];
Nfc = length(F0);

for i=1:Nfc
    oneThirdOctaveFilterBank{i} = octaveFilter('FilterOrder', N, ...
        'CenterFrequency', F0(i), 'Bandwidth', '1/3 octave',...
        'SampleRate', fs);
end


%% Init of PlayRec
AutoPlayrecInit(fs)

%% Generate Pink Noise
T=2;
pinkDuration = fs*T;
pinkNoise = dsp.ColoredNoise(1,pinkDuration,1);
rng default;
pinkPlay = pinkNoise();
stimuli=pinkPlay.*0.5;

%% Play pink noise and record response
out=playrec('playrec',[stimuli stimuli],[1 2],length(stimuli)+1500,1);
playrec('block',out);
x = playrec('getRec',out);

pause(1);

%% Apply 1/3 octave band filters on signal x
for i=1:Nfc
        L(:,i) = oneThirdOctaveFilterBank{i}(x);
end

%% Convert to dB SPL re 20 uPa
LdB = 20*log10(rms(L)/20e-6);

LdB'