diff --git a/example/KHU_example.m b/example/KHU_example.m index 46aea67..f8a02ab 100644 --- a/example/KHU_example.m +++ b/example/KHU_example.m @@ -7,11 +7,11 @@ %load boundary voltages of baseline recording, cleaning the data and %plotting -[BV_baseline, KHUs_baseline]=KHU_Load('..\resources\data\KHU\B1',1,1); +[BV_baseline, KHUs_baseline]=KHU_Load('../resources/data/KHU/B1',1,1); %% load Perturbation %load boundary voltages of baseline recording, cleaning the data but not %plotting -[BV_perturbation, KHUs_perturbation]=KHU_Load('..\resources\data\KHU\P1',1,0); +[BV_perturbation, KHUs_perturbation]=KHU_Load('../resources/data/KHU/P1',1,0); diff --git a/example/ScouseTom_example.m b/example/ScouseTom_example.m index 5aa028f..38bdf4f 100644 --- a/example/ScouseTom_example.m +++ b/example/ScouseTom_example.m @@ -7,14 +7,14 @@ %load boundary voltages of resistor phantom recording. The log files %generated by the ScouseTom ssytem are located in same directory -[BVstruc_TD]=ScouseTom_Load('..\resources\data\ScouseTom\forreal_2_clean.bdf'); +[BVstruc_TD]=ScouseTom_Load('../resources/data/ScouseTom/forreal_2_clean.bdf'); %% load Multifrequency data %load boundary voltages of resistor phantom recording using multiple %frequecies. The log files generated by the ScouseTom system is located in %the same directory -[BVstruc_MF]=ScouseTom_Load('..\resources\data\ScouseTom\MultiFreq.bdf'); +[BVstruc_MF]=ScouseTom_Load('../resources/data/ScouseTom/MultiFreq.bdf'); %% load contact impedance check @@ -22,4 +22,4 @@ % electrodes. This code is a bit different as we want to estimate the % contact impedance and plot it quickly. The protocol is also independent % of that set by the ExpSetup -[Zout]=ScouseTom_Load('..\resources\data\ScouseTom\S1a_Z4.bdf'); +[Zout]=ScouseTom_Load('../resources/data/ScouseTom/S1a_Z4.bdf'); diff --git a/example/SwissTom_example.m b/example/SwissTom_example.m index 4ad9676..53e8d02 100644 --- a/example/SwissTom_example.m +++ b/example/SwissTom_example.m @@ -7,4 +7,4 @@ %load boundary voltages of resistor phantom recording. Removing injected %channels, saturated channels, and plotting data -[STout]=SwissTom_Load('..\resources\data\Swisstom\r-phantom.eit',1,1,1); +[STout]=SwissTom_Load('../resources/data/Swisstom/r-phantom.eit',1,1,1); diff --git a/src/matlab/KHU_Load.m b/src/matlab/KHU_Load.m index 65a4fe4..6186da8 100644 --- a/src/matlab/KHU_Load.m +++ b/src/matlab/KHU_Load.m @@ -359,7 +359,7 @@ KHU.info=info; %get the name for the new files -k=strfind(dirname,filesep); +k=strfind(dirname,'/'); newnamestr=dirname(k(end)+1:end); diff --git a/test/S3b_MF1_log.mat b/test/S3b_MF1_log.mat deleted file mode 100644 index 220b888..0000000 Binary files a/test/S3b_MF1_log.mat and /dev/null differ diff --git a/tests/Test_ScouseTom.m b/tests/Test_ScouseTom.m new file mode 100644 index 0000000..ef17743 --- /dev/null +++ b/tests/Test_ScouseTom.m @@ -0,0 +1,163 @@ +% test code is functioning on real examples. Requires files from https://doi.org/10.5281/zenodo.4783547 + +datadir='I:/Load_data_testing/'; +%datadir='data/'; +%% load bdfs + +somethingwentwrong=0; + + +try + BVs=ScouseTom_Load([datadir 'Baseline1.bdf']); +catch err + fprintf(2,'error on initial test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'Baseline1_missingstart.bdf']); +catch err + fprintf(2,'error on missing start test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'Baseline1_missingend.bdf']); +catch err + fprintf(2,'error on missing end test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'Baseline1_missingboth.bdf']); +catch err + fprintf(2,'error on missing both test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); +%% multifreq + +try + BVs=ScouseTom_Load([datadir 'MultiFreq.bdf']); +catch err + fprintf(2,'error on Multifreq test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'MultiFreq_missingstart.bdf']); +catch err + fprintf(2,'error on Multifreq missing start test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'MultiFreq_missingend.bdf']); +catch err + fprintf(2,'error on Multifreq missing end test\n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'MultiFreq_missingboth.bdf']); +catch err + fprintf(2,'error on Multifreq missing both test\n'); + somethingwentwrong=1; +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +%% Z check +try + Z=ScouseTom_Load([datadir 'firstinnir_zcheck_1.bdf']); +catch err + fprintf(2,'error on Z load \n'); + somethingwentwrong=1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + + +%% AC check + +try + BVs=ScouseTom_Load([datadir 'rptest2.eeg']); +catch err + fprintf(2,'error on Actichamp test 1\n'); + somethingwentwrong=1; +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + +try + BVs=ScouseTom_Load([datadir 'rptest3.eeg']); +catch err + fprintf(2,'error on Actichamp test 2\n'); + somethingwentwrong=1; +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + + +%% Big file + +try + BVs=ScouseTom_Load([datadir 'RealThing_NIRFACE_2.bdf']); +catch err + fprintf(2,'error on Multifreq big file test\n'); + somethingwentwrong=1; +end + +disp('================================================='); +disp('================================================='); +disp('================================================='); + + + +%% did it work? + +if somethingwentwrong + fprintf(2,'Oh No! SOMETHING WENT WRONG!\n'); +else + fprintf('YAY! IT WORKED :)\n'); +end + diff --git a/tests/data/readme.md b/tests/data/readme.md new file mode 100644 index 0000000..e446dff --- /dev/null +++ b/tests/data/readme.md @@ -0,0 +1,2 @@ +# Test dataset +Put files from [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4783527.svg)](https://doi.org/10.5281/zenodo.4783527) in this directory to run tests diff --git a/tests/filters/check_acc.m b/tests/filters/check_acc.m new file mode 100644 index 0000000..fb08a3b --- /dev/null +++ b/tests/filters/check_acc.m @@ -0,0 +1,121 @@ +function [ Amp_error, Phase_error,Vsig,Vsigdemod,Filt,trim_demod] = check_acc( Fc,InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetInj, Padsec, Fs,decimate_factor) + %UNTITLED Summary of this function goes here + % Detailed explanation goes here + + %% + inj = [ 2 4]; + + if exist('DCoffset','var') == 0 || isempty(DCoffset) + DCoffset=0; + end + + if exist('DCoffsetInj','var') == 0 || isempty(DCoffsetInj) + DCoffsetInj=0; + end + + + BW=100; + chn=5; + + if exist('Padsec','var') == 0 || isempty(Padsec) + Padsec=4; + end + + if exist('Fs','var') == 0 || isempty(Fs) + Fs=16384; + end + + if exist('decimate_factor','var') == 0 || isempty(decimate_factor) + decimate_factor=0; + end + + + %% Create ideal values, and voltages + + AmpActual= repmat(Amp_Meas,chn,1); + AmpActual(inj)=Amp_Inj; + + a=MeasPhaseDiff; + + %force normal valus of deg + a = mod(a,360); % [0 2pi) + + % shift + jj = a > 180; + a(jj) = a(jj) - 360; + jj = a < 0 - 180; + a(jj) = a(jj) + 360; + + MeasPhaseDiff_corr = a; + + MeasPhase = InjPhase + MeasPhaseDiff_corr; %deg + PhaseActual= repmat(MeasPhaseDiff_corr,chn,1); + PhaseActual(inj)=0; + + Totaltime=ceil(InjTime + 2*Padsec); + t = 0:1/Fs:Totaltime-1/Fs; + %make sin wave + + v_m= Amp_Meas*sin(2*pi*Fc*t+(pi*MeasPhase/180))+DCoffset; + + % change amplitude + v_i=Amp_Inj*sin(2*pi*Fc*t+(pi*InjPhase/180))+DCoffsetInj; + + + %% + V=repmat(v_m,chn,1); + + V(inj(1),:) = v_i; + V(inj(2),:) = v_i; + + V=V'; + + %% decimate now if we want to + + if decimate_factor + for iChn = 1 :size(V,2) + Vtmp(:,iChn) = decimate(V(:,iChn),decimate_factor); + end + Fs=Fs/decimate_factor; + V=Vtmp; + end + + %% + + + %pad with a second of data either side, so the hilbert is more realistic + datastart = round(Padsec*Fs); + dataend = round((Padsec+InjTime)*Fs); + + V(1:datastart,:)=V(1:datastart,:)*0.1; + V(dataend:end,:)=V(dataend:end,:)*0.1; + + InjectionWindows =[datastart dataend]; + + + %% + %find the corresponding filter settings + [trim_demod,Filt,Fc_found]=ScouseTom_data_GetFilterTrim(V(datastart:dataend,inj(1)),Fs,[],[],[]); + + Vsig = V(datastart:dataend,inj(1)); + + [ Vdata_demod,Pdata_demod ] = ScouseTom_data_DemodHilbert( V,Filt); + Vsigdemod = Vdata_demod(datastart:dataend,inj(1)); + + %% + [Vmag,PhaseRaw,VmagSTD,PhaseRawSTD]=ScouseTom_data_getBV(Vdata_demod,Pdata_demod,trim_demod,InjectionWindows); + % [Vmag,PhaseRaw,VmagSTD,PhaseRawSTD]=ScouseTom_data_getBVChunk(Datafilt,trim_demod,InjectionWindows); + + + [Phase]=ScouseTom_data_PhaseEst(PhaseRaw,inj,1); + + Phase_deg=(180/pi) * (Phase); + + Phase_error = PhaseActual - Phase_deg'; + + Amp_error = AmpActual - Vmag'; + + fprintf('Amp error : %.6f, Phase error : %.6f\n',mean(Amp_error),mean(Phase_error)); + end + + \ No newline at end of file diff --git a/tests/filters/test_accuracy.m b/tests/filters/test_accuracy.m new file mode 100644 index 0000000..6251ee1 --- /dev/null +++ b/tests/filters/test_accuracy.m @@ -0,0 +1,811 @@ +% script to test the accuracy of the processing code + +Failed = 0; +Amp_tolerance = 5e-3; %how much error is ok for amp use about 10e-5 smaller than amp +Phase_tolerance = 5e-3; + +plotflag = 1; + +%% test 1 +% fixed amplitude, time and freq. Sweep through phase + +disp('Test 1'); + +Fc = 1800; + +InjTime=2; + +Amp_Inj = 500; +Amp_Meas = 150; + +InjPhase = -95:20:95; +InjPhaseNum = size(InjPhase,2); +MeasPhaseDiff = -355:40:355; +PhaseNum = size(MeasPhaseDiff,2); + +Amp_error_tot=zeros(InjPhaseNum,PhaseNum); +Phase_error_tot=zeros(InjPhaseNum,PhaseNum); + +for iInjPhase = 1:InjPhaseNum + curInjPhase = InjPhase(iInjPhase); + for iPhase = 1:PhaseNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc,InjTime,Amp_Inj,Amp_Meas,curInjPhase,MeasPhaseDiff(iPhase) );'); + Amp_error_totT1(iInjPhase,iPhase) = mean(Amp_error); + Phase_error_totT1(iInjPhase,iPhase) = mean(Phase_error); + + end +end +%% +if plotflag + figure + + subplot(2,1,1) + plot(MeasPhaseDiff,Amp_error_totT1) + ylabel('Error'); + xlabel('Phase difference beetween inj and meas'); + title('T1: Amp error for different starting phase'); + + subplot(2,1,2) + plot(MeasPhaseDiff,Phase_error_totT1) + ylabel('Error'); + xlabel('Phase difference'); + title('T1: Phase error for different starting phase'); + +end + +Test1_AmpError=max(max(abs(Amp_error_totT1))); +Test1_PhaseError=max(max(abs(Phase_error_totT1))); + + +try + assert( Test1_AmpError < Amp_tolerance, 'Test1 Amp error failed'); +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test1_PhaseError < Phase_tolerance, 'Test1 Phase error failed'); +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +%% test 2 + +clear Amp_error_tot Phase_error_tot + +%changing amplitude +disp('Test 2'); + +InjTime=2; +Fc = 1800; +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +Amp_sc=0.5:0.5:10; +AmpNum = size(Amp_sc,2); + +for iAmp = 1:AmpNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc,InjTime,Amp_Inj*Amp_sc(iAmp),Amp_Meas*Amp_sc(iAmp),InjPhase,MeasPhaseDiff);'); + Amp_error_totT2(iAmp) = mean(Amp_error); + Phase_error_totT2(iAmp) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Amp_Meas*Amp_sc,Amp_error_totT2) + ylabel('Error'); + xlabel('Amplitude'); + title('T2: Amp error for different Amplitudes'); + + subplot(2,1,2) + plot(Amp_Meas*Amp_sc,Phase_error_totT2) + ylabel('Error'); + xlabel('Amplitude'); + title('T2: Phase error for different Amplitudes'); + +end + +Test2_AmpError=max(max(abs(Amp_error_totT2))); +Test2_PhaseError=max(max(abs(Phase_error_totT2))); + + +try + assert( Test2_AmpError < Amp_tolerance, 'Test2 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test2_PhaseError < Phase_tolerance, 'Test2 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +%% test 3 + +clear Amp_error_tot Phase_error_tot + +%changing freq - large fixed inj time +disp('Test 3 - large fixed injection time 10s'); + +InjTime=10; + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff);'); + Amp_error_totT3(iFreq) = mean(Amp_error); + Phase_error_totT3(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT3) + ylabel('Error'); + xlabel('Frequency'); + title('T3: Amp error 10s injection'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT3) + ylabel('Error'); + xlabel('Frequency'); + title('T3: Phase error 10s injection'); + +end + +Test3_AmpError=max(max(abs(Amp_error_totT3))); +Test3_PhaseError=max(max(abs(Phase_error_totT3))); + +try + assert( Test3_AmpError < Amp_tolerance, 'Test3 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test3_PhaseError < Phase_tolerance, 'Test3 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + + +%% test 4 + +clear Amp_error_tot Phase_error_tot + +%changing freq - large fixed inj time +disp('Test 4 - small fixed injection time 0.5s'); + +InjTime=.5; + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff);'); + Amp_error_totT4(iFreq) = mean(Amp_error); + Phase_error_totT4(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT4) + ylabel('Error'); + xlabel('Frequency'); + title('T4: Amp error 0.5s injection'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT4) + ylabel('Error'); + xlabel('Frequency'); + title('T4: Phase error 0.5s injection'); + +end + +Test4_AmpError=max(max(abs(Amp_error_totT4))); +Test4_PhaseError=max(max(abs(Phase_error_totT4))); + +try + assert( Test4_AmpError < Amp_tolerance, 'Test4 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test4_PhaseError < Phase_tolerance, 'Test4 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +%% test 5 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 32 +disp('Test 5 - fixed cycles of 32'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 32; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff);'); + Amp_error_totT5(iFreq) = mean(Amp_error); + Phase_error_totT5(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT5) + ylabel('Error'); + xlabel('Frequency'); + title('T5: Amp error 32 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT5) + ylabel('Error'); + xlabel('Frequency'); + title('T5: Phase error 32 cycles'); + +end + +Test5_AmpError=max(max(abs(Amp_error_totT5))); +Test5_PhaseError=max(max(abs(Phase_error_totT5))); + +try + assert( Test5_AmpError < Amp_tolerance, 'Test5 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test5_PhaseError < Phase_tolerance, 'Test5 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +%% test 6 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 3 +disp('Test 6 - fixed cycles of 3'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 3; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff);'); + Amp_error_totT6(iFreq) = mean(Amp_error); + Phase_error_totT6(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT6) + ylabel('Error'); + xlabel('Frequency'); + title('T6: Amp error 3 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT6) + ylabel('Error'); + xlabel('Frequency'); + title('T6: Phase error 3 cycles'); + +end + +Test6_AmpError=max(max(abs(Amp_error_totT6))); +Test6_PhaseError=max(max(abs(Phase_error_totT6))); + +try + assert( Test6_AmpError < Amp_tolerance, 'Test6 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test6_PhaseError < Phase_tolerance, 'Test6 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + + +%% test 7 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 64 +disp('Test 7 - fixed cycles of 64'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 64; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff);'); + Amp_error_totT7(iFreq) = mean(Amp_error); + Phase_error_totT7(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT7) + ylabel('Error'); + xlabel('Frequency'); + title('T7: Amp error 64 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT7) + ylabel('Error'); + xlabel('Frequency'); + title('T7: Phase error 64 cycles'); + +end + +Test7_AmpError=max(max(abs(Amp_error_totT7))); +Test7_PhaseError=max(max(abs(Phase_error_totT7))); + +try + assert( Test7_AmpError < Amp_tolerance, 'Test7 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test7_PhaseError < Phase_tolerance, 'Test7 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +%% test 8 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 128 +disp('Test 8 - fixed cycles of 128'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 128; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff);'); + Amp_error_totT8(iFreq) = mean(Amp_error); + Phase_error_totT8(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT8) + ylabel('Error'); + xlabel('Frequency'); + title('T8: Amp error 128 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT8) + ylabel('Error'); + xlabel('Frequency'); + title('T8: Phase error 128 cycles'); + +end + +Test8_AmpError=max(max(abs(Amp_error_totT8))); +Test8_PhaseError=max(max(abs(Phase_error_totT8))); + +try + assert( Test8_AmpError < Amp_tolerance, 'Test8 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test8_PhaseError < Phase_tolerance, 'Test8 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + + +%% test 9 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 128 +disp('Test 9 - fixed cycles of 128 with dc offset'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 128; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; +DCoffset = 100; +DCoffsetinj = 400; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj );'); + Amp_error_totT9(iFreq) = mean(Amp_error); + Phase_error_totT9(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT9) + ylabel('Error'); + xlabel('Frequency'); + title('T9: Amp error 128 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT9) + ylabel('Error'); + xlabel('Frequency'); + title('T9: Phase error 128 cycles'); + +end + +Test9_AmpError=max(max(abs(Amp_error_totT9))); +Test9_PhaseError=max(max(abs(Phase_error_totT9))); + +try + assert( Test9_AmpError < Amp_tolerance, 'Test9 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test9_PhaseError < Phase_tolerance, 'Test9 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end +%% test 10 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 32 +disp('Test 10 - fixed cycles of 32 with dc offset'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 32; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; +DCoffset = 100; +DCoffsetinj = 400; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj );'); + Amp_error_totT10(iFreq) = mean(Amp_error); + Phase_error_totT10(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT10) + ylabel('Error'); + xlabel('Frequency'); + title('T10: Amp error 32 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT10) + ylabel('Error'); + xlabel('Frequency'); + title('T10: Phase error 32 cycles'); + +end + +Test10_AmpError=max(max(abs(Amp_error_totT10))); +Test10_PhaseError=max(max(abs(Phase_error_totT10))); + +try + assert( Test10_AmpError < Amp_tolerance, 'Test10 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test10_PhaseError < Phase_tolerance, 'Test10 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +%% test 11 + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 3 +disp('Test 11 - fixed cycles of 3 with dc offset'); + +Fc = [5:5:150 200:200:2000]; +FreqNum = size(Fc,2); + +Cycles = 3; +T=(1./Fc); %Period in s +InjTime=(T.*Cycles); + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; +DCoffset = 100; +DCoffsetinj = 400; + +for iFreq = 1:FreqNum + + evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj );'); + Amp_error_totT11(iFreq) = mean(Amp_error); + Phase_error_totT11(iFreq) = mean(Phase_error); + +end + +%% +if plotflag + figure + + subplot(2,1,1) + plot(Fc,Amp_error_totT11) + ylabel('Error'); + xlabel('Frequency'); + title('T11: Amp error 3 cycles'); + + subplot(2,1,2) + plot(Fc,Phase_error_totT11) + ylabel('Error'); + xlabel('Frequency'); + title('T11: Phase error 3 cycles'); + +end + +Test11_AmpError=max(max(abs(Amp_error_totT11))); +Test11_PhaseError=max(max(abs(Phase_error_totT11))); + +try + assert( Test11_AmpError < Amp_tolerance, 'Test11 Amp error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +try + assert( Test11_PhaseError < Phase_tolerance, 'Test11 Phase error failed') +catch err + Failed = 1; + fprintf(2, '%s\n', getReport(err, 'extended')); +end + +% SNR TEST + + +%% + +figure +hold on +plot(Fc,Amp_error_totT5); +plot(Fc,Amp_error_totT7); +plot(Fc,Amp_error_totT8); +hold off +ylabel('Error'); +xlabel('Frequency'); +legend('32','64','128') + + + + +%% Check a given ExpSetup + +load('S3b_MF1_log.mat','ExpSetup'); + + + +clear Amp_error_tot Phase_error_tot + +%changing freq - fixed cycles of 3 +disp('Test 12 - ExpSetup'); + +Fc = ExpSetup.Freq'; +FreqNum = size(Fc,2); + + +InjTime=ExpSetup.MeasurementTime/1000'; + + +Amp_Inj = 9500; +Amp_Meas = 1200 ; +InjPhase=0; +MeasPhaseDiff=-30; +DCoffset = 100; +DCoffsetinj = 400; + +for iFreq = 1:FreqNum + + % evalc('[Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj );'); + [Amp_error, Phase_error] = check_acc( Fc(iFreq),InjTime(iFreq),Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj ); + + + Amp_error_totT12(iFreq) = mean(Amp_error); + Phase_error_totT12(iFreq) = mean(Phase_error); + +end + +%% + +if plotflag + figure + + subplot(2,1,1) + plot(Fc,100*(Amp_error_totT12/Amp_Inj)) + ylabel('Error %'); + xlabel('Frequency'); + title('T12: Amp error 3 cycles'); + + subplot(2,1,2) + + ylabel('Error'); + xlabel('Frequency'); + title('T12: Phase error 3 cycles'); + +end +%% + +figure +hold on +plot(Fc/1000,100*(Amp_error_totT12/Amp_Inj)) +plot(Fc/1000,100*(Phase_error_totT12/MeasPhaseDiff)) + +ylabel('Error %'); +xlabel('Frequency (kHz)'); + +legend('Amplitude','Phase'); + + + +%% + + + + + + + + + + + + + + + + + + + + diff --git a/tests/filters/test_decimate.m b/tests/filters/test_decimate.m new file mode 100644 index 0000000..b3e8a05 --- /dev/null +++ b/tests/filters/test_decimate.m @@ -0,0 +1,55 @@ +Fcur = 1000; +FreqNum = size(Fcur,2); + +% Cycles = 2*6000; +% T=(1./Fcur); %Period in s +% InjTime=(T.*Cycles); + +InjTime=2; + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=0; +DCoffset = 0; +DCoffsetinj = 0; + +Fs=100000; + + +[Amp_error1, Phase_error1,V1,Vd1,Filt1,tr1] = check_acc( Fcur,InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj,[],Fs); +%% +dec_factor=2; + +Fs2=Fs/dec_factor; + +[Amp_error2, Phase_error2,V2,Vd2,Filt2,tr2] = check_acc( Fcur,InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj,[],Fs2); +%% +[Amp_error3, Phase_error3,V3,Vd3,Filt3,tr3] = check_acc( Fcur,InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj,[],Fs,dec_factor); + + + +%% + +t1 = (0:1:length(Vd1)-1)/Fs; +t2 = (0:1:length(Vd2)-1)/Fs2; +t3 = (0:1:length(Vd3)-1)/Fs2; + +figure; +hold on +plot(t1,Vd1); +plot(t2,Vd2); +plot(t3,Vd3); +hold off +legend('Original','Ideal downsample','Decimated'); + +figure; +hold on +plot(t1,V1); +plot(t2,V2); +plot(t3,V3); +hold off +legend('Original','Ideal downsample','Decimated'); +xlim([1, 1.01]) + diff --git a/tests/filters/test_detrend.m b/tests/filters/test_detrend.m new file mode 100644 index 0000000..9e561fb --- /dev/null +++ b/tests/filters/test_detrend.m @@ -0,0 +1,106 @@ + +Fc=20; + +Fcur = 15; +FreqNum = size(Fcur,2); + + +NumInj= 5; + +Cycles = 128; +T=(1./Fcur); %Period in s +InjTime=(T.*Cycles)*NumInj; +InjSamples=(T.*Cycles)*Fs; + +% InjTime=10; + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; + +inj = [ 2 4]; + +DCoffset=10+10*(1:NumInj); +DCoffsetInj=100 + 50*(1:NumInj); + + +Fs=16384; +BW=100; +chn=5; + +Padsec=4; + +%% Create ideal values, and voltages + +AmpActual= repmat(Amp_Meas,chn,1); +AmpActual(inj)=Amp_Inj; + +a=MeasPhaseDiff; + +%force normal valus of deg +a = mod(a,360); % [0 2pi) + +% shift +jj = a > 180; +a(jj) = a(jj) - 360; +jj = a < 0 - 180; +a(jj) = a(jj) + 360; + +MeasPhaseDiff_corr = a; + +MeasPhase = InjPhase + MeasPhaseDiff_corr; %deg +PhaseActual= repmat(MeasPhaseDiff_corr,chn,1); +PhaseActual(inj)=0; + +Totaltime=ceil(InjTime + 2*Padsec); +t = 0:1/Fs:Totaltime-1/Fs; +%make sin wave + +v_m= Amp_Meas*sin(2*pi*Fc*t+(pi*MeasPhase/180)); + +% change amplitude +v_i=Amp_Inj*sin(2*pi*Fc*t+(pi*InjPhase/180)); + +%% +V=repmat(v_m,chn,1); + +V(inj(1),:) = v_i; +V(inj(2),:) = v_i; + +V=V'; + +%pad with a second of data either side, so the hilbert is more realistic +datastart = round(Padsec*Fs); +dataend = round((Padsec+InjTime)*Fs); + +V(1:datastart,:)=V(1:datastart,:)*0.1; +V(dataend:end,:)=V(dataend:end,:)*0.1; + +InjectionWindows(:,1) = (datastart + ceil(InjSamples * (0:NumInj-1)))'; +InjectionWindows(:,2) = [InjectionWindows(2:end,1) ; dataend]; + + +DCoffsetVec = ones(chn,1)*DCoffset; +DCoffsetVec(inj(1),:) = DCoffsetInj; +DCoffsetVec(inj(2),:) = DCoffsetInj; + +for iInj = 1:NumInj + +V(InjectionWindows(iInj,1):InjectionWindows(iInj,2),:)=bsxfun(@plus,V(InjectionWindows(iInj,1):InjectionWindows(iInj,2),:),DCoffsetVec(:,iInj)'); + +end + + +%% + +Voff =mean(V(datastart:dataend,:)); + +plot(bsxfun(@minus,V,Voff)) + + + + + + + diff --git a/tests/filters/test_filters.m b/tests/filters/test_filters.m new file mode 100644 index 0000000..5a6b8c7 --- /dev/null +++ b/tests/filters/test_filters.m @@ -0,0 +1,209 @@ +%% + +% Fs=16384; +Fs=100000; + +Fc_all = 5:100:550; +FreqNum = length(Fc_all); + +BWtarget =100; + +Fstopomin=0.5; +Fcomin = 1.5; + + +for iFreq = 1: FreqNum + + + Fc=Fc_all(iFreq); + + disp(Fc); + + if Fc == 25; + disp('stop'); + end + + %% old filter + + BWold = BWtarget; + Fco(2) = Fc+(BWold/2); + Fco(1) = max([(Fc -(BWold/2)), Fstopomin]); + [B,A] = butter(3,Fco./(Fs/2)); + + stold(iFreq)= isstable(B,A); + imp_lenold(iFreq)=impzlength(B,A,0.001); + imp_len_msold(iFreq) = (imp_lenold(iFreq) / Fs)*(1000); + len_cyc = (1/Fc)*1000; + imp_len_cycold(iFreq)=imp_len_msold(iFreq)/len_cyc; + + [Hold,W] =freqz(B,A,Fco(1):0.2:Fco(2),Fs); + fc_idx = find (W > Fc,1); + Gainold(iFreq) =abs(Hold(fc_idx)); + + + %% New Bandpass + + BWbp = BWtarget; + + StopBandDiffbp = 350; + + Fpass1bp = Fc-BWbp/2; + + Astop1 = 60; + Astop2= 60; + Apass = 0.5; + + + if ~(Fpass1bp > 0) + Fpass1bp = max([Fc/2 Fcomin]); + Astop1 = 15; + end + + Fstop1bp = Fpass1bp -StopBandDiffbp; + Fstop1bp = max([Fstop1bp Fstopomin]); + + Fpass2bp = Fc+BWbp/2; + Fstop2bp = Fpass2bp + StopBandDiffbp; + + + + hbp = fdesign.bandpass('fst1,fp1,fp2,fst2,ast1,ap,ast2', Fstop1bp, Fpass1bp, ... + Fpass2bp, Fstop2bp, Astop1, Apass, Astop2, Fs); + + dbp = design(hbp, 'butter', 'MatchExactly', 'passband', 'SOSScaleNorm', 'Linf'); + + % dbp = designfilt('bandpassiir', ... % Response type + % 'FilterOrder', 20, ... + % 'HalfPowerFrequency1',Fpass1bp, ... % Frequency constraints + % 'HalfPowerFrequency2',Fpass2bp, ... + % 'SampleRate',Fs) ; % Sample rate + % + + + + stbp(iFreq)= isstable(dbp); + imp_lenbp(iFreq)=impzlength(dbp,0.001); + imp_len_msbp(iFreq) = (imp_lenbp(iFreq) / Fs)*(1000); + len_cyc = (1/Fc)*1000; + imp_len_cycbp(iFreq)=imp_len_msbp(iFreq)/len_cyc; + + [Hbp,W] =freqz(dbp,Fstop1bp:0.2:Fstop2bp,Fs); + fc_idx = find (W > Fc,1); + Gainbp(iFreq) =abs(Hbp(fc_idx)); + + + %% New Cascaded filter + + BWc = BWtarget; + + StopBandDiffc = 350; + + Fpass1c = Fc-BWc/2; + + Astop1 = 60; + Astop2= 60; + Apass = 0.5; + + if ~(Fpass1c > 0) + Fpass1c = max([Fc/2 Fcomin]); + Astop1 =30; + end + + Fstop1c = Fpass1c -StopBandDiffc; + Fstop1c = max([Fstop1c Fstopomin]); + + Fpass2c = Fc+BWc/2; + Fstop2c = Fpass2c + StopBandDiffc; + + + + % hchp = fdesign.highpass('fst,fp,ast,ap', Fstop1c, Fpass1c, Astop1, Apass, Fs); + % + % dchp = design(hchp, 'butter', 'MatchExactly', 'passband', 'SOSScaleNorm', 'Linf'); + + % hclp = fdesign.lowpass('fp,fst,ap,ast', Fpass2c, Fstop2c, Apass, Astop2, Fs); + % + % dclp = design(hclp, 'butter', 'MatchExactly', 'passband', 'SOSScaleNorm', 'Linf'); + + dclp = designfilt('lowpassiir', ... % Response type + 'FilterOrder', 10, ... + 'HalfPowerFrequency',Fpass2c, ... % Frequency constraints + 'SampleRate',Fs) ; % Sample rate + + + dchp = designfilt('highpassiir', ... % Response type + 'FilterOrder', 10, ... + 'HalfPowerFrequency',Fpass1c, ... % Frequency constraints + 'SampleRate',Fs) ; % Sample rate + + + stc(iFreq)= all([isstable(dchp) isstable(dclp)]); + imp_lenclp(iFreq)=impzlength(dclp,0.001) ; + imp_lenchp(iFreq) = impzlength(dchp,0.001); + imp_len_msclp(iFreq) = (imp_lenclp(iFreq) / Fs)*(1000); + imp_len_mschp(iFreq) = (imp_lenchp(iFreq) / Fs)*(1000); + len_cyc = (1/Fc)*1000; + imp_len_cycclp(iFreq)=imp_len_msclp(iFreq)/len_cyc; + imp_len_cycchp(iFreq)=imp_len_mschp(iFreq)/len_cyc; + [Hclp,W] =freqz(dclp,Fstop1c:0.2:Fstop2c,Fs); + [Hchp,W] =freqz(dchp,Fstop1c:0.2:Fstop2c,Fs); + fc_idx = find (W > Fc,1); + Gainclp(iFreq) = abs(Hclp(fc_idx)) ; + + Gainchp(iFreq) = abs(Hchp(fc_idx)); + + %% Freq Info + +end + + +%% + +figure; +hold on +plot(Fc_all,10*log10(Gainold)); +plot(Fc_all,10*log10(Gainbp)); +plot(Fc_all,10*log10(Gainclp)); +plot(Fc_all,10*log10(Gainchp)); +hold off +legend('Old','Bandp','Casclp','Caschp') +% legend('Bandp','Casclp','Caschp') +title('gain') +xlabel('Freq') +ylabel('Gain db'); + + +figure; +hold on +plot(Fc_all,imp_len_msold); +plot(Fc_all,imp_len_msbp); +plot(Fc_all,imp_len_msclp); +plot(Fc_all,imp_len_mschp); + +hold off +legend('Old','Bandp','Casclp','Caschp') +% legend('Bandp','Casclp','Caschp') +title('timedelay ms') +ylabel('Settling time ms') +xlabel('Freq') + + + +figure; +hold on +plot(Fc_all,imp_len_cycold); +plot(Fc_all,imp_len_cycbp); +plot(Fc_all,imp_len_cycclp); +plot(Fc_all,imp_len_cycchp); + +hold off +legend('Old','Bandp','Casclp','Caschp') +% legend('Bandp','Casclp','Caschp') +title('timedelay cycles') +ylabel('Settling time cycles') +xlabel('Freq') + + + + + diff --git a/tests/filters/test_singlecase.m b/tests/filters/test_singlecase.m new file mode 100644 index 0000000..1c8f25a --- /dev/null +++ b/tests/filters/test_singlecase.m @@ -0,0 +1,90 @@ +Fcur = 5; +FreqNum = size(Fcur,2); + +Cycles = 64; +T=(1./Fcur); %Period in s +InjTime=(T.*Cycles); + +% InjTime=10; + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; +DCoffset = 0; +DCoffsetinj = 0; + +% Fs=100000; +Fs=16384; + +[Amp_error1, Phase_error1,V1,Vd1,Filt1,tr1] = check_acc( Fcur,InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj,[],Fs); + +trim_demod=200; + + +t1 = (0:1:length(V1)-1)/Fs; + + +%% + +Fcur = 5; +FreqNum = size(Fcur,2); + +Cycles = 32; +T=(1./Fcur); %Period in s +InjTime=(T.*Cycles); + +% InjTime=10; + + +Amp_Inj = 500; +Amp_Meas = 150; +InjPhase=0; +MeasPhaseDiff=-30; +DCoffset = 100; +DCoffsetinj = 500; + +% Fs2=100000; +Fs2=16384; + + +[Amp_error2, Phase_error2,V2,Vd2,Filt2,tr2] = check_acc( Fcur,InjTime,Amp_Inj,Amp_Meas,InjPhase,MeasPhaseDiff,DCoffset,DCoffsetinj,[],Fs2); + + +t2 = (0:1:length(V2)-1)/Fs2; + +% + +%% + +fvtool(Filt1,Filt2) +%% +[H1,F1]=freqz(Filt1,Fcur-500:Fcur+500,Fs); +[H2,F2]=freqz(Filt2,Fcur-500:Fcur+500,Fs2); + +figure; +hold on +plot(F1,10*log10(abs(H1))); +plot(F2,10*log10(abs(H2))); +hold off + +%% + +figure +hold on +plot(t1,V1) +plot(t2,V2) +% plot(V8) +hold off + + +figure +hold on +% plot(Vd1(tr1:end-tr1)) +plot(t1,Vd1) +plot(t2,Vd2) +% plot(Vd2(tr2:end-tr2)) +% plot(Vd8) +hold off +% ylim(Amp_Inj + [-1 +1]) diff --git a/tests/readme.md b/tests/readme.md new file mode 100644 index 0000000..8f5515a --- /dev/null +++ b/tests/readme.md @@ -0,0 +1,6 @@ +# Functional Tests + +`Test_ScouseTom.m` tests require the files from [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4783547.svg)](https://doi.org/10.5281/zenodo.4783547) to be placed in the data directory as the EEG system files get very big very quickly. These tests only check that the code is functional and are by no means unit tests. + +### Filter tests +Checks to characterise frequency dependence of filtering method. Historical curiosity more than anything else.