close all; clear
load akustikdata
fs = 16000;
Ts = 1/fs;
N = length(u);
t = Ts*(1:N)';
% y1 = detrend(y1);
% y2 = detrend(y2);
% u = detrend(u);

%%
figure(1)
plot(t,[u y1 y2])
legend('input','large room', 'small office','FontSize',14)
ylabel('Time[s]')
set(gca,'FontSize',16)

% soundsc(u,fs)
% soundsc(y1,fs)
% soundsc(y2,fs)

%figure(13)
% Hmm, large 60Hz noise 
%uend = u(30000:end);
%[Pxx,f] = pwelch(uend,[],[],[],fs);
%loglog(f,Pxx)

%%
figure(2)
M = ? % Try values in the range 200-3000
na = ? % Try values in the range 100-1000
tv = Ts*(0:M)';
data1 = iddata(y1,u,Ts);
g1hat = cra(data1,M,na,0);
plot([g1(1:M+1)])
grid on
hold on

figure(3)
data2 = iddata(y2,u,Ts);
g2hat = cra(data2,M,na,0);
plot([g2hat g2(1:M+1)])
grid on
hold on

%%
figure(4)
y1hat = filter(g1hat,1,u);
plot([y1 y1hat])
soundsc(y1,fs)
pause(5)
soundsc(y1hat,fs)
SNR1 = 20*log10(std(y1)/std(y1-y1hat))

figure(5)
y2hat = filter(g2hat,1,u);
plot([y2 y2hat])
soundsc(y2,fs)
pause(5) 
soundsc(y2hat,fs)
SNR2 = 20*log10(std(y2)/std(y2-y2hat))




