clear all
%close all
warning off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Extract and concatenate sample paths of raw and filtered data for one channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%trigger = textread('./data_trig/trig_166.txt');
%trigger = textread('./data_trig/R0370/trig_160.txt');

trigger = textread('E:\data_trig\R0376\trig_178.txt');%31.5Hz

%trigger = textread('E:\data_trig\R0380\trig_178.txt');%31.5Hz
%trigger = textread('E:\data_trig\R0380\trig_160.txt');%1.5Hz
%trigger = textread('E:\data_trig\R0380\trig_165.txt');%3.5Hz

%trigger = textread('E:\data_trig\R0380\trig_168.txt');%7.5Hz
%trigger = textread('E:\data_trig\R0380\trig_175.txt');%15.5Hz trig 16

am_freq = 31.5;%15.5;%7.5;%3.5;%1.5;%31.5;
len = length(trigger);
pretrigger = 0;
posttrigger = 0;
index = 0;
stim_duration = 2;
chanx = 66;%45;%117;%100;%120;%66;%120;%66;%104;%
%120;%best%97;%excell%128;%vgood%138;3.5Hz
%97;%122;%1.5Hz
%98;%63;%64;%66;%97;%31.5Hz
%66;%31.5 Hz
channel = [chanx 158 159 160];
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_chan66_11around.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_contourmidline.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_frontal.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_Lcontour.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_midline.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_Rcontour.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_viusal.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_6chanrnd66.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_25_42.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_0_23-filtered-mref_occi12.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_0_23-filtered-mref_front111.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_84_100.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_chan0.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_0_23-filtered-mref_front.sqd';
%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_0_23-filtered-mref_chan66_circ2.sqd';

% srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128.sqd';%Sqd file to be ported
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128.sqd';%Sqd file to be ported

dest = 'E:\sqdfile\R0376 AM 2 sec SFN.10.05.04.run2.AdC300.sqd'

srce = 'E:\sqdfile\R0376 AM 2 sec SFN.10.05.04.run2-filtered_128.sqd'


%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2.sqd';%Sqd file to be ported

%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered-mref-pca256taps.sqd'
%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered-self.sqd';%Sqd file to be ported

%dest = 'F:\Nayef\CorexAntiSymmetry\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref_0_23-filtered-mref_front.sqd';
%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered-mref.sqd'
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_ff99.sqd';%Sqd
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2.sqd';%Sqd file to be ported
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-qrd_lsl.sqd';%Sqd file to be ported
%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-lpf-filtered.sqd'

%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2.sqd';%Sqd file to be ported
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered-svdsc.sqd';%Sqd file to be ported
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-calm.sqd';%Sqd file to be ported

%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered-mref.sqd';%Sqdfile to be ported
%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-qrd_lsl.sqd';%Sqd file to be ported

%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128.sqd';%Sqd file to be ported
%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_256.sqd';%Sqd file to be ported

%compare neural denoising
%srce = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128.sqd';%Sqd file to be ported
%dest = 'E:\sqdfile\R0380.AM.2.sec.SFN.10.6.04.Run2-filtered_128-filtered.sqd';%Sqd file to be ported


% srce = 'E:\sqdfile\R0370.AM.2sec.SFN.9.24.04.Run2.sqd';%Sqd file to be ported
% dest = 'E:\sqdfile\R0370.AM.2sec.SFN.9.24.04.Run2-filtered.sqd';%Sqd file to be ported

%srce = 'E:\sqdfile\R0292.ModBand.1.30.04.Run2.sqd';%Sqd file to be ported
%dest = 'E:\sqdfile\R0292.ModBand.1.30.04.Run2-filtered.sqd';%Sqd file to be ported

%[data,info] = sqdread(srce,'Channels',[0:0],'Samples',[1 1]);%extract info
info = sqdread(srce,'info');%extract info
sf = info.SampleRate;%sampling frequency in KHZ
num_samples = info.SamplesAvailable;%total number of samples
for chan = channel,
    index = index+1;
    data = sqdread(srce,'Channels',[chan-1:chan-1],'Samples',[1 num_samples]);
    D = sqdread(dest,'Channels',[chan-1:chan-1],'Samples',[1 num_samples]);

    %eval(['load ./data_out/data_flms',num2str(chan),'.mat']);
    %eval(['load ./data_in/data',num2str(chan),'.mat data']);
    len_data = length(D);
    data_proc = [];
    data_raw = [];
    
    for j = 1:len,
        lowtrig = trigger(j)*sf - pretrigger*sf + 1;%trigger -(pretrigger) + 1(matlab starts at 1 vs meg160 at 0)
        hightrig = lowtrig + pretrigger*sf + stim_duration*sf + posttrigger*sf -1;%  (pretrigger) + 1s recorded data +  (posttrigger)
        data_proc = [data_proc; D(round(lowtrig):round(hightrig),1)];
        data_raw = [data_raw; data(round(lowtrig):round(hightrig),1)];
    end
    data_filt(:,index) = data_proc;
    data_orig(:,index) = data_raw;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
