#Post-process the time-domain output power to obtain the RF response
#
#INPUTS:
#box_filter_window - time window for the box filtering
#maxFreq - maximum frequency for the RF response

function rf_bandwidth_postprocess(box_filter_window,maxFreq){
    
    step_excitation = getnamed('STEP_1','delay');
    
    osc = getresult("OOSC_1","mode 1/signal");
    power=osc.getattribute(osc.getattribute);
    t = osc.getparameter("time");
    
    #plot(t,power);
    #setplot("y max",0.1);
    
    #Downsample.
    t = t(1:10:end);
    power = power(1:10:end);
    
    #filter out higher frequency oscillations
    dt = t(2)-t(1);
    box_filter_window = floor(box_filter_window/dt);
    power_filtered = power;
    for(i=1; i<=length(t); i=i+1){
        window_range = (i-floor(box_filter_window/2)):(i+floor(box_filter_window/2));
        window_range = window_range(find(window_range>=1));
        window_range = window_range(find(window_range<=length(t)));
        power_filtered(i) = mean(power(window_range));
    }
    power = power_filtered;
    for(i=1; i<=length(t); i=i+1){#repeat
        window_range = (i-floor(box_filter_window/2)):(i+floor(box_filter_window/2));
        window_range = window_range(find(window_range>=1));
        window_range = window_range(find(window_range<=length(t)));
        power_filtered(i) = mean(power(window_range));
    }
    power = power_filtered;
    for(i=1; i<=length(t); i=i+1){#repeat
        window_range = (i-floor(box_filter_window/2)):(i+floor(box_filter_window/2));
        window_range = window_range(find(window_range>=1));
        window_range = window_range(find(window_range<=length(t)));
        power_filtered(i) = mean(power(window_range));
    }
    
    #plot(t,power_filtered,"time [s]","power [W]","","linewidth=2");
    
    #Select only times after turning on small signal step excitation
    t_step = t(find(t>step_excitation));
    power_step = power_filtered(find(t>step_excitation));
    
    ## Calculate the 3dB bandwidth by using FFT to get the frequency response
    I=power_step;
    t=t_step;
    N= length(t);
    th = t;
    Nh = N;
    x=t;
    y=I;
    nuderiv;
    dIdt = dydx;
    
    #plotxy(t*1e12,dIdt/max(dIdt),t*1e12,I/max(I),"time (ps)","Normalized amplitude (a.u) ");
    
    # Interpolate the impulse response on a dense grid
    t_interp = th(1):0.1e-12:th(Nh);  #uniform time grid
    dIdt_interp = interp(dIdt, th,t_interp);
    #plotxy(t_interp*1e12,dIdt_interp/max(dIdt_interp),t*1e12,I/max(I),"time (ps)","Normalized amplitude (a.u) ");
    #legend("impulse response","step response");
    
    # take fft of the original impulse response to get the frequency response
    RF = fft(dIdt_interp,2,2^12);
    RF = 10*log10(abs(RF/RF(1)));
    w = fftw(t_interp - t_interp(1),2,2^12);
    freq = w/2/pi;
    RF = RF(find(freq<=maxFreq));
    freq = freq(find(freq<=maxFreq));

    output = struct;
    output.bandwidth = freq(find(RF,-3));
    output.rf_response = RF;
    output.freq = freq;
    return output;
}