clear all
close all 

run /users/swinkels/deploy/MatlabVirgoTools/trunk/startup.m

gps_s=1187963502;
gps_calm=1187965203 ;
dur=100; 
fs = 20e3; %frequenza di campionamento
nav = 100; %medie
npoints = fs*dur;
npt=2^floor(log(2*npoints/nav)/log(2)); 

data_n1='Hrec_hoft_20000Hz';
data_n2='ENV_EDB_MIC';

datan1=getChannel(data_n1,gps_s,dur);
datan1_calm=getChannel(data_n1,gps_calm,dur);
datan2=getChannel(data_n2,gps_s,dur);
datan2_calm=getChannel(data_n2,gps_calm,dur);

data1=decimate(datan1,1); 
data1_calm=decimate(datan1_calm,1); 
data2=decimate(datan2,1); 
data2_calm=decimate(datan2_calm,1); 

[d1,f] = pwelch(data1,hann(npt),[],npt,fs);
[d1_calm,f] = pwelch(data1_calm,hann(npt),[],npt,fs);
[d2,f] = pwelch(data2,hann(npt),[],npt,fs);
[d2_calm,f] = pwelch(data2_calm,hann(npt),[],npt,fs);
data_1=sqrt(d1);
data_1_calm=sqrt(d1_calm);
data_2=sqrt(d2);
data_2_calm=sqrt(d2_calm);

noise= data_2_calm.*data_1./data_2;
upp_lim= noise;
diff=d1./d1_calm;

mains_line = [50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450];
idx = diff>= 2.5;   
IDX = diff<2.5;

upp_lim(idx)=NaN;
noise(IDX)=NaN;

idx_mains = ones(1,length(f));
for k =1 :length(mains_line)
    [g,i] = min(abs(f-mains_line(k)));
    idx_mains(i) = 0;
end
idx_mains = logical(idx_mains);
idx = idx & idx_mains';

figure(1)
subplot(3,1,1)
loglog(f,d1, f, d1_calm )
xlabel('freq [Hz]')
ylabel(sprintf('%s',data_n1 ))
xlim([30 2000])
subplot(3,1,2)
loglog(f,d2, f, d2_calm )
xlabel('freq [Hz]')
ylabel(sprintf('%s',data_n2 ))
xlim([30 2000])
subplot(3,1,3)
loglog(f, data_2,f(idx), data_2(idx),'.')
xlabel('freq [Hz]')
xlim([30 2000])

figure(3)
loglog(f,data_1_calm,f, upp_lim,'b',f,noise,'r')
grid on
ylabel('[1/sqrt(Hz)]')
xlabel('Frequency [Hz]')
legend('Hrec (quiet)','proj: ENV EDB MIC upper limit','proj.: ENV EDB MIC ')
title(' DET LAB Acoustic noise (0.001) projection on Hrec');
xlim([30 2000])


% fileID = fopen('%s.txt',data_n2,'w');
% fprintf(fileID,'freq Hrec_calm \n \n %f %f', f, data_1_calm);
% fclose(fileID);

data=[f, data_1_calm, data_1, data_2_calm, data_2]
csvwrite('csvlist.dat',data)