am_bin = am_freq*stim_duration*len+1;
E = data_filt(:,1);
u = data_orig(1:length(E),2:4);
%u = sqrt(mean(data_orig(:,2:4).^2,2));
fft_E1 = abs(fft(data_orig(:,1))).^2;
fft_E2 = abs(fft(E)).^2;
fft_u2 = mean(abs(fft(u)).^2,2);

snr_filt = 10*log10(sum(fft_E2)./sum(fft_u2));
lb = round(length(fft_E2)*2/sf);
ub = round(length(fft_E2)*4/sf);
snr_2to4hz_filt = 10*log10(sum(fft_E2(lb:ub,:))./sum(fft_u2(lb:ub,:)));
snr_1to10hz_filt = 10*log10(sum(fft_E2(lb:2.5*ub,:))./sum(fft_u2(lb:2.5*ub,:)));
snr_2to20hz_filt = 10*log10(sum(fft_E2(lb:5*ub,:))./sum(fft_u2(lb:5*ub,:)));
lb = round(length(fft_E2)*175/sf);
ub = round(length(fft_E2)*185/sf);
snr_175to185hz_filt = 10*log10(sum(fft_E2(lb:ub,:))./sum(fft_u2(lb:ub,:)));
snr_am_bin_filt = 10*log10(fft_E2(am_bin)./fft_u2(am_bin,:));

snr_raw = 10*log10(sum(fft_E1)./sum(fft_u2));
lb = round(length(fft_E1)*2/sf);
ub = round(length(fft_E1)*4/sf);
snr_2to4hz_raw = 10*log10(sum(fft_E1(lb:ub,:))./sum(fft_u2(lb:ub,:)));
snr_1to10hz_raw = 10*log10(sum(fft_E1(lb:2.5*ub,:))./sum(fft_u2(lb:2.5*ub,:)));
snr_2to20hz_raw = 10*log10(sum(fft_E1(lb:5*ub,:))./sum(fft_u2(lb:5*ub,:)));
lb = round(length(fft_E1)*175/sf);
ub = round(length(fft_E1)*185/sf);
snr_175to185hz_raw = 10*log10(sum(fft_E1(lb:ub,:))./sum(fft_u2(lb:ub,:)));
snr_am_bin_raw = 10*log10(fft_E1(am_bin)./fft_u2(am_bin,:));


data_orig_fft = fft(data_orig(:,1));
data_filt_fft = fft(data_filt(:,1));
orig_am_bin = [abs(data_orig_fft(am_bin)),180/pi*angle(data_orig_fft(am_bin))];%amplitude/phase before filtering
filt_am_bin = [abs(data_filt_fft(am_bin)),180/pi*angle(data_filt_fft(am_bin))];%amplitude/phase after filtering

infoc.nepoch = len;
infoc.chan = chanx;
infoc.am = am_freq;

infoc.snr_raw = snr_raw;
infoc.snr_2to4hz_raw = snr_2to4hz_raw;
infoc.snr_1to10hz_raw = snr_1to10hz_raw;
infoc.snr_2to20hz_raw = snr_2to20hz_raw;
infoc.snr_175to185hz_raw = snr_175to185hz_raw;
infoc.snr_am_bin_raw = snr_am_bin_raw;

infoc.snr_filt = snr_filt;
infoc.snr_2to4hz_filt = snr_2to4hz_filt;
infoc.snr_1to10hz_filt = snr_1to10hz_filt;
infoc.snr_2to20hz_filt = snr_2to20hz_filt;
infoc.snr_175to185hz_filt = snr_175to185hz_filt;
infoc.snr_am_bin_filt = snr_am_bin_filt;

infoc.orig_am_bin = orig_am_bin;
infoc.filt_am_bin = filt_am_bin;

infoc
%plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = ([1:length(data_filt)]-1)/length(data_filt)*sf;

