% Eric LaForest, April 2010
% laforest AT eecg.utoronto.ca

function [results] = project3(time)
% Wiener-Bose nonparametric model of FHN parametric model

    function output = driven_FHN(t, input, stimulus, time_step) 
        V = input(1);
        W = input(2);

        % Pre-set FitzHugh-Nagumo System driven by time-series
        index = round((t / time_step) + 1);
        Vdot = V - (V^3 / 3) - W + stimulus(index);
        Wdot = 0.08 * (V + 0.7 - 0.8 * W);
        
        output = [Vdot; Wdot;];
    end

    function  [T,R] = solve_driven_FHN(time, initial_conditions, stimulus, time_step)
        this_FHN = @(t,input) driven_FHN(t, input, stimulus, time_step);
        [T,R] = ode15s(this_FHN, time, initial_conditions);
    end

    function [GWN] = generate_GWN(time, amplitude)
        GWN_raw = rand(size(time)) - 0.5; % for a zero mean
        GWN = GWN_raw * amplitude * 2;
    end

    % Based on:
    % http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect9.html
    % http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml
    function [bandwidth, magnitude] = measure_bandwidth(time, time_series, time_step)
        padded_length = 2^nextpow2(length(time));
        spectrum = fft(time_series, padded_length) / length(time);
        spectrum = spectrum(1:length(spectrum)/2+1);
        magnitude = abs(spectrum);
        sampling_frequency = 1/time_step;
        bandwidth = sampling_frequency/2 * linspace(0,1,padded_length/2+1);
    end

    function plot_bandwidth(time)
        time_step = time(2) - time(1);
        GWN = generate_GWN(time, 2);
        [T,R] = solve_driven_FHN(time, [0,0], GWN, time_step);
        [bandwidth, magnitude] = measure_bandwidth(T, R, time_step); 
        plot(bandwidth, magnitude);
    end


    % Not used here, but if you need it: a memoized recursive Nth-order Wiener-Bose filter
    function output_series = discrete_WB_linear_filter(order, alpha, time_series)
        len = length(time_series);
        output_series = zeros(1,len);
        
        discrete_WB_linear_filter_1_output = repmat(NaN,1,len);
        discrete_WB_linear_filter_n_inner_output = repmat(NaN,order,len);
        
        filter_previous_1 = repmat(NaN,1,len);
        filter_previous_n = repmat(NaN,order,len);
        filter_current_lower_order_n = repmat(NaN,order,len);
        filter_previous_lower_order_n = repmat(NaN,order,len);
        
        function output = discrete_WB_linear_filter_n(order, alpha, time_series, index)
            if (index == 0)
                output = 0;
                return;
            end

            if (order < 1)
                error('Filter order must be >= 1');
            end

            function output = discrete_WB_linear_filter_1(alpha, time_series, index)
                if (index == 0)
                    output = 0;
                    return;
                end

                if (~isnan(discrete_WB_linear_filter_1_output(index)))
                    output = discrete_WB_linear_filter_1_output(index);
                    return;
                end
                
                if (isnan(filter_previous_1(index)))
                    filter_previous_1(index) = discrete_WB_linear_filter_1(alpha, time_series, index - 1);
                end

                filter_previous = filter_previous_1(index);
                current_point = time_series(index);
                discrete_WB_linear_filter_1_output(index) = (sqrt(alpha) * filter_previous) + (sqrt(1 - alpha) * current_point);
                output = discrete_WB_linear_filter_1_output(index);
            end

            function output = discrete_WB_linear_filter_n_inner(order, alpha, time_series, index)
                if (index == 0)
                    output = 0;
                    return;
                end

                if (~isnan(discrete_WB_linear_filter_n_inner_output(order, index)))
                    output = discrete_WB_linear_filter_n_inner_output(order, index);
                    return;
                end
                
                if (order == 1)
                    output = discrete_WB_linear_filter_1(alpha, time_series, index);
                    return;
                end

                if (isnan(filter_previous_n(order, index)))
                    filter_previous_n(order, index) = discrete_WB_linear_filter_n_inner(order, alpha, time_series, index - 1);
                end
                filter_previous = filter_previous_n(order, index);

                if (isnan(filter_current_lower_order_n(order, index)))
                    filter_current_lower_order_n(order, index) = discrete_WB_linear_filter_n_inner(order - 1, alpha, time_series, index);
                end
                filter_current_lower_order = filter_current_lower_order_n(order, index);

                if (isnan(filter_previous_lower_order_n(order, index)))
                    filter_previous_lower_order_n(order, index) = discrete_WB_linear_filter_n_inner(order - 1, alpha, time_series, index - 1);
                end
                filter_previous_lower_order = filter_previous_lower_order_n(order, index);

                discrete_WB_linear_filter_n_inner_output(order, index) = sqrt(alpha) * (filter_previous + filter_current_lower_order) - filter_previous_lower_order;
                output = discrete_WB_linear_filter_n_inner_output(order, index);
            end

            output = discrete_WB_linear_filter_n_inner(order, alpha, time_series, index);
        end
    
        for index = 1:len
            output_series(index) = discrete_WB_linear_filter_n(order, alpha, time_series, index); 
        end
    end
