
% Systemidentifikation, Basic 

%% define parameters
playdev=2;  % ID of the output device, see audiodevinfo
recdev=0; % ID of the input device

fs=48000; % sampling rate (Hz)
reverbExp=round(1*fs);  % expected reverberation trail
preOff=round(0.05*fs);  % expected latency  

f1=50;  % Start-Frequenz des Sweeps (Hz)
f2=18000; % End-Frequenz des Sweeps (Hz)
T=1; % Dauer des Sweeps (s)
fadeIn=round(0.01*fs); % fade-in des Sweeps (samples)
fadeOut=round(0.01*fs); % fade-out des Sweeps (samples)

%% check devices
disp(['Play: ' audiodevinfo(0,playdev)]);
disp(['Record: ' audiodevinfo(1,recdev)]);

%% create figure
figure('Position',[0 0 1200 800]);

%% create excitation signal
[sw,isw]=expsweep(f1,f2,T,fs);
sw=fade(sw,length(sw)+reverbExp, fadeIn , fadeOut);

subplot(2,3,1);
specgram(sw,256,fs);
title('Excitation signal');

%% Play and record
p=audioplayer(sw,fs,24,playdev);
r=audiorecorder(fs,24,1,recdev);
record(r);
playblocking(p);
stop(r);
y=getaudiodata(r);

subplot(2,3,2);
plot(y);
set(gca,'Ylim',[-1.1 1.1]);
title('Recorded signal');

subplot(2,3,3);
specgram(y,256,fs);
title('Recorded signal');

%% System identification
nges=length(sw)+length(isw);
h=real(ifft(fft(y,nges).*fft(isw,nges)));
h0=real(ifft(fft(sw,nges).*fft(isw,nges)));

subplot(2,3,4);
tax=0:1/fs:(length(h0)-1)/fs;
plot(tax, etc(h0));
hold on; box on;
plot(tax, etc(h),'r');
xlabel('Samples');
title('Impulse responses');
legend('control','measured');

%% Latency calculation
[a0,t0]=max(etc(h0));
[a,t]=max(etc(h));
lat=(t-t0)/fs;
att=a0-a;

%% crop to the reverb
h0s=h0(round(T*fs)-preOff:t0+reverbExp);
hs=h(round(T*fs)-preOff:t0+reverbExp);

subplot(2,3,5); 
hold on; box on;
tax=0:1/fs:(length(hs)-1)/fs;
plot(tax, etc(h0s));
plot(tax, etc(hs),'r');
xlabel('Time (s)');
title('Cropped impulse responses');
legend('control','measured');

subplot(2,3,6);
spect(fft(h0s),fs,'a');
spect(fft(hs),fs,'a','r');
title('Amplitude spectra');
legend('control','measured','Location','SouthWest');

%% display results
disp(['Latency: ' num2str(lat) ' s']);
disp(['Attenuation: ' num2str(att) ' dB']);