% %Black and white Paper
% figure;
% plot(f,abs(fft(data_orig(:,1))),'Color',[0.5 0.5 0.5]);
% %plot(f,abs(fft(data_orig(:,1))));
% hold
% plot(f,abs(fft(data_filt(:,1))),'k');
% axis([0 200 -inf 800000]);
% xlabel('Hz');

% figure;
% subplot(2,1,1);plot(f,abs(fft(data_orig(:,1))),'k');
% axis([0 200 -inf 800000]);
% title('DFT Raw Neuronal Channel ');
% xlabel('Hz');
% subplot(2,1,2);plot(f,abs(fft(data_orig(:,2:4))),'k');
% axis([0 200 -inf 800000]);
% xlabel('Hz');
% title('Mean Square DFT Reference Channels');

figure;
plot(f,abs(fft(data_orig(:,1))));
%plot(f,abs(fft(data_orig(:,1))));
hold
plot(f,abs(fft(data_filt(:,1))),'r');
axis([0 200 -inf 800000]);
xlabel('Hz');
legend('original data','filtered data');

figure;
plot(f,abs(fft(data_orig(:,1))),'b');
%plot(f,abs(fft(data_orig(:,1))));
hold
plot(f,abs(fft(data_filt(:,1))),'r');
axis([0 60 -inf 800000]);
xlabel('Frequency (Hz)');
ylabel('PSD');
legend('original data','filtered data');

% figure;
% plot(f,abs(fft(data_orig(:,1))),'b');
% %plot(f,abs(fft(data_orig(:,1))));
% hold
% plot(f,abs(fft(data_filt(:,1))),'r');
% axis([14 17 -inf 400000]);

figure;
plot(f,abs(fft(data_orig(:,1))));
%plot(f,abs(fft(data_orig(:,1))));
hold
plot(f,abs(fft(data_filt(:,1))),'r');
axis([0 10 -inf 800000]);
xlabel('Hz');

figure;
bar(f,abs(fft(data_filt(:,1))));
axis([am_freq-1.20 am_freq+1.20 -inf 200000]);
xlabel('Hz');

figure;
%plot(f,abs(fft(data_orig(:,1))),'Color',[0.5 0.5 0.5]);
plot(f,abs(fft(data_filt(:,1))),'r');
hold
plot(f,abs(fft(data_orig(:,1))));
axis([0 200 -inf 800000]);
xlabel('Hz');
legend('filt. data','original data');

figure;
plot(f,abs(fft(data_filt(:,1))));
axis([0 60 -inf 560000]);
xlabel('Hz');

figure;
subplot(2,1,1);plot(f,abs(fft(data_orig(:,1))));
axis([0 200 -inf 800000]);
%title('DFT Raw Neuronal Channel ');
xlabel('Hz');
subplot(2,1,2);plot(f,abs(fft(data_orig(:,2:4))));
axis([0 200 -inf 800000]);
xlabel('Hz');
%title('Mean Square DFT Reference Channels');


pmtm_coef= 3/2;
[pm_rawdata F] = pmtm(data_orig(:,1),pmtm_coef,[],sf); 
[pm_filtdata F] = pmtm(data_filt(:,1),pmtm_coef,[],sf); 
for k = 1:3,
    [pm_noise(:,k) F] = pmtm(data_filt(:,1+k),pmtm_coef,[],sf); 
end

[psdvalOrig,psdfreq]=psd(data_orig(:,1),2048*4,500);
[psdvalFilt,psdfreq]=psd(data_filt(:,1),2048*4,500);
figure;plot(psdfreq,20*log10(psdvalFilt),'r',psdfreq,20*log10(psdvalOrig),'b')


% figure;
% subplot(3,1,1);semilogy(F,pm_rawdata);
% axis([0 200 -inf inf]);
% title('pmtm raw data');
% xlabel('Hz');
% ylabel('db');
% subplot(3,1,2);semilogy(F,pm_noise);
% axis([0 200 -inf inf]);
% title('pmtm noise');
% xlabel('Hz');
% ylabel('db');
% subplot(3,1,3);semilogy(F,pm_filtdata);
% axis([0 200 -inf inf]);
% title('pmtm error');
% xlabel('Hz');
% ylabel('db');
% 
% %Test Read/Write filtering successful
% info = sqdread(srce,'info');%extract info
% num_samples = info.SamplesAvailable;%total number of samples
% samp = num_samples;%1000;%
% datar = sqdread(srce,'Channels',[0:191],'Samples',[samp samp]);
% dataw = sqdread(dest,'Channels',[0:191],'Samples',[samp samp]);
% [datar' dataw']