% GJLFP, eLife 2026. Richa Sirmaur and Rishikesh Narayanan.
% Figure 5: Patterned inputs: waveform, fourier, wavelet

%%
fs = 40000;
freqs = [1 4 16 64];

% bandpass ranges
fc1_all = [0.5 3 10 60]; % cutoff setting for 1, 4, 16, 64 Hz respectively
fc2_all = [1.5 5 20 70]; % cutoff setting for 1, 4, 16, 64 Hz respectively

colIdx = 24; % electrode (≈97 µm)

time = 0:0.025:2000;
T = time(1:80000).';

for fidx = 1:length(freqs)

    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);

    % LOAD DATA (diff frequency folders)
    cd(sprintf('../chemsyn/ACTIVE/Pattern/nosod/Extra/%dHz', f0))
    mat_act = load('allMATS.mat');

    cd(sprintf('../chemsyn/PASSIVE/Pattern/nosod/Extra/%dHz', f0))
    mat_pass = load('allMATS.mat');

    Wn = [fc1 fc2] / (fs/2);
    [b,a] = butter(1, Wn, 'bandpass');
    MatA = mat_act.Mat3;
    MatP = mat_pass.Mat3;
    MatA(1:2001,:) = repmat(MatA(2002,:),2001,1); % remove initial (50 ms) transient
    MatP(1:2001,:) = repmat(MatP(2002,:),2001,1);
    dataA = filtfilt(b,a,MatA(:,colIdx));
    dataP = filtfilt(b,a,MatP(:,colIdx));

    % time window (10 cycles each) for example waveform 
    period = 1000 / f0;   % ms
    win    = 10 * period; % showing 10 cycles for waveform plots
    tStart = 500;
    tEnd   = tStart + win;
    idx = T >= tStart & T <= tEnd;

    figure('Name', sprintf('%d Hz', f0));

    % 1) WAVEFORM
    subplot(3,1,1)
    plot(T(idx), dataA(idx), 'r', 'LineWidth', 1.7); hold on;
    plot(T(idx), dataP(idx), 'b', 'LineWidth', 1.7);
    xlabel('Time (ms)');
    ylabel('Amplitude (nV)');
    title(sprintf('LFPs Waveforms (%d Hz)', f0));
    legend('Active','Passive');
    box off;

    % 2) FOURIER 
    L = length(dataA);
    YA = fft(dataA);
    YP = fft(dataP);
    P1A = abs(YA/L);
    P1P = abs(YP/L);
    P1A = P1A(1:L/2+1);
    P1P = P1P(1:L/2+1);
    P1A(2:end-1) = 2*P1A(2:end-1);
    P1P(2:end-1) = 2*P1P(2:end-1);
    f = fs*(0:(L/2))/L;

    subplot(3,1,2)
    plot(f, P1A, 'r', 'LineWidth', 2); hold on;
    plot(f, P1P, 'b', 'LineWidth', 2);
    set(gca,'XScale','log');

    xlim([1e-1 1e3]); % f0 == 1 % change frequency-dependent limits
    xlim([1e0 1e3]); % others
    xlabel('Frequency (Hz)');
    ylabel('Power (a.u)');
    title('Fourier Transform (%d Hz)', f0);
    legend('Active','Passive');
    box off;

    % 3) WAVELET 
    subplot(3,1,3)

    padA = [zeros(200000,1); dataA; zeros(200000,1)]; % zeropadding for active
    % padP = [zeros(200000,1); dataP; zeros(200000,1)];
    cwt(padA,'morse',fs);
    colormap(jet);
    colorbar;
    axis xy;
    title('Wavelet (Active- %d Hz)', f0);
    xlabel('Time');
    ylabel('Frequency');

end