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

%% this is  from default models- chemical synapse. same for no sodium models- chemical synapses
cd '/chemsyn/active/LFP'
load('allMATS.mat');

fs = 40000;
fc1 = 0.5;
fc2 = 300;

Wn = [fc1 fc2] / (fs/2);
n = 2;
[b,a] = butter(n, Wn, 'bandpass');

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

locations = {'proximal','intermediate','distal'};

locMap = {
    'Mat5', 24;   % intermediate; 97.15um
    'Mat3', 24;   % proximal ; 55 um
    'Mat7', 38;    % distal;152.19   
};

matStruct = struct();
matStruct.Mat5 = Mat5;
matStruct.Mat3 = Mat3;
matStruct.Mat7 = Mat7;

for l = 1:numel(locations)

    locName = locations{l};

    matName = locMap{l,1};
    colIdx  = locMap{l,2};

    dataRaw = matStruct.(matName)(:, colIdx);

    data = filtfilt(b,a,dataRaw);

    % 1) WAVEFORM
    figure;
    subplot(3,1,1)
    plot(T, data, 'm', 'LineWidth', 1.7);
    xlim([400 1400]);
    ylim([-2000 1000]);
    xlabel('Time (ms)');
    ylabel('Amplitude (nV)');
    title([locName ' - Waveform']);
    box off;

    % 2) FOURIER
    Y = fft(data);  % FFT by column of Filteredfielddata array
    L = length(data);
    P2 = abs(Y/L); % 2-sided PSD, power
    P1 = P2(1:L/2+1); % 1-sided
    P1(2:end-1) = 2*P1(2:end-1); % scale, don't double DC and Fmax components

    f = fs*(0:(L/2))/L; % compute frequency

    subplot(3,1,2)
    plot(f, P1, 'LineWidth', 2);
    xlim([0 1000]);
    ylim([0 100]);
    xlabel('Frequency (Hz)');
    ylabel('Power (a.u.)');
    title([locName ' - Fourier Transform']);
    box off;

    
    % 3) WAVELET
    subplot(3,1,3)
    disp ('zero padding')
    padData = [zeros(40001,1); data; zeros(40001,1)];

    cwt(padData, 'morse', fs);
    colormap(jet);
    colorbar;
    axis xy; % Flip the y-axis to have lower scales at the bottom

    title([locName ' - Wavelet']);
    xlabel('Time');
    ylabel('Frequency');
    xlim([1.5 2.1])
    minValue=0;
    maxValue=1200;
    caxis([minValue, maxValue]);

end

%% same for gap junctions..