% Scope:
% This script runs the Method Validation (Recorded Data) for the article
% G. Fischer, J. Haueisen, D. Baumgarten & M. Kofler
% "Converting spectral evoked-to-background ratio into time-domain
% signal-to-noise ratio – validation for high frequency oscillations"
% MethodsX 15 (2025) 103617; https://doi.org/10.1016/j.mex.2025.103617
%
% (c) Gerald Fischer 2025/09/22
% -------------------------------------------------------------------------
% subject 47e male; 187 cm; right median nerve stimulation
% preprocessed signal data (stim-spikes removed); data in structure
load StSignal_MED47e_C395.mat; 
fsample = StSignal.fsample;               % assign sample rate
ntrig=size(StSignal.nTrigger,2);          % get number of triggers
% do not use first and last trigger (may contain corrupted data)
StSignal.nTrigger=StSignal.nTrigger(2:ntrig-1);
Nsweep=[90:15:1005]; Nsweep(1,62)=1004;   % sweep count steps; max 1004
% -------------------------------------------------------------------------
% compute and expected SNR 
% -> load computed NFTA data structure for same recording (MethodsX 2023)
load StNFTA_MED_47e_C395.mat;
TWRdB = 10*log10(1/(0.01*3.96));    % 10 ms window in 3.96 Hz cycle
EBRdB = getEBRdB_HFO_01(StNFTA);    % compute EBRdB in subroutine
SNRe=TWRdB+EBRdB+10*log10(Nsweep);  % compute expected SNR
% -------------------------------------------------------------------------
% compute HFO using the zero-phase IIR filter
[b,a]=butter(2,[400 800]/4800,'bandpass');  % filter coefficients
filtdata=filtfilt(b,a,StSignal.data);       % zero-phase IIR-filter
% perform trial averaging and get SNR
[avgDataIIR] = trialAverage(filtdata',StSignal.nTrigger);
[SNRiir] = getSNR4MethX(avgDataIIR,Nsweep,fsample);
% -------------------------------------------------------------------------
% compute HFO using two Ricker kernels in series
[kRick] = getRickerKernel(585,fsample); % kernel for 585 Hz center freq.
kRick2=conv(kRick,kRick,'full');        % 2 kernels in series
filtdata=conv(StSignal.data,kRick2,'same');  % filter by convolution
% perform trial averaging and get SNR
[avgDataRic] = trialAverage(filtdata',StSignal.nTrigger);
[SNRric] = getSNR4MethX(avgDataRic,Nsweep,fsample);
% -------------------------------------------------------------------------
% check SNR
disp(['ZP-IIR vs expected: cc ' ...
    num2str(corr(SNRiir',SNRe')) ' rms ' num2str(rms(SNRiir-SNRe))]);
disp(['ZP-IIR vs 2xRicker: cc ' ...
    num2str(corr(SNRiir',SNRric')) ' rms ' num2str(rms(SNRiir-SNRric))]);
% -------------------------------------------------------------------------
% plot results
figure('units','centimeters','InnerPosition',[1 1 18.5 8]);
subplot(1,2,2);
time = avgDataIIR.time*1e3;    % time vector in [ms]
HFOiir=avgDataIIR.mtx(1004,:); % last IIR waveform
HFOric=avgDataRic.mtx(1004,:); % last double Ricker waveform
hold on;
plot(time,-HFOiir,'color',[0.4 0 .9],'LineWidth',1.5); % negativity up
plot(time,-HFOric,'color',[0.9 0 .4],'LineWidth',1.5); % negativity up
plot([15 15],[-0.5 0.5],'k--'); % response window left boarder
plot([25 25],[-0.5 0.5],'k--'); % response window right boarder
legend('zero-phase IIR','2\times Ricker','location','NorthEast');
yticks([-0.3 -0.2 -0.1 0 0.1 0.2 0.3]);
yticklabels({'0.3','0.2','0.1','0.0','-0.1','-0.2','-0.3'})
xlim([0 50]);
ylim([-0.35 0.35]);
ylabel('amplitude [\muV]');
xlabel('time [ms]');
% plot SNR
subplot(1,2,1);
hold on;
plot(Nsweep,SNRiir(1,:),'color',[0.4 0 .9],'LineWidth',1.5);
plot(Nsweep,SNRric(1,:),'color',[0.9 0 .4],'LineWidth',1.5);
plot(Nsweep,SNRe,'k--','LineWidth',1);
legend('zero-phase IIR','2\times Ricker','expected','location','NorthWest');
ylim([0 12]);
xlim([0 1005]);
xlabel('sweep count');
ylabel('SNR [dB]');
grid on;