%{
tic;
plot(discrete_WB_linear_filter(10, 0.2, [1,zeros(1,10000)]));
toc;
%}  

    function data = read_lysis_file(filename)
        fid = fopen(filename, 'r');
        % skip the header, as MATLAB can't digest it worth a @#$%
        for i = 1:7
            line = fgetl(fid);
        end
        data = cell2mat(textscan(fid, '%f%*s', 'delimiter', '\r\n'));
        fclose(fid);
    end

    function write_lysis_header(fid, filename, comment)
       fprintf(fid, ' %s\r\n', filename);
       fprintf(fid, ' %s\r\n', comment);
       % ***** THESE ARE HARDCODED ****
       fprintf(fid, 'TS\r\n');
       fprintf(fid, ' 0.0000000E+00\r\n');
       fprintf(fid, '   99.95000\r\n');
       fprintf(fid, '  5.0000001E-02\r\n');
       fprintf(fid, '        2000\r\n');        
    end

    function write_lysis_file(filename, input_filename, data)
       fid = fopen(filename, 'w');
       comment = sprintf('Output derived from %s', input_filename);
       write_lysis_header(fid, filename, comment);
       fprintf(fid, ' %f\r\n', data);
       fclose(fid);
    end

    function write_cos_file(filename, amplitude, frequency)
        ts = 0.05;
        time = 0:ts:99.95;
        data = amplitude * cos(time * frequency);
        comment = sprintf('%f * cos(t * %f)', amplitude, frequency);
        write_lysis_file(filename, comment, data);
    end

    function write_const_file(filename, amplitude)
        ts = 0.05;
        time = 0:ts:99.95;
        data = ones(size(time)) * amplitude;
        comment = sprintf('constant amplitude %f', amplitude);
        write_lysis_file(filename, comment, data);        
    end

ts = 0.05;
time = 0:ts:99.95;

%input_filename = 'const0';
%output_filename = strcat('output-', input_filename);

%stimulus = read_lysis_file(input_filename);
%stimulus = [ones(1,2)*8, zeros(1,2000)];
%stimulus = generate_GWN(time, 0.9);
%[T,R] = solve_driven_FHN(time, [-1.1994,-0.6243], stimulus, ts);
%plot(T,R)
%[bandwidth, magnitude] = measure_bandwidth(T, R, ts); 
%[bandwidth, magnitude] = measure_bandwidth(time, stimulus, ts); 
%plot(bandwidth, magnitude);

%write_lysis_file(output_filename, input_filename, R(:,1));

%{
write_cos_file('cos11', 1, 1);
write_cos_file('cos3333', 0.33, 0.33);
write_cos_file('cos331', 0.33, 1);
write_cos_file('cos145533', 1.455, 0.33);
write_cos_file('cos14551', 1.455, 1);

write_const_file('const33', 0.33);
write_const_file('const1', 1);
write_const_file('const1455', 1.455);
write_const_file('const1minus', -1);
write_const_file('const0', 0);
%}

%{
input_filename = 'gwn1-10.0';
%output_filename = strcat('output-', input_filename);
output_filename = 'output-gwn1-10.0';
system = read_lysis_file(input_filename);
model = read_lysis_file(output_filename);

plot(time, system, 'b', time, model, 'r');
xlabel('Time');
ylabel('Voltage');
title('Suprathreshold FHN System (GWN1 stimulus, +/-10V peak)');
legend('GWN1 Stimulus', 'Response');
%}

%{
input_filename = 'lsk1sup';
data = read_lysis_file(input_filename);

plot(data, 'b');
xlabel('points (512)');
ylabel('Magnitude');
title('First-Order Kernel of High-Power Wiener-Bose Model (Lee-Schetzen Estimate)');
%legend('GWN1 Stimulus', 'Response');
%}

%{
numel = 9;
input_filename = 'ltc2sub95-512';
data = read_lysis_file(input_filename);
k2 = zeros(numel, numel);
for i = 1:numel
    index = (i + ((i-1)*numel));
    range = index:(index+(numel-1));
    k2(1:numel,i) = data(range);
end
bar3(k2, 0.5, 'b');
grid on
xlabel('X Points');
ylabel('Y Points');
zlabel('Z Magnitude');
title('Coefficients of Second-Order Kernel of Low-Power Wiener-Bose Model');
%}

%{
figures = importdata('figures_list');
for i = 1:length(figures)
    figure = cell2mat(figures(i));
    handle = hgload(figure);
    print('-depsc2', handle, strcat(figure, '.eps'));
    close all;
end
%}

%{ 
input_filename = 'output-gwn1-10.0';
stimulus = read_lysis_file(input_filename);
[bandwidth, magnitude] = measure_bandwidth(time, stimulus, ts); 
semilogy(bandwidth, magnitude);
%xlim([0 2]);
%ylim([0 0.07]);
xlabel('Frequency (Hz)');
ylabel('Power');
title('Frequency Spectrum of Suprathreshold FHN System Response');
%}

end